Acdreamer博客数论学习 Day One 等比数列求和. 求解$S_n = (a+a^2+…+a^n)modM$
采用二分.n%2==0时 $ S_n = (1+a^{(\frac{n}{2})} )S_{\frac{n}{2}} $ ,n%2==1时,$S_n = (1+a^{\frac{n+1}{2}})S_{\frac{n-1}{2}} + a^{\frac{n+1}{2}}$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 #include <cstdio> #include <iostream> #define ll long long using namespace std ;const int mod = 7 ;ll pow (int a, int k) { ll ans = 1 ; while (k){ if (k&1 )ans=ans*a%mod; a=a*a%mod; k>>=1 ; } return ans; } ll S (int a, int n) { if (n==1 )return a%mod; if (n&1 ){ return ((1 +pow (a,(n+1 )/2 ))%mod*S(a,(n-1 )/2 )%mod + pow (a,(n+1 )/2 )%mod)%mod; }else { return (1 +pow (a,n/2 ))%mod*S(a,n/2 )%mod; } } int main () { int a,n; while (~scanf ("%d%d" ,&a,&n)){ printf ("%lld\n" ,S(a,n)); } return 0 ; }
题目:https://vjudge.net/problem/POJ-3233.
计算矩阵1-n次方和
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 int n,k,m;int mod = 1e9 +7 ;#include <cstdio> #include <iostream> #include <cstring> using namespace std ;struct Matrix { int mp[35 ][35 ]; Matrix(){ memset (mp,0 ,sizeof mp); } Matrix(int v){ memset (mp,0 ,sizeof mp); for (int i = 1 ; i <= n; i++){ mp[i][i] = v; } } }; Matrix multi (Matrix a, Matrix b) { Matrix c; for (int i = 1 ; i <= n; i++){ for (int k = 1 ; k <= n; k++){if (a.mp[i][k]) for (int j = 1 ; j <= n; j++){if (b.mp[k][j]) c.mp[i][j] = c.mp[i][j] + a.mp[i][k]*b.mp[k][j]; c.mp[i][j] %= mod; } } } return c; } Matrix pls (Matrix a, Matrix b) { Matrix c; for (int i = 1 ; i <= n; i++){ for (int j = 1 ; j <= n; j++){ c.mp[i][j] = a.mp[i][j]+b.mp[i][j]; c.mp[i][j] %= mod; } } return c; } Matrix powmul (Matrix a, int k) { Matrix c = Matrix(1 ); for (int i = 1 ; i <= n; i++) c.mp[i][i] = 1 ; while (k){ if (k&1 )c = multi(c,a); a = multi(a,a); k>>=1 ; } return c; } Matrix S (Matrix a, int k) { Matrix temp = Matrix(1 ); if (k==1 )return a; if (k&1 ){ Matrix temp2 = powmul(a,(k+1 )/2 ); return pls(multi(pls(temp, temp2), S(a,(k-1 )/2 )), temp2); }else { Matrix temp2 = powmul(a,k/2 ); return multi(pls(temp, temp2), S(a,k/2 )); } } int main () { while (~scanf ("%d%d%d" ,&n,&k,&m)){ Matrix a; for (int i = 1 ; i <= n; i++){ for (int j = 1 ; j <= n; j++){ scanf ("%d" ,&a.mp[i][j]); } } mod = m; a = S(a,k); for (int i = 1 ; i <= n; i++){ for (int j = 1 ; j <= n; j++){ printf ("%d" ,a.mp[i][j]); if (j!=n)printf (" " ); } puts ("" ); } } return 0 ; }
素数定理 $$ \pi(n)为小于等于n的素数的个数,有lim_{n→∞}\frac{\pi(n)}{n/lnn}=1 $$
定理$gcd(a^n-1,a^m-1)=a^{gcd(n,m)}-1$ 在a>1,m,n>0时成立.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 #include <cstdio> #include <iostream> #define ll long long using namespace std ;int mod = 1e9 +7 ;ll gcd (ll a, ll b) { return b?gcd(b,a%b):a; } ll pow (ll a, ll k) { ll ans = 1 ; while (k){ if (k&1 )ans=ans*a%mod; a = a*a%mod; k>>=1 ; } return ans; } int main () { int a,m,n,k1,t;scanf ("%d" ,&t); while (t--){ scanf ("%d%d%d%d" ,&a,&m,&n,&k1); mod = k1; ll ans = pow (a,gcd(n,m))-1 +mod; printf ("%lld\n" ,ans%mod); } return 0 ; }
定理扩展$a>b,gcd(a,b)=1, gcd(a^m-b^m,a^n-b^n)=a^{gcd(m,n)}-b^{gcd(m,n)}$ 定理$G=gcd(C_n^1,C_n^2,C_n^3,…,Cn^{n-1})$ n为素数,答案为n; n有多个素因子,答案为1; n只有一个素因子,答案就是该素因子(其实包括了第一种情况 题目:https://vjudge.net/problem/HDU-2582
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 #include <cstdio> #include <iostream> #include <cstring> #define ll long long #define maxn 1000100 using namespace std ;bool judge[maxn];int prime[maxn/10 ];ll biao[maxn]; int tot = 0 ;void init () { memset (prime,0 ,sizeof prime); memset (judge,0 ,sizeof judge); for (int i = 2 ; i < maxn; i++){ if (!judge[i]){ for (int j = i+i; j < maxn; j+=i){ judge[j]=1 ; } prime[tot++]=i; } } fill(biao,biao+maxn,1 ); for (int i = 0 ; i < tot; i++){ for (ll j = prime[i]; j < maxn; j*=prime[i]){ biao[j]=prime[i]; } } for (int i = 4 ; i < maxn; i++){ biao[i]+=biao[i-1 ]; } } int main () { int n; init(); while (~scanf ("%d" ,&n)){ printf ("%lld\n" ,biao[n]); } return 0 ; }
定理$F_n 为斐波那契额数列 ,gcd(F_m,F_n)=F_{gcd(m,n)}$ http://acm.nyist.net/JudgeOnline/problem.php?pid=468
此处需要学会一种新的判断素数的方式叫做Miller_Rabin素数测试.在DayTwo中有.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 #include <stdio.h> #include <algorithm> #include <iostream> using namespace std ;long long multi (long long a, long long b, long long mod) { long long temp = a, sum = 0 ; while (b){ if (b & 1 ) sum = (sum + temp) % mod; temp = (temp + temp) % mod; b = b >> 1 ; } return sum; } long long Modular_exponent (long long a, long long x, long long mod) { long long t = a % mod, r = 1 ; while (x){ if (x & 1 ) r = multi(r, t, mod); t = multi(t, t, mod); x = x >> 1 ; } return r; } bool Miller_Rabin (long long n) { long long time = 20 ; if (n < 2 ) return false ; if (n == 2 ) return true ; bool flag = false ; for (long long k = 0 ; k <= time; ++k){ flag = false ; long long d = n - 1 , r = 0 , t, a = rand() % (n - 2 ) + 2 ; while ((d & 1 ) == 0 ){ d = d>>1 ; r++; } t = Modular_exponent(a, d, n); if (t == 1 || t == n-1 ) {flag = true ; continue ;} for (long long i = 1 ; i < r; i++){ t = multi(t, t, n); if (t == 1 ) {flag = false ; return flag;} if (t == n-1 ) {flag = true ; break ;} } if (!flag) break ; } return flag; } int main () { long long n; while (scanf ("%lld" ,&n)!=EOF) puts ((n == 4 || (n != 2 && Miller_Rabin(n))) ? "Yes" : "No" ); return 0 ; }
定理:给定两个互素的正整数A和B,那么他们最大不能组合的数为AB-A-B,不能组合的个数为num=$\frac{(A-1)(B-1)}{2}$ ax+by=c的存在非负解的条件就是c>AB-A-B.
HDU - 1792
1 2 3 4 5 6 7 8 9 10 #include <cstdio> #include <iostream> using namespace std ;int main () { int a,b; while (~scanf ("%d%d" ,&a,&b)){ printf ("%d %d\n" ,a*b-a-b,(a-1 )*(b-1 )/2 ); } return 0 ; }
定理$\sum{gcd(i,N)=\sum{d\phi(\frac{N}{d})}}$ 数论专题里也有一题如是说.
https://vjudge.net/problem/POJ-2480
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 #include <cstdio> #include <cstring> #include <iostream> #define ll long long #define maxn 50000 using namespace std ;bool judge[maxn];int prime[maxn];int tot;void init () { tot = 0 ; memset (judge,0 ,sizeof judge); memset (prime, 0 , sizeof prime); for (int i = 2 ; i < maxn; i++){ if (!judge[i]){ for (int j = i+i; j < maxn; j+=i){ judge[j] = 1 ; } prime[tot++]=i; } } } int main () { init(); ll n; while (~scanf ("%lld" ,&n)){ ll ans = n; for (int i = 0 ; i < tot && prime[i]*prime[i] <= n; i++){ if (n%prime[i]==0 ){ int cot = 0 ; while (n%prime[i]==0 ){cot++;n/=prime[i];} ans += ans * cot * (prime[i]-1 )/prime[i]; } } if (n!=1 )ans += ans*(n-1 )/n; printf ("%lld\n" ,ans); } }
剩下没有题目的定理 $$ (n+1)lcm(C_n^0,C_n^1,C_n^2,…,C_N^{n-1},C_n^N)=lcm(1,2,…,n+1) $$
$$ 定理:任何n个连续的正整数的乘积均可被n!整除 $$
结论: $$ 如果p是素数,那么C_p^1,C_p^2,C_p^3,…C_p^{p-1},均能被p整除\ 如果p是素数,那么有(x+y+..+w)^p≡x^p+y^p+…+w^p(mod p)\ $$
Day Two 卢卡斯-莱默检验法 检验梅森素数 $$ 梅森素数M_p =2^p-1,p为素数.,定义序列{s_i},所有的i ≥ 0, \ si= \begin{cases} 4, i=0 \ s_{i-1}^{2}-2,i>0 \end{cases} M_p是素数当且仅当s_{p-2} \equiv 0(modM_p) $$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 #include <cstdio> #include <iostream> #define ll long long using namespace std ;ll T,p; ll multi (ll a, ll b, ll m) { ll ans = 0 ; while (b){ if (b&1 )ans = (ans+a)%m; a = (a+a)%m; b>>=1 ; } return ans; } int main () { scanf ("%lld" ,&T); while (T--){ scanf ("%lld" ,&p); ll n = ((ll) 1 <<p)-1 ; ll r = 4 ; if (p==2 )puts ("yes" ); else { while ((--p)!=1 )r=(multi(r,r,n)-2 +n)%n; if (r%n==0 )puts ("yes" ); else puts ("no" ); } } }
Miller-Rabin素数测试 费马小定理 :对于素数p和任意整数a,有$$a^p ≡ a(mod p)$$(同余)。反过来,满足$$a^p ≡ a(mod p)$$,p也几乎 一定是素数。
伪素数 :如果n是一个正整数,如果存在和n互素的正整数a满足 $a^{n-1} ≡ 1(mod n)$,我们说n是基于a的伪素数。如果一个数是伪素数,那么它几乎肯定是素数。
序列: a^m%n; a^(2m)%n; a^(4m)%n; …… ;a^(m*2^q)%n把上述测试序列叫做Miller测试,关于Miller测试,有下面的定理:
定理:若n是素数,a是小于n的正整数,则n对以a为基的Miller测试,结果为真.Miller测试进行k次,将合数当成素数处理的错误概率最多不会超过4^(-k).
Miller-Rabin测试 :不断选取不超过n-1的基b(s次),计算是否每次都有$$b^{n-1} ≡ 1(mod n)$$,若每次都成立则n是素数,否则为合数。
二次探测定理 :p为素数,0<x<p,则方程$x^2 ≡ 1(mod p)的解为1或者p-1$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 #include <cstdio> #include <cmath> #include <cstdlib> #include <iostream> #define ll long long using namespace std ;const int Times = 10 ;ll multi (ll a, ll b ,ll m) { ll ans = 0 ; a%=m; while (b){ if (b&1 )ans = (ans+a)%m; a = (a+a)%m; b >>= 1 ; } return ans; } ll quick_mod (ll a, ll b, ll m) { ll ans = 1 ; a %= m; while (b){ if (b&1 )ans = multi(ans, a, m); a = multi(a,a,m); b>>=1 ; } return ans; } bool Miller_Rabin (ll n) { if (n==2 )return true ; if ((n<2 ) || !(n&1 ))return false ; ll m = n-1 ; int k = 0 ; while ((m&1 )==0 ){ k++; m>>=1 ; } for (int i = 0 ; i < Times; i++){ ll a = rand()%(n-1 )+1 ; ll x = quick_mod(a,m,n); ll y = 0 ; for (int j = 0 ; j < k; j++){ y = multi(x,x,n); if (y==1 && x!=1 && x!= n-1 )return false ; x = y; } if (y!=1 )return false ; } return true ; } int main () { int T; scanf ("%d" ,&T); while (T--){ ll n; scanf ("%lld" ,&n); if (Miller_Rabin(n))puts ("Yes" ); else puts ("No" ); } return 0 ; }
例题:[http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1186
java写大数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 import java.math.BigInteger;import java.util.Random;import java.util.Scanner;public class test { public static final int Times = 10 ; public static BigInteger quick_mod (BigInteger a, BigInteger b, BigInteger m) { BigInteger ans = BigInteger.ONE; a = a.mod(m); while (!b.equals(BigInteger.ZERO)){ if (b.mod(BigInteger.valueOf(2 )).equals(BigInteger.ONE)){ ans = (ans.multiply(a)).mod(m); b = b.subtract(BigInteger.ONE); } b = b.divide(BigInteger.valueOf(2 )); a = (a.multiply(a)).mod(m); } return ans; } public static boolean Miller_Rabin (BigInteger n) { if (n.equals(BigInteger.valueOf(2 )))return true ; if (n.equals(BigInteger.ONE))return false ; if ((n.mod(BigInteger.valueOf(2 ))).equals(BigInteger.ZERO))return false ; BigInteger m = n.subtract(BigInteger.ONE); BigInteger y = BigInteger.ZERO; int k = 0 ; while ((m.mod(BigInteger.valueOf(2 ))).equals(BigInteger.ZERO)){ k++; m = m.divide(BigInteger.valueOf(2 )); } Random d = new Random(); for (int i = 0 ; i < k; i++){ int t = 0 ; if (n.compareTo(BigInteger.valueOf(10000 ))==1 ){ t = 10000 ; }else { t = n.intValue()-1 ; } int a = d.nextInt(t) + 1 ; BigInteger x = quick_mod(BigInteger.valueOf(a), m, n); for (int j = 0 ; j < k; j++){ y = (x.multiply(x)).mod(n); if (y.equals(BigInteger.ONE) && !(x.equals(BigInteger.ONE)) && !(x.equals(n.subtract(BigInteger.ONE)))){ return false ; } x = y; } if (!(y.equals(BigInteger.ONE)))return false ; } return true ; } public static void main (String[] args) { Scanner in = new Scanner(System.in); while (in.hasNextBigInteger()){ BigInteger n = in.nextBigInteger(); if (Miller_Rabin(n))System.out.println("Yes" ); else System.out.println("No" ); } } }
扩展欧几里得 二元一次方程$$ax_1+by_1=gcd(a,b)$$;
a>b时,b==0时,x=1,y=0, 否则设$bx_2+(a\ mod\ b)y_2 = gcd(b, a\ mod \ b)=gcd(a,b)$
又有$ax_1 +by_1=bx_2+(a-[\frac{a}{b}]b)y_2$
则有$x_1=y_2,y_1=x_2-[\frac{a}{b}]y_2$, 用递推实现.
1 2 3 4 5 6 7 8 9 void e_gcd (int a, int b, int &d, int &x, int &y) { if (!b){ d=a,x=1 ,y=0 ; return ; }else { e_gcd(b,a%b,d,y,x); y -= (a/b)*x; } }
代表的是ax+by=gcd(a,b)的解x,y,和d为gcd(a,b);,得到的x,y中a>b则x为正.在gcd(a,b)==1的时候,可以领y=(y%a+a)%a把y变成正的,再根据等式把x变成负的就行.
例题:https://vjudge.net/problem/POJ-2142
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 #include <cstdio> #include <iostream> #include <cmath> using namespace std ;int gcd (int a, int b) { return b?gcd(b,a%b):a; } void e_gcd (int a, int b, int &d, int &x, int &y) { if (!b){ d=a,x=1 ,y=0 ; return ; }else { e_gcd(b,a%b,d,y,x); y -= (a/b)*x; } } int main () { int a,b,c; while (~scanf ("%d%d%d" ,&a,&b,&c) && (a+b+c)){ int x, y, d, tx, ty; int gd = gcd(a,b); a /= gd; b/= gd; c/= gd; e_gcd(a,b,d,x,y); tx = x*c; tx = (tx%b + b)%b; ty = (c-tx*a)/b; if (ty < 0 )ty = -ty; y = y*c; y = (y%a + a)%a; x = (c-b*y)/a; if (x<0 )x =-x; if (x+y>tx+ty){ x=tx; y=ty; } printf ("%d %d\n" ,x,y); } return 0 ; }
https://vjudge.net/problem/HIT-2815
题目没读清,特判就写错了,扩展欧几里得和上面一道题一样的用.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 #include <cstdio> #include <iostream> #define ll long long using namespace std ;ll gcd (ll a, ll b) { return b?gcd(b, a%b):a; } void e_gcd (ll a, ll b, ll &d, ll &x, ll &y) { if (!b){ d = a; x = 1 ; y = 0 ; return ; }else { e_gcd(b,a%b,d,y,x); y -= (a/b)*x; } } int main () { ll t,a,b,d,x,y;scanf ("%lld" ,&t); while (t--){ scanf ("%lld%lld" ,&a,&b); if (a*b==0 ){ if (a==1 || b==1 ){ puts ("1" ); }else { puts ("-1" ); } }else if (a==1 || b==1 ){ if (max(a,b)==2 ) puts ("1" ); else puts ("2" ); } else if (gcd(a,b)!=1 ){ puts ("-1" ); }else { e_gcd(a,b,d,x,y); ll tx,ty; tx = (x%b+b)%b; ty = (1 -a*x)/b; if (ty<0 )ty = -ty; y = (y%a+a)%a; x = (1 -b*y)/a; if (x<0 )x=-x; if (x+y>tx+ty){ x = tx; y = ty; } printf ("%lld\n" ,x+y-1 ); } } return 0 ; }
毕达哥拉斯三元组的解 $x^2+y^2=z^2$,只考虑本原解.不妨假设x为偶数,y,z为奇数.这样就必有$x=2mn,y=m^2-n^2,z=m^2+n^2$.
例题:https://vjudge.net/problem/POJ-1305
输出两个数据,一个是z不大于n的本原解的个数,和本原解中没有用到的数的个数.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 #include <cstdio> #include <cstring> #include <iostream> #define maxn 1000111 using namespace std ;int biao[maxn];bool flag[maxn];int gcd (int a, int b) { return b?gcd(b,a%b):a; } int main () { int n; memset (flag,0 ,sizeof flag); for (int i = 2 ; i*i < maxn; i++){ for (int j = 1 ; j < i; j++){ if (gcd(i,j)!=1 )continue ; if ((i&1 )&&(j&1 ))continue ; if (i*i+j*j<maxn)biao[i*i+j*j]++; } } for (int i = 1 ; i < maxn; i++){ biao[i]+=biao[i-1 ]; } while (~scanf ("%d" ,&n)){ memset (flag,0 ,sizeof flag); int ans = 0 ; for (int i = 2 ; i*i < maxn; i++){ for (int j = 1 ; j < i; j++){ if (gcd(i,j)!=1 )continue ; if ((i&1 )&&(j&1 ))continue ; for (int k = 1 ; k*(i*i+j*j)<=n; k++){ flag[2 *i*j*k]=flag[k*(i*i-j*j)]=flag[k*(i*i+j*j)]=1 ; } } } for (int i = 1 ; i <= n; i++) if (flag[i])ans++; printf ("%d %d\n" ,biao[n],n-ans); } return 0 ; }
例题二:
利用毕达哥拉斯三元组+欧拉函数+容斥原理
https://vjudge.net/problem/HDU-3939?? ?
寻找满足的i,j满足gcd(i,j)=1,dfs()容斥原理来找有多少互斥的数.
这里需要判断i<=j和i>j的情况.其中,i<=j时,i为偶数,则直接加phi[i],i为奇数时,由于j必须为偶数,不能用phi[i],所以用容斥原理,一共i>>1个数(晒去了2的倍数).大于时同理.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 #include <cstdio> #include <iostream> #include <cstring> #include <cmath> #define ll long long #define maxn 1001000 using namespace std ;bool judge[maxn];int prime[maxn/10 ];int phi[maxn];int check[40 ];int tot,num;ll ans; void init () { tot = 0 ; memset (judge, 0 , sizeof judge); for (int i = 2 ; i < maxn; i++){ if (!judge[i]){ for (int j = i+i; j < maxn; j+=i){ judge[j]=1 ; } prime[tot++]=i; } } for (int i = 1 ; i < maxn; i++)phi[i]=i; for (int i = 2 ; i < maxn; i+=2 )phi[i]>>=1 ; for (int i = 3 ; i < maxn; i+=2 )if (phi[i]==i){ for (int j = i; j < maxn; j+=i){ phi[j] = phi[j] - phi[j]/i; } } } void prime_check (int n) { num = 0 ; if (!judge[n]){ check[num++]=n; return ; } for (int i = 0 ; i < tot && n > 1 ; i++){ if (n%prime[i]==0 ){ check[num++]=prime[i]; while (n%prime[i]==0 ) n/=prime[i]; if (n>1 && !judge[n]){ check[num++]=n; return ; } } } } void dfs (int k, int r, int s, int n) { if (k==num){ if (r&1 ) ans -= n/s; else ans+= n/s; return ; } dfs(k+1 , r, s, n); dfs(k+1 , r+1 , s*check[k], n); } ll L; int main () { init(); int t;scanf ("%d" ,&t); while (t--){ ans = 0 ; scanf ("%lld" ,&L); int m = (int )sqrt (1.0 *L); for (int i = m; i > 0 ; i--){ int p = (int )sqrt (L - (ll)i*i); if (i<=p){ if (i&1 ){ prime_check(i); dfs(0 ,0 ,1 ,i>>1 ); }else ans+=phi[i]; }else { prime_check(i); if (i&1 )dfs(0 ,0 ,1 ,p>>1 ); else dfs(0 ,0 ,1 ,p); } } printf ("%lld\n" ,ans); } }
欧拉函数与欧拉定理 定理一:$$m与n互素,\phi(mn)=\phi(m)\phi(n)$$, 所以n为奇数时,我们有$\phi(2n)=\phi(n)$
定理二: p是素数,有$\phi(p^a)=p^a-p^{a-1}$,容斥原理.
定理三: $n=p_1^{\alpha_1}p_2^{\alpha_2}…p_k^{\alpha_k}时,\phi(n)=n(1-\frac{1}{p_1})(1-\frac{1}{p_2})…(1-\frac{1}{p_k})$ 定理1+3
定理四: $\sum_{d|n}{\phi(d)}= n$ 莫比乌斯反演.
定理五: (a,m)=1,则$x \equiv ba^{\phi(m)-1} 是同于方程 ax \equiv b(mod \ m)$
定理六: n>2时,$\phi(n)$为偶数
定理三求解欧拉函数 :
1 2 3 4 5 6 7 8 9 10 11 int phi (int n) { int ans = n; for (int i = 2 ; i*i <= n; i++){ if (n%i==0 ){ ans = ans - ans/i; while (n%i==0 ) n/=i; } } if (n>1 )ans = ans-ans/n; return ans; }
容斥原理,筛法去减不满足的.
1 2 3 4 5 6 7 8 9 10 11 12 13 #define maxn 100000 int p[maxn];void phi_table () { memset (p,0 ,sizeof p); for (int i = 1 ; i < maxn; i++)p[i]=i; for (int i = 2 ; i < maxn; i+=2 )p[i]>>=1 ; for (int i = 3 ; i < maxn; i+=2 ){ if (p[i]==i){ for (int j = i; j < maxn; j+=i) p[j]=p[j]-p[j]/i; } } }
二次筛法 这道题有类似的在kuangbin的数论基础专题里有出现过.
例题:https://vjudge.net/problem/POJ-2689
注意1 2这个样例...还有晒的时候不能把质数也晒了....
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 #include <cstdio> #include <iostream> #include <cstring> #define maxn 50000 int prime[maxn];bool judge2[21 *maxn];bool judge1[maxn];using namespace std ;int tot=0 ;void init () { memset (judge1,0 ,sizeof judge1); for (int i = 2 ; i<maxn; i++){ if (!judge1[i]){ prime[tot++] = i; for (int j = i+i; j < maxn; j+=i){ judge1[j]=true ; } } } } int main () { init(); unsigned int u,v; while (~scanf ("%d%d" ,&u,&v)){ if (u==1 )u=2 ; memset (judge2,0 ,sizeof judge2); for (int i = 0 ; i < tot; i++){ if (prime[i]>v)break ; unsigned int j=u/prime[i]*prime[i]; if (u%prime[i]==0 && u!=prime[i])judge2[0 ]=1 ; while (j<=u)j+=prime[i]; if (j==prime[i])j+=j; for (;j<=v;j+=prime[i]){ judge2[j-u]=1 ; } } int st,stlen=maxn,ed,edlen=0 ; int cot=0 ,len=0 ; for (int i = 0 ; i <= v-u; i++){ if (!judge2[i]){ if (len &&len<stlen){ st = i-len; stlen = len; } if (len && len>edlen){ ed = i-len; edlen = len; } cot++; len = 0 ; } len ++ ; } if (cot!=2 && st==ed){ ed = st+stlen; } if (cot>1 )printf ("%d,%d are closest, %d,%d are most distant.\n" ,u+st,u+st+stlen,u+ed,u+ed+edlen); else puts ("There are no adjacent primes." ); } return 0 ; }
逆元 正整数a,m,如果有$ax \equiv 1(mod \ m)$,这个同于方程中x的最小正整数解叫做a模m的逆元.
逆元一般用扩展欧几里得算法来求.m为素数时,由费马小定理得到逆元为$a^{m-2}modm$
逆元常见问题
已知b|a,求解$ans = \frac{a}{b}mod \ m$
常见解法:扩展欧几里得.m是素数时,直接用费马小定理.但他们都需要a与m互素.
除法取模求逆元的通用方法: $ans =\frac{ a \ mod(mb)}{b}$ $$ \frac{a}{b}=km+x(x<m),a=kbm+bx,a\ mod(bm)=bx,x = \frac{a\ mod(bm)}{b} $$ 例题:https://vjudge.net/problem/POJ-1845
求$A^B$的所有因子和对9901取余的结果.
考虑$$A=p_1^{\alpha_1}p_2^{\alpha_2}…p_k^{\alpha_k} ,A^B=p_1^{\alpha_1 B}p_2^{\alpha_2B}…p_k^{\alpha_kB}$$
所有的因子和:$$sum=\prod_{1\le i \le k} \sum_{0 \le j \le \alpha_iB} p_i^{j}$$
可以用二分求等比数列和,也可以用等比数列公式(需要用逆元)直接求.分别写在了Method_one和Method_two中.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 #include <cstdio> #include <iostream> #include <cstring> #include <cmath> #define ll long long #define maxn 10000 using namespace std ;ll mod =9901 ; bool judge[maxn];ll check[maxn][2 ]; int prime[maxn/2 ];int tot,num;ll a,b; void init () { tot=0 ; memset (judge, 0 , sizeof judge); for (int i = 2 ; i < maxn; i++){ if (!judge[i]){ for (int j = i+i; j < maxn; j+=i){ judge[j]=1 ; } prime[tot++]=i; } } } ll multi (ll a, ll b, ll m) { ll ans = 0 ; a %= m; while (b){ if (b&1 )ans = (ans+a)%m; a = (a+a)%m; b >>= 1 ; } return ans; } ll largequickmod (ll p, ll k, ll m) { ll ans = 1 ; while (k){ if (k&1 )ans = multi(ans,p,m); p = multi(p,p,m); k>>=1 ; } return ans; } ll quickmod (ll p, ll k) { ll ans=1 ; p%=mod; while (k){ if (k&1 )ans = ans*p%mod; p = (p*p)%mod; k>>=1 ; } return ans; } ll S (ll p, ll k) { if (p==0 )return 0 ; if (k==0 )return 1 ; ll temp = quickmod(p,(k>>1 ))%mod; if (k&1 ){ return S(p,k>>1 )%mod*(1 +temp*p%mod)%mod; }else { return (S(p,(k>>1 )-1 )%mod*(1 +temp*p%mod)%mod+temp)%mod; } } void Method_one () { ll ans = 1 ; for (int i = 0 ; i < num; i++){ ans = (ans * S(check[i][0 ], b*check[i][1 ]))%mod; } if (a==0 )ans=0 ; printf ("%lld\n" ,ans); } void Method_two () { ll ans = 1 ; for (int i = 0 ; i < num; i++){ ll M = mod*(check[i][0 ]-1 ); ans *= (largequickmod(check[i][0 ], check[i][1 ]*b+1 , M)+M-1 )/(check[i][0 ]-1 ); ans %= mod; } printf ("%lld\n" ,ans); } int main () { init(); while (~scanf ("%lld%lld" ,&a,&b)){ ll temp = a; num=0 ; memset (check, 0 , sizeof check); for (int i = 0 ; i < tot; i++){ if (temp<=1 )break ; if (temp%prime[i]==0 ){check[num++][0 ]=prime[i];} while (temp%prime[i]==0 ){ temp/=prime[i]; check[num-1 ][1 ]++; } } if (temp>1 ){ check[num][0 ]=temp; check[num++][1 ]=1 ; } Method_two(); } return 0 ; }
1→p模p的所有逆元,这里p为奇质数.用快速幂的复杂度O(plog(p)),用递推式可以优化到O(p): $$ inv[i] = (M-M/i)inv[M\%i] \%M \ 推导:t=M/i,k=M\%i,有t i+k \equiv 0 (modM), -ti \equiv k(modM)\ 两边同时除i k记作:-tinv[k] \equiv invi \ 即 inv[i] = (M- M/i) inv[M\%i]\%M $$ 初始化inv[1]=1,就可以通过递推法求出1→p模奇素数的所有逆元.
http://www.lydsy.com/JudgeOnline/problem.php?id=2186
例题:1-N!中与M!互质的数的个数,M<=N;M!|N!, .当n是m的倍数时,1-n中与m互素的个数为$\frac{n}{m}\phi(m)$
题目中 $\frac{N!}{M!}\phi(M!) = \frac{N!}{M!}M!(1-\frac{1}{p_1})(1-\frac{1}{p_2})…(1-\frac{1}{p_k})$.
时间卡的特别紧.bool开的就超时了.inv[i]算逆元.multi[maxn]算阶乘取模,要求p与m互质,否则需要用除法取模.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 #include <cstdio> #include <bitset> #define ll long long #define maxn 10000005 using namespace std ;ll mod; bitset <maxn> judge;void init () { judge.set (); for (int i = 4 ; i < maxn; i+=2 )judge[i]=0 ; for (int i = 3 ; i < maxn; i+=2 ){ if (judge[i]){ for (int j = i+i; j < maxn; j+=i){ judge[j]=0 ; } } } } ll multi[maxn]; int inv[maxn];ll ans[maxn]; int main () { init(); int T,R; scanf ("%d%d" ,&T,&R); mod = R; multi[0 ]=1 ; for (int i = 1 ; i < maxn; i++){ multi[i]=multi[i-1 ]*i%mod; } inv[1 ]=1 ; for (int i = 2 ; i < mod && i < maxn; i++){ inv[i]=(mod-mod/i)*inv[mod%i]%mod; } ans[1 ]=1 ; for (int i = 2 ; i < maxn; i++){ if (judge[i]){ ans[i] = ans[i-1 ]*(i-1 )%mod; ans[i] = ans[i]*inv[i%mod]%mod; }else { ans[i] = ans[i-1 ]; } } int n,m; while (T--){ scanf ("%d%d" ,&n,&m); ll anss = 1L L*multi[n]*ans[m]%mod; printf ("%lld\n" ,anss); } return 0 ; }
欧拉筛 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 void getphi () { int i,j; phi[1 ]=1 ; for (i=2 ;i<=N;i++) { if (!mark[i]) { prime[++tot]=i; phi[i]=i-1 ; } for (j=1 ;j<=tot;j++) { if (i*prime[j]>N) break ; mark[i*prime[j]]=1 ; if (i%prime[j]==0 ) { phi[i*prime[j]]=phi[i]*prime[j];break ; } else phi[i*prime[j]]=phi[i]*(prime[j]-1 ); } } }