func_code_string: #include <bits/stdc++.h>   using namespace std;   struct mat {     long double x[51][51];   } mi, ma;   int n, p;   long double pa, push_back, aa[51][51], bb[51][51], dp[51][51], ta, tb, ans;   mat tim(mat a, mat b) {     mat re;     for (int c = 0; c < 51; c++)       for (int d = 0; d < 51; d++) {         re.x[c][d] = 0;         for (int e = 0; e < 51; e++) re.x[c][d] += a.x[c][e] * b.x[e][d];       }     return re;   }   mat po(int p) {     if (p == 0) return mi;     if (p & 1) return tim(ma, po(p - 1));     mat tt = po(p / 2);     return tim(tt, tt);   }   int main() {     scanf( %d%d , &n, &p);     pa = p;     pa /= 1000000007 - 7;     push_back = 1 - pa;     for (int a = 0; a <= 50; a++) aa[a][0] = bb[a][0] = 1;     aa[1][1] = pa;     aa[1][2] = push_back;     bb[1][2] = 1;     for (int a = 2; a <= 50; a++) {       aa[a][1] = pa;       aa[a][2] = pa * pa + push_back;       bb[a][2] = 1;       for (int b = 3; b <= 50; b++) {         aa[a][b] = aa[a][b - 1] * aa[a - 1][b - 1];         bb[a][b] = bb[a][b - 1] * aa[a - 1][b - 1];       }     }     for (int a = 1; a <= 50; a++) dp[1][a] = a;     for (int a = 2; a <= 50; a++)       for (int b = 1; b <= 50; b++) {         ta = tb = 0;         if (b == 1) {           for (int c = 2; c <= 50; c++) {             ta += dp[a - 1][c] * bb[a - 1][c] * (1 - aa[a - 2][c]);             tb += bb[a - 1][c] * (1 - aa[a - 2][c]);           }         } else {           for (int c = 1; c < b; c++) {             ta += dp[a - 1][c] * aa[a - 1][c] * (1 - aa[a - 2][c]);             tb += aa[a - 1][c] * (1 - aa[a - 2][c]);           }         }         dp[a][b] = b + ta / tb;       }     if (n <= 50) {       for (int a = 1; a <= 50; a++)         ans += aa[n][a] * (1 - aa[n - 1][a]) * dp[n][a];     } else {       for (int a = 0; a < 51; a++) mi.x[a][a] = 1;       ma.x[0][0] = 1;       for (int a = 1; a <= 50; a++) {         ma.x[a][0] = a;         if (a == 1) {           ta = 0;           for (int b = 2; b <= 50; b++) {             ma.x[a][b] += bb[50][b] * (1 - aa[50][b]);             ta += bb[50][b] * (1 - aa[50][b]);           }           for (int b = 2; b <= 50; b++) ma.x[a][b] /= ta;         } else {           ta = 0;           for (int b = 1; b < a; b++) {             ma.x[a][b] += aa[50][b] * (1 - aa[50][b]);             ta += aa[50][b] * (1 - aa[50][b]);           }           for (int b = 1; b < a; b++) ma.x[a][b] /= ta;         }       }       ma = po(n - 50);       for (int a = 1; a <= 50; a++) {         ta = ma.x[a][0];         for (int b = 1; b <= 50; b++) ta += ma.x[a][b] * dp[50][b];         ans += aa[50][a] * (1 - aa[50][a]) * ta;       }     }     printf( %.20lf n , (double)ans);     return 0;   }   