Pollard_rho大数质因数分解+拉格朗日四平方和定理（bzoj 2904: 平方和）

2019-04-14 12:32发布

站内文章 / 模拟电子
11673 0 946

2904: 平方和

Time Limit: 10 Sec  Memory Limit: 128 MB
Submit: 160  Solved: 73
[Submit][Status][Discuss]

Description

给定一个整数N，求N最少可以拆成多少个完全平方数的和

Input

第一行一个整数TEST，表示数据组数
接下来TEST行，每行一个整数N

Output

TEST行，每行一个正整数，表示N最少可以拆成多少个完全平方数的和。  
0≤N≤10^18 ,TEST<=1000  Sample Input 5 413357 567226 766291 597007 289611

Sample Output

3 3 3 4 3

拉格朗日四平方和定理：每个正整数均可表示为4个整数的平方和 所以这道题的答案不会超过4 ①答案为1：直接对n开更号判断即可 ②答案为2：将n进行质因数分解（因为较大所以要用Pollard_rho算法），所有对4取模为3的质因数次数都是偶数 性质：假设n为两个答案为2的数相乘，那么Ans[n]还是等于2 ③答案为3：如果不满足其他所有条件，那么答案就为3了hh ④答案为4：n可以表示成4^x(8m+7)(x和m都为整数)的形式
#include #include #include using namespace std; #define LL long long int t, cnt, pri[50005], a[50005] = {1,1}; LL fat[101]; LL Multi(LL a, LL b, LL mod) { LL t; a %= mod, b %= mod; t = a*b-(LL)((double)a/mod*b+0.5)*mod; t = (t+mod)%mod; return t; } /* LL Multi(LL a, LL b, LL mod) //和上面一样，但比上面慢很多 { LL ans = 0; a %= mod; while(b) { if(b%2==1) ans = (ans+a)%mod, b--; else a = (a+a)%mod, b /= 2; } return ans; } */ LL Pow(LL a, LL b, LL mod) { LL ans = 1; a %= mod; while(b) { if(b&1) ans = Multi(ans, a, mod), b--; else a = Multi(a, a, mod), b /= 2; } return ans; } LL Gcd(LL a, LL b) { if(b==0) return a; return Gcd(b, a%b); } int Miller_Rabin(LL n) { int i, j, k; LL a, x, y, mod; if(n==2) return 1; if(n<2 || n%2==0) return 0; k = 0, mod = n-1; while(mod%2==0) { k++; mod /= 2; } for(i=1;i<=10;i++) { a = rand()%(n-1)+1; x = Pow(a, mod, n); y = 0; for(j=1;j<=k;j++) { y = Multi(x, x, n); if(y==1 && x!=1 && x!=n-1) return 0; x = y; } if(y!=1) return 0; } return 1; } LL Divi(LL n) { LL i, k, x, y, p, c; if(n==1) return 1; k = 2, p = 1; y = x = rand()%n, c = rand()%(n-1)+1; for(i=1;p==1;i++) { x = (Multi(x, x, n)+c)%n; p = abs(x-y); p = Gcd(n, p); if(i==k) y = x, k *= 2; } return p; } void Pollard_rho(LL n) { LL p; if(n==1) return; if(Miller_Rabin(n)) fat[++cnt] = n; else { p = Divi(n); Pollard_rho(n/p); Pollard_rho(p); } } int Jud2(LL n) { int i, sum; cnt = 0; Pollard_rho(n); sort(fat+1, fat+1+cnt); for(i=1;i<=cnt;i++) { if(fat[i]%4==3) { sum = 1; while(fat[i]==fat[i+1]) i++, sum++; if(sum%2==1) return 0; } } return 1; } int main(void) { LL T, n, i, j, t; for(i=2;i<50000;i++) { if(a[i]==0) pri[++cnt] = i; for(j=i*i;j<=50000;j+=i) a[j] = 1; } scanf("%lld", &T); while(T--) { scanf("%lld", &n); t = sqrt(n); if(t*t==n) printf("1 "); else if(Jud2(n)) printf("2 "); else { while((n&3)==0) n /= 4; if(n%8==7) printf("4 "); else printf("3 "); } } return 0; }

热门文章