数论算法

数论

素数问题、同余问题、中国剩余定理、Nim积、高斯消元法求线性方

程组解、Pell方程、polya计数、矩阵二分快速幂、伪素数、基本数

值计算方法（定积分求解，多项式求根，周期性方程）

1、与整数除法运算相关的算法

整除问题：

欧几里得算法及利用其求最小公倍数：

拓展欧几里得算法及利用其解线性同余方程：

a^b%c几种计算方法

中国剩余定理

#include

using namespace std;

int ext_gcd(int a, int b, int &x, int &y)

{

int tmp,d;

if(b == 0)

{

x = 1;

y = 0;

return a;

}

d = ext_gcd(b, a % b, x, y);

tmp = x;

x = y;

y = tmp - a/b*y;

return d;

}

int China_theory(int a[],int b[],int n)

{

int res = 0,m,*m1,M = 1,temp;

m1 = new int[n];

memset(m1, 0, sizeof(m1));

for(int i = 0; i < n; i++)

M *= a[i];

for(int i = 0; i < n; i++)

{

m = M / a[i];

ext_gcd(m, a[i], m1[i],temp);

res = res % M + (m * m1[i] * b[i]) % M;--

res =(res + M) % M;

}

delete m1;

return res;

}

int main()

{

int *a, *b, n;

while(scanf("%d",&n) && n != 0)

{

a = new int[n];

b = new int[n];

for(int i = 0; i < n; i++)

scanf("%d",&a[i]);

for(int i = 0; i < n; i++)

scanf("%d",&b[i]);

printf("%d\n",China_theory(a,b,n));

delete a;

delete b;

}

return 0;

}

2、素数相关问题

素数无穷性证明：

//欧几里德的精彩反证:

/*

假设结论不成立,即只有有限多的素数p1,p2,...,pn。

令m=1+TT(1<=i<=n)pi,即所有素数的乘积,再加一。

因为这个m比所有的素数都大,因此一定是合数。换句换说,某个素数能够整除它。

哪一个素数能整除它呢？不难知道,m不能被p1整除,因为m除以p1的余数为1。

同样的,m也不能被p2整除,因为余数也是1。事实上,当m除以任何一个pi时,余数总是1。

如果m真是合数,唯一的可能是:p1,p2,...,pn并没有包括所有的素数。

但这又和假设矛盾。

结论:素数有无穷多个！

*/

/*

素数不仅是无穷多的,而且还不算特别少。不超过x的素数大约有x/lnx个。

换句话说:大约每lnx个整数中就有一个是素数。

*/

筛法求素数：

#include

#include

#include

using namespace std;

void prime(bool *&b,int n)

{

int i,j;

b=new bool[n];

memset(b,1,sizeof(b));

b[0]=0;

for(i=2;i<=sqrt(n)+1;i++)

if(b[i-1])

for(j=1;i-1+j*i

b[i-1+j*i]=0;

for(i=0;i

if(b[i])

printf("%d\n",i+1);

}

int main()

{

//freopen("out.txt","w",stdout);

int n;

bool *b;

scanf("%d",&n);

prime(b,n);

return 0;

}

基本素数判定方法及几步改进：费马小定理：

自己的想法及错误：

米勒-勒宾测试：

整数因子分解pollard-rho算法：1、比较傻的方法（打素数表）

#include

#include

#define MAX 10000000

using namespace std;

struct power{

int x;

int y;

};

power a[MAX];

int grid(power a[],int n)

{

int i,j,k=0,p=n-1;

bool *b=new bool[n];

memset(b,1,sizeof(b));

b[0]=0;

for(i=2;i<=sqrt(double(n))+1;i++)

if(b[i-1])

for(j=1;i-1+j*i

b[i-1+j*i]=0;

for(i=0;i

if(b[i])

a[k++].x=i+1;

return k;

}

void devide(power a[],int n,int h)

{

int i;

for(;n>1;)

{

for(i=0;i

{

if(n%a[i].x==0)

{

n/=a[i].x;

a[i].y++;

break;

}

}

}

}

int main()

{

int n,h;

while(scanf("%d",&n)!=EOF)

{

for(int j=0;j

{

a[j].x=0;

a[j].y=0;

}

h=grid(a,n);

devide(a,n,h);

for(int i=0;i

{

if(a[i].y!=0)

printf("%d %d\n",a[i].x,a[i].y);

}

}

return 0;

}

2、比较优化的方法（不用打素数表）

int divide(__int64 n)

{

memset(a,0,sizeof(a));

__int64 i,k = 0;

for(i = 2; i <= sqrt(double(n)); i++)

{

if(n % i == 0)

{

while(n % i == 0)

{

n /= i;

a[k]++;

}

k++;

}

}

if(n != 1) a[k++]++;

return k;

}

N!末尾0的个数

用N一直除5，知道商为0，将各次商相加即得答案N！的素因子分解：

#include

#include

using namespace std;

void prime(bool *&b,int n){

int i,j;

b=new bool[n];

memset(b,1,sizeof(b));

b[0]=0;

for(i=2;i<=sqrt(n)+1;i++){

if(b[i-1])

for(j=1;i-1+j*i

b[i-1+j*i]=0;

}

int count_mod(int a,int n){

int sum=0,i=a;

for(;a<=n;a*=i)

sum+=n/a;

return sum;

}

int main(){

int i,n;

bool *b;

while(scanf("%d",&n)!=EOF){

prime(b,n);

for(i=0;i

if(b[i])

printf("%d %d\n",i+1,count_mod(i+1,n));

}

return 0;

}

欧拉公式：

欧拉φ函数：φ(n)是所有小于n的正整数里，和n互素的整数的个数。n是一个正整数。欧拉证明了下面这个式子：如果n的标准素因子分解式是p1^a1*p2^a2*……*pm^am，其中众pj(j=1,2,……,m)都是素数，而且两两不等。则有φ(n)=n(1-1/p1)(1-1/p2)……(1-1/pm) 利用容斥原理可以证明它。*写程序的时候，公式貌似不可以用

不借助质数表

int euler(int n){

int m = n, i;

if(n % 2 == 0) m /= 2;

while(n % 2 == 0) n /= 2;

for(i = 3; n != 1; i += 2){

if(n % i == 0) m -= m / i;

while(n % i == 0) n /= i;

}

return m;

}

借助质数表

const int maxn = 40000;

bool pf[maxn];

int p[maxn], numOfP = 0;

void init(){

memset(pf, 0, sizeof(pf));

for(int i = 2; i < maxn; i++){

if(!pf[i]){

p[numOfP++] = i;

for(int j = i * i; j < maxn; j += i)

pf[j] = true;

}

}

}

int euler1(int n){

int sum = 1;

int tmp;

for(int i = 0; i < numOfP; i++){

if(n == 1) break;

if(n % p[i] == 0){

tmp = 0;

while(n % p[i] == 0){

n /= p[i];

tmp++;

}

for(int j = 1; j < tmp; j++)

sum *= p[i];

sum *= p[i] - 1;

}

}

if(n == 1)

return sum;

else

return sum * (n - 1);

}

int euler2(int n)

{

int i;

for(i = 2; i <= sqrt((double)n); i++){

if(n % i == 0 && !pf[i]){

n /= i;

if(n % i ==0)

return i * euler2(n);

else

return (i - 1) * euler2(n);

}

}

return n-1;

}

反素数：

数论四大经典定理：

威尔逊定理:

若p为质数，则p可整除(p-1)!+1。

欧拉定理（也称费马-欧拉定理）

若n,a为正整数，且n,a互素，(a,n) = 1，则a^φ(n) ≡ 1 (mod n)

孙子定理（又称中国剩余定理）

公元前后的《孙子算经》中有“物不知数”问题：“今有物不知其数，三三数之余二，五五数之余三，七七数之余二，问物几何？”答为“23”。明朝程大位用歌谣给出了该题的解法：“三人同行七十稀，五树梅花廿一枝，七子团圆月正半，除百零五便得知。” 以现代的说法，是找出三个关键数70，21，15。解法的意思就是用70乘3除所得的余数，21乘5除所得的馀数，15乘7除所得的余数，然後总加起来，除以105的余数就是答案。

费马小定理

假如p是质数，且(a,p)=1，那么a^(p-1) ≡１（mod p）。假如p是质数，且a,p互质，那么a的(p-1)次方除以p的余数恒等于1。

Poj1730

Perfect Pth Powers

Time Limit: 1000MS Memory Limit: 10000K

Total Submissions: 10870 Accepted: 2429

Description

We say that x is a perfect square if, for some integer b, x = b2. Similarly, x is a perfect cube if, for some integer b, x = b3. More generally, x is a perfect pth power if, for some integer b, x = b p. Given an integer x you are to determine the largest p such that x is a perfect p th power.

Input

Each test case is given by a line of input containing x. The value of x will have magnitude at least 2 and be within the range of a (32-bit) int in C, C++, and Java. A line containing 0 follows the last test case.

Output

For each test case, output a line giving the largest integer p such that x is a perfect

p th power.

Sample Input

17

1073741824

25

Sample Output

1

30

2

程序

#include

#include

#include

#include

#define max 10000

using namespace std;

__int64 a[max];

inline

int divide(__int64 n)

{

memset(a,0,sizeof(a));

__int64 i,k = 0;

for(i = 2; i <= sqrt(double(n)); i++)

{

if(n % i == 0)

{

while(n % i == 0)

{

n /= i;

a[k]++;

}

k++;

}

}

if(n != 1) a[k++]++;

return k;

}

inline

int gcd(int a,int b)

{

if(b == 0) return a;

else return gcd(b,a%b);

}

int main()

{

__int64 n;

while(scanf("%I64d",&n) && n != 0)

{

__int64 h,m,g=0;

__int64 p;

if(n < 0) p = -n;

else p = n;

m = divide(p);

for(h = 0; h < m; h++)

{

g = gcd(g,a[h]);

}

if(n < 0)

while(g % 2 == 0)

g /= 2;

printf("%I64d\n",g);

}

return 0;

}