#include <bits/stdc++.h>
using namespace std;
const double PI = acos(0) * 2;
const double EPS = 1e-8;
const long long MOD = 51123987;
const int MAXN = 2e6 + 5;
const int oo = 1e9;
const double foo = 1e30;
template <class T>
int getbit(T s, int i) {
  return (s >> i) & 1;
}
template <class T>
T onbit(T s, int i) {
  return s | (T(1) << i);
}
template <class T>
T offbit(T s, int i) {
  return s & (~(T(1) << i));
}
template <class T>
int cntbit(T s) {
  return __builtin_popcounll(s);
}
inline void addmod(long long& a, long long val, int p = MOD) {
  if ((a = (a + val)) >= p) a -= p;
}
inline void submod(long long& a, long long val, int p = MOD) {
  if ((a = (a - val)) < 0) a += p;
}
inline int mult(int a, int b, int p = MOD) { return (long long)a * b % p; }
const int MAXC = 26;
struct palindromicTree {
  int s[MAXN], fail[MAXN], len[MAXN], num[MAXN], cnt[MAXN];
  int n, p, last;
  vector<pair<int, int> > nxt[MAXN];
  int newnode(int l) {
    nxt[p].clear();
    cnt[p] = fail[p] = num[p] = 0;
    len[p] = l;
    return p++;
  }
  void init() {
    p = 0;
    newnode(0);
    newnode(-1);
    n = last = 0;
    s[n] = -1;
    fail[0] = 1;
  }
  int getfail(int x) {
    while (s[n - len[x] - 1] != s[n]) x = fail[x];
    return x;
  }
  int find(int cur, int c) {
    for (int i = 0; i < ((int)(nxt[cur]).size()); i++) {
      if (nxt[cur][i].first == c) return nxt[cur][i].second;
    }
    return 0;
  }
  int add(char c) {
    c -= 'a';
    s[++n] = c;
    int cur = getfail(last);
    if (!find(cur, c)) {
      int now = newnode(len[cur] + 2);
      fail[now] = find(getfail(fail[cur]), c);
      nxt[cur].push_back(make_pair(c, now));
      num[now] = num[fail[now]] + 1;
    }
    last = find(cur, c);
    cnt[last]++;
    return num[last];
  }
  void count() {
    for (int i = p - 1; p >= 0; p--) {
      cnt[fail[i]] += cnt[i];
    }
  }
} PA;
int n;
char s[2000100];
long long finish[2000100], start[2000100];
int main() {
  scanf("%d", &n);
  scanf("%s", s);
  PA.init();
  for (int i = 0; i < n; i++) {
    finish[i] = PA.add(s[i]);
  }
  PA.init();
  for (int i = 0; i < n; i++) {
    start[n - i - 1] = PA.add(s[n - i - 1]);
  }
  long long ans = 0;
  for (int i = 0; i < n; i++) {
    ans = (ans + start[i]) % MOD;
  }
  ans = (ans * ((ans - 1ll + MOD) % MOD)) / 2ll;
  ans %= MOD;
  long long tot = 0;
  for (int i = 0; i < n; i++) {
    long long tmp = (tot * start[i]) % MOD;
    ans = (ans - tmp + MOD) % MOD;
    tot += finish[i];
  }
  cout << ans % MOD;
}
