#include <bits/stdc++.h>
using namespace std;
const int maxn = 2e6 + 5;
const int mod = 51123987;
char s[maxn];
struct Pam {
  vector<pair<int, int> > next[maxn];
  int fail[maxn];
  int len[maxn];
  int num[maxn];
  int S[maxn];
  long long a[maxn], b[maxn];
  int last, n, p;
  int newNode(int l) {
    next[p].clear();
    len[p] = l;
    return p++;
  }
  void init() {
    n = last = p = 0;
    newNode(0);
    newNode(-1);
    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 u, int c) {
    vector<pair<int, int> >& x = next[u];
    int sz = x.size();
    for (int i = 0; i < sz; ++i) {
      if (x[i].first == c) return x[i].second;
    }
    return 0;
  }
  int add(int c) {
    S[++n] = c;
    int cur = getFail(last);
    int x = find(cur, c);
    if (!x) {
      int now = newNode(len[cur] + 2);
      x = now;
      fail[now] = find(getFail(fail[cur]), c);
      next[cur].push_back(make_pair(c, now));
      num[now] = num[fail[now]] + 1;
    }
    last = x;
    return num[last];
  }
  void solve(int len) {
    memset(a, 0, sizeof(a));
    memset(b, 0, sizeof(b));
    init();
    for (int i = 0; i < len; i++) {
      a[i] = add(s[i] - 'a');
    }
    init();
    for (int i = len - 1; i >= 0; i--) {
      b[i] = (b[i + 1] + add(s[i] - 'a')) % mod;
    }
    long long ans = (b[0] * (b[0] - 1) / 2) % mod;
    for (int i = 0; i < len; i++) {
      ans = ((ans - (a[i] * b[i + 1]) + mod) % mod) % mod;
    }
    printf("%lld\n", ans);
  }
} pam;
int main() {
  int n;
  while (~scanf("%d", &n)) {
    scanf("%s", s);
    pam.solve(n);
  }
  return 0;
}
