**Warning: this post is very technical. You may find it very hard to understand, or boring.**

One of my friend asked me about this problem.

The description of this problem is simple: given a string pattern, find the expected position of first occurrence of such pattern in random strings.

The problem looks like dynamic programming, because there seems exists some recursive relations between its subproblems.

Let’s first assume there exists a random string generator, which generate random strings letter by letter. Left F[n] be the **number of expected letters to be generated before the pattern occur, given that most recently generated n letters matches first n letters of the pattern**. For example, assume the pattern is *ABCDE*, and the random string generator generates *….ABC*, then F[3] would be the expected number of letters needed before *ABCDE* is generated.

We can consider “N-prefix” of the pattern as “state” for DP(described by N), and F[0] would be the answer we need: expected number of letters needed before pattern is observed, given no prefix match(random string generator hasn’t generated anything yet).

Now let’s see how one state transfer into another state. Given that most recent generated N letters match N-prefix of patter, if one more letters is generated, how would be the suffix-prefix match goes? In the best case, it will become N+1, if the newly generated letter matches the N+1th letter in the pattern. In other cases, the state would transfer to some number no larger than N, and that number can be found easily: find the largest prefix-suffix match between generated string and pattern.

Now the recursive relationship between states becomes clear:** F[i] = 1+\sum_{c=’A’}^{‘Z’}F[next(i,c)]/26**. Base case would be **F[pattern.length]=0**.

In the original problem, number of letters is a variable, rather than fixed 26. **next(i,c)** is state transition function.

Now comes another problem: this DP cannot be solved because states graph is not DAG, and no topological order can be found. Actually, this is not a DP problem. After listing out all equations of form ** F[i] = 1+\sum_{c=’A’}^{‘Z’}F[next(i,c)]/26**, we can already solve the linear system: using Gaussian-Elimination.

By now, the problem is mostly solved. However, when I submit the code to the judge, I got “wrong answer”. I guess this is because of precision problem: I use double data type, since the expected length can be real number. However, the required output is integer, thus I cast the result into “long long” data type. Real number data type can store 15-16 significant digits. However, long long data type can represent numbers with almost 20 significant digits. Thus, in extreme case, converting double may lose necessary precision when the answer is very large.

I fixed this problem by using “long double” instead of “double”. To my best knowledge, “long double” is some system dependent data type in C. It has higher precision than normal double: on my machine(gcc), sizeof(long double) gives 12, and actually 80 bits are used to store data, 16 more bits are for padding. I am still not sure how many significant digits “long double” can have, but this problem shows that at least it is useful in some cases.

#include <vector>
#include <iostream>
#include <algorithm>
#include <string>
using namespace std;
const long double EPI=1e-6;
string s;
int N;
long double abs(long double a){return a>0?a:-a;}
int g(int i,char c){
string t=s.substr(0,i)+c;
for(int j=0;j<i+1;j++){
if(t==s.substr(0,t.size()))
return i+1-j;
t=t.substr(1);
}
return 0;
}
void sub_from(vector<long double> & a,vector<long double> &b, long double k){
// a=a-bk
for(int i=0;i<a.size();i++)
a[i]-=b[i]*k;
}
vector<long double> solve_linear_system(vector<vector<long double> >& a){
for(int i=0;i<a.size();i++){
int j;
for(j=i;j<a.size();j++)
if(abs(a[j][i]) > EPI)
break;
swap(a[j],a[i]);
for(j=i+1;j<a.size();j++)
sub_from(a[j],a[i],a[j][i]/a[i][i]);
}
for(int i=a.size()-1;i>=0;i--){
for(int j=i-1;j>=0;j--)
sub_from(a[j],a[i],a[j][i]/a[i][i]);
}
vector<long double> ans;
for(int i=0;i<a.size();i++)
ans.push_back(a[i].back()/a[i][i]);
return ans;
}
long double solve(){
int m=s.size();
vector<vector<long double> > a;
for(int i=0;i<m;i++){
a.push_back(vector<long double>(m+2,0));
a[i][i]=N;
a[i][m+1]=N;
for(int j=0;j<N;j++)
a[i][g(i,'A'+j)]-=1;
}
a.push_back(vector<long double>(m+2,0));
a[m][m]=1; //xm=0
vector<long double> ans=solve_linear_system(a);
return ans[0]; //x0
}
int main(int argc, const char *argv[])
{
int times;
cin>>times;
for(int i=0;i<times;i++){
cin>>N>>s;
if(i)cout<<endl;
cout<<"Case "<<i+1<<":"<<endl<<(long long)solve()<<endl;
}
return 0;
}