1 /* 用于求整数解得方程组. */
  2
  3 #include <iostream>
  4 #include<cstdio>
  5 #include <cstring>
  6 #include <cmath>
  7
  8 #include<map>
  9
 10 #define mp make_pair
 11 using namespace std;
 12
 13 const int maxn = 333;
 14
 15 int equ, var; // 有equ个方程，var个变元。增广阵行数为equ, 分别为0到equ - 1，列数为var + 1，分别为0到var.
 16 double a[maxn][maxn];
 17 double x[maxn]; // 解集.
 18 bool free_x[maxn]; // 判断是否是不确定的变元.
 19 int free_num;
 20
 21 map<pair<int,int>,int> S;
 22 int lit[maxn];
 23 int big[maxn];
 24
 25 void Debug(void)
 26 {
 27     int i, j;
 28     for (i = 0; i < equ; i++)
 29     {
 30         for (j = 0; j < var + 1; j++)
 31         {
 32             printf("%6.2lf",a[i][j]);
 33         }
 34         cout << endl;
 35     }
 36     cout << endl;
 37 }
 38
 39 //inline int gcd(int a, int b)
 40 //{
 41 //    int t;
 42 //    while (b != 0)
 43 //    {
 44 //        t = b;
 45 //        b = a % b;
 46 //        a = t;
 47 //    }
 48 //    return a;
 49 //}
 50 //
 51 //inline int lcm(int a, int b)
 52 //{
 53 //    return a * b / gcd(a, b);
 54 //}
 55
 56 // 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)
 57 int Gauss(void)
 58 {
 59     int i, j, k;
 60     int max_r; // 当前这列绝对值最大的行.
 61     int col; // 当前处理的列.
 62     double ta, tb;
 63     int LCM;
 64     double teS;
 65     int free_x_num;
 66     int free_index;
 67     // 转换为阶梯阵.
 68     col = 0; // 当前处理的列.
 69     for (k = 0; k < equ && col < var; k++, col++)
 70     {
 71         // 枚举当前处理的行.
 72         // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
 73         max_r = k;
 74         for (i = k + 1; i < equ; i++)
 75         {
 76             if (fabs(a[i][col]) > fabs(a[max_r][col])) max_r = i;
 77         }
 78         if (max_r != k)
 79         {
 80             // 与第k行交换.
 81             for (j = k; j < var + 1; j++) swap(a[k][j], a[max_r][j]);
 82         }
 83         if (a[k][col] == 0)
 84         {
 85             // 说明该col列第k行以下全是0了，则处理当前行的下一列.
 86             k--;
 87             continue;
 88         }
 89         for (i = k + 1; i < equ; i++)
 90         {
 91             // 枚举要删去的行.
 92             if (a[i][col] != 0)
 93             {
 94 //                LCM = lcm(fabs(a[i][col]), fabs(a[k][col]));
 95 //                ta = LCM / fabs(a[i][col]), tb = LCM / fabs(a[k][col]);
 96                 ta=1.0;
 97                 tb=fabs(a[i][col])/fabs(a[k][col]);
 98                 if (a[i][col] * a[k][col] < 0) tb = -tb; // 异号的情况是两个数相加.
 99                 for (j = col; j < var + 1; j++)
100                 {
101                     a[i][j] = a[i][j] * ta - a[k][j] * tb;
102                 }
103             }
104         }
105     }
106     //Debug();
107     // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
108     for (i = k; i < equ; i++)
109     {
110         // 对于无穷解来说，如果要判断哪些是自由变元，那么初等行变换中的交换就会影响，则要记录交换.
111         if (a[i][col] != 0) return -1;
112     }
113     // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
114     // 且出现的行数即为自由变元的个数.
115     if (k < var)
116     {
117         // 首先，自由变元有var - k个，即不确定的变元至少有var - k个.
118         for (i = k - 1; i >= 0; i--)
119         {
120             // 第i行一定不会是(0, 0, ..., 0)的情况，因为这样的行是在第k行到第equ行.
121             // 同样，第i行一定不会是(0, 0, ..., a), a != 0的情况，这样的无解的.
122             free_x_num = 0; // 用于判断该行中的不确定的变元的个数，如果超过1个，则无法求解，它们仍然为不确定的变元.
123             for (j = 0; j < var; j++)
124             {
125                 if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
126             }
127             if (free_x_num > 1) continue; // 无法求解出确定的变元.
128             // 说明就只有一个不确定的变元free_index，那么可以求解出该变元，且该变元是确定的.
129             teS = a[i][var];
130             for (j = 0; j < var; j++)
131             {
132                 if (a[i][j] != 0 && j != free_index) teS -= a[i][j] * x[j];
133             }
134             x[free_index] = teS / a[i][free_index]; // 求出该变元.
135             free_x[free_index] = 0; // 该变元是确定的.
136         }
137         return var - k; // 自由变元有var - k个.
138     }
139     // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
140     // 计算出Xn-1, Xn-2 ... X0.
141     for (i = var - 1; i >= 0; i--)
142     {
143         teS = a[i][var];
144         for (j = i + 1; j < var; j++)
145         {
146             if (a[i][j] != 0) teS -= a[i][j] * x[j];
147         }
148         //if (teS % a[i][i] != 0) return -2; // 说明有浮点数解，但无整数解.
149         x[i] = teS / a[i][i];
150     }
151     return 0;
152 }
153
154 int main(void)
155 {
156     int i, j,k;
157     int lvln=21;//21
158     int fin;
159     double P,Q;
160     k=0;
161     for(i=0; i<lvln-1; i++)
162         for(j=0; j<=i; j++)
163         {
164             lit[k]=j;
165             big[k]=i;
166             S[mp(j,i)]=k;
167             k++;
168         }
169     lit[k]=lvln-2;
170     big[k]=lvln-1;
171     S[mp(lit[k],big[k])]=k;
172     k++;
173     equ=k;
174     var=k;
175     while (scanf("%lf", &P) != EOF)
176     {
177         Q=1.0-P;
178         memset(a, 0, sizeof(a));
179         memset(x, 0, sizeof(x));
180         memset(free_x, 1, sizeof(free_x)); // 一开始全是不确定的变元.
181 //        for (i = 0; i < equ; i++)
182 //            for (j = 0; j < var + 1; j++)
183 //                scanf("%d", &a[i][j]);
184         a[k-1][k-1]=1;
185         a[k-1][var]=0;
186         for(i=0; i<k-1; i++)
187         {
188             //cout<<i<<'.'<<lit[i]<<','<<big[i]<<endl;
189             if(lit[i]!=big[i]) a[i][S[mp(lit[i]+1,big[i])]]+=P;
190             else a[i][S[mp(big[i],lit[i]+1)]]+=P;
191             a[i][S[mp(max(lit[i]-2,0),big[i])]]+=Q;
192             //cout<<max(lit[i]-2,0)<<','<<big[i]<<','<<S[mp(max(lit[i]-2,0),big[i])]<<endl;
193             a[i][var]-=1.0;
194             a[i][i]-=1.0;
195         }
196         //Debug();
197         free_num = Gauss();
198         //Debug();
199 //        if (free_num == -1) printf("无解!\n");
200 //        else if (free_num == -2) printf("有浮点数解，无整数解!\n");
201 //        else if (free_num > 0)
202 //        {
203 //            printf("无穷多解! 自由变元个数为%d\n", free_num);
204 //            for (i = 0; i < var; i++)
205 //            {
206 //                if (free_x[i]) printf("x%d 是不确定的\n", i + 1);
207 //                else printf("x%d: %d\n", i + 1, x[i]);
208 //            }
209 //        }
210 //        else
211 //        {
212 //            for (i = 0; i < var; i++)
213 //            {
214 //                printf("x%d %d,%d: %lf\n", i + 1,lit[i],big[i], x[i]);
215 //            }
216 //        }
217 //        printf("\n");
218         printf("%.6lf\n",x[0]);
219     }
220     return 0;
221 }