妙妙题。
首先可以有一个 \(O(kn^2)\) 的 dp，但是显然不行。
但是，发现其中的大多数转移都浪费在自环上了，所以考虑不要这个东西。
这个 dp 一共有三种转移：
-
左右端点一起向内移动一格；
-
左端点或右端点单独移动；
-
左右端点都不动。
所以考虑加一维 \(k\) 表示走了 \(k\) 次转移 1，这样可以计算出转移 2 的次数。
最后对于转移 1 走了 \(i\) 次，转移 2 走了 \(j\) 次的方案数一起计算贡献：
\[g_{i,j}=\sum\limits_{x=0}^{m-i-j}26^x\sum\limits_{y=0}^{m-i-j-x}24^y\times 25^{m-i-j-x}\times \binom{k+j-1}{j-1}\times \binom{m-j-x-1}{i-1}
\]
如果使用生成函数的话是可以解决的，但是想想细节太多，不想写。
不过值得注意的是，\(g_{i,j}=25g_{i,j-1}-24g_{i-1,j}\)。
所以换一种思路，这种东西也是可以矩阵乘法得到的，复杂度为 \(O(n^4\log k)\)，过不去。
但是发现其实可以把这些所有的 \(i,j\) 放在同一个矩阵里面转移，所以复杂度变成 \(O(n^3\log k)\)，可以过去。
注意有个细节，就是如果最后长度为奇数，那么最后一步转移不能是转移 3，把贡献减去即可。
代码
#include<bits/stdc++.h>
using namespace std;
using ll=long long;
const int N=2e2+10,M=N*1.5,mod=1e4+7;
int n,m,f[N][N][N];
char a[N];
int cnt,p[N],q[N];
struct matrix{
int a[M][M];
matrix(){
memset(a,0,sizeof a);
}
matrix operator * (const matrix &x)const{
matrix b;
for(int k=0;k<cnt;k++)
for(int i=0;i<=k;i++)
for(int j=k;j<cnt;j++)
b.a[i][j]=(b.a[i][j]+a[i][k]*x.a[k][j])%mod;
return b;
}
}base;
struct vec{
int a[M];
vec(){
memset(a,0,sizeof a);
}
vec operator * (const matrix &x)const{
vec b;
for(int i=0;i<cnt;i++)
for(int j=i;j<cnt;j++)
b.a[j]=(b.a[j]+a[i]*x.a[i][j])%mod;
return b;
}
}ans;
int main(){
freopen(".in","r",stdin);
//freopen(".out","w",stdout);
scanf("%s%d",a+1,&m),n=strlen(a+1);
bool flag=(m+n)&1;
m=(m+n+1)/2;
f[1][n][0]=1;
for(int l=1;l<=n;l++){
for(int r=n;r>=l;r--){
for(int i=0;i<=n;i++){
if(!f[l][r][i])continue;
if(a[l]==a[r]){
(f[l+1][r-1][i+1]+=f[l][r][i])%=mod;
}
else{
(f[l+1][r][i]+=f[l][r][i])%=mod;
(f[l][r-1][i]+=f[l][r][i])%=mod;
}
}
}
}
for(int i=0;i<=n;i++)p[i]=cnt++;
for(int i=(n+1)/2;i>=0;i--)q[i]=cnt++;
for(int l=1;l<=n+1;l++){
for(int r=0;r<l;r++){
for(int i=0;i<=n;i++){
if(!f[l][r][i])continue;
int j=(l-1)+(n-r)-i*2;
// fprintf(stderr,"%d %d %d\n",i,j,f[l][r][i]);
base.a[p[j]][q[i]]=(base.a[p[j]][q[i]]+f[l][r][i])%mod;
}
}
}
for(int i=1;i<=n;i++)base.a[p[i]][p[i]]=24;
for(int i=0;i<n;i++)base.a[p[i]][p[i+1]]=1;
for(int i=1;i<=(n+1)/2;i++)base.a[q[i]][q[i]]=25,base.a[q[i]][q[i-1]]=1;
base.a[q[0]][q[0]]=26;
ans.a[p[0]]=1;
for(int y=m+1;y;base=base*base,y>>=1)if(y&1)ans=ans*base;
int res=ans.a[q[0]];
if(flag){
base=matrix(),ans=vec();
for(int l=0;l<=n+1;l++){
for(int r=0;r<l;r++){
for(int i=1;i<=n;i++){
if(!f[l][r][i]||l!=r+1)continue;
int j=(l-1)+(n-r)-i*2;
base.a[p[j]][q[i]]=(base.a[p[j]][q[i]]+f[l][r][i])%mod;
}
}
}
for(int i=1;i<=n;i++)base.a[p[i]][p[i]]=24;
for(int i=0;i<n;i++)base.a[p[i]][p[i+1]]=1;
for(int i=1;i<=(n+1)/2;i++)base.a[q[i]][q[i]]=25,base.a[q[i]][q[i-1]]=1;
base.a[q[0]][q[0]]=26;
ans.a[p[0]]=1;
for(int y=m;y;base=base*base,y>>=1)if(y&1)ans=ans*base;
(res+=mod-ans.a[q[1]])%=mod;
}
cout<<res;
return 0;
}