Split Divisibilities (Project Euler 598)

题目大意：

求将$100!$ 拆成$a*b$的方案数,其中$a<=b$并且它们的约数个数一样多。

 

思路：

先将$100!$质因数分解, 结果如图：

首先想到一个暴力DP, dp[i][j][k]表示考虑完前i个质数, 目前a有j个约数,b有k个约数的方案数。 用map保存状态。

答案就是sum(dp[25][j][j]).    

但是状态数会很多(大概有1e8个状态),所以考虑 中途相遇法。 对前3个质数做一次DP, 然后对后面22个质数做一次DP。 

最后答案就是 sum (dp1[3][i1][j1] * dp2[22][i2][j2])   条件是 i1 * i2 = j1 * j2.  即  i1 / j1 =  j2 / i2 .

一个优化是只保存 j和k互质的状态。 然后 最后 答案的时候 枚举 i1,j1,  在 dp2中 查找 j2 / i2 = i1 / j1的点 。

因为a<=b,所以最后答案还需要除以2.

代码：

 

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <algorithm>
  4 #include <cmath>
  5 #include <set>
  6 #include <cstring>
  7 #include <map>
  8 #include <queue>
  9 using namespace std;
 10 
 11 typedef long long ll;
 12 #define N 10000000
 13 #define M 1100
 14 typedef pair<int,int> pii;
 15 
 16 
 17 bool flag[N];
 18 int p[N],phi[N];
 19 
 20 struct node
 21 {
 22     ll x,y;
 23     bool operator < (const node &t)const
 24     {
 25         return y*t.x<x*t.y;
 26     }
 27     node (ll _x = 0, ll _y = 0){x = _x; y = _y;}
 28 };
 29 
 30 map<node, ll> mp;
 31 
 32 ll Gcd(ll x, ll y)
 33 {
 34     ll tmp;
 35     while (y)
 36     {
 37         tmp = x % y;
 38         x = y, y = tmp;
 39     }
 40     return x;
 41 }
 42 
 43 void Get_Primes(int lim)
 44 {
 45     phi[1]=1;
 46     for (int i=2;i<=lim;i++)
 47     {
 48         if (!flag[i]) p[++p[0]]=i,phi[i]=i-1;
 49         for (int j=1;j<=p[0] && i*p[j]<=lim;j++)
 50         {
 51             flag[i*p[j]]=true;
 52             if (i%p[j]==0) 
 53             {
 54                 phi[i*p[j]]=phi[i]*p[j];
 55                 break;
 56             }
 57             else phi[i*p[j]]=phi[i]*(p[j]-1);
 58         }
 59     }
 60 }
 61 
 62 map<pair<ll,ll>, ll> f[6], g[23];
 63 int cnt[30];
 64 
 65 int main()
 66 {
 67     freopen("in.in","r",stdin);
 68     freopen("out.out","w",stdout);
 69 
 70     int n = 100;
 71     Get_Primes(n); 
 72     for (int i = 1; i <= p[0]; ++i)
 73     {
 74         int  x = p[i];
 75         while (x <= n) cnt[i] += n / x, x *= p[i];
 76     }
 77     
 78     f[0][make_pair(1,1)] = 1;
 79     for (int i = 1; i <= 3; ++i)
 80     {
 81         for (map<pair<ll,ll>, ll>::iterator it = f[i - 1].begin(); it != f[i - 1].end(); it++)
 82         {
 83             ll k1 = (*it).first.first, k2 = (*it).first.second;
 84             for (int j = 0; j <= cnt[i]; ++j)
 85             {
 86                 ll kx = k1 * (j + 1), ky = k2 * (cnt[i] - j + 1) , d = Gcd(kx, ky);
 87                 f[i][make_pair(kx / d, ky / d)] += (*it).second;
 88             }
 89         }
 90     }
 91     g[0][make_pair(1,1)] = 1;
 92     for (int i = 1; i <= p[0] - 3; ++i)
 93     {
 94         for (map<pair<ll,ll>, ll>::iterator it = g[i - 1].begin(); it != g[i - 1].end(); it++)
 95         {
 96             ll k1 = (*it).first.first, k2 = (*it).first.second;
 97             for (int j = 0; j <= cnt[i + 3]; ++j)
 98             {
 99                 ll kx = k1 * (j + 1), ky = k2 * (cnt[i + 3] - j + 1), d = Gcd(kx, ky);
100                 g[i][make_pair(kx / d, ky / d)] += (*it).second;
101             }
102         }
103     }
104     for (map<pair<ll,ll>, ll>::iterator it = g[p[0] - 3].begin(); it != g[p[0] - 3].end(); it++)
105     {
106         ll x = (*it).first.first, y = (*it).first.second;
107         mp[node(x, y)] += (*it).second;
108     }
109     
110     ll res = 0;
111     for (map<pair<ll,ll>, ll>::iterator it = f[3].begin(); it != f[3].end(); it++)
112     {
113         ll x = (*it).first.first, y = (*it).first.second;
114         res += (*it).second * mp[node(y, x)];
115     }
116     cout << res / 2 << endl;
117     return 0;
118 }
View Code

答案：543194779059

posted @ 2017-04-26 13:32  lzw4896s  阅读(435)  评论(0)  编辑  收藏  举报