C语言old阶乘,欧拉计划549题:阶乘的整除性

最新推荐文章于 2021-05-21 09:59:39 发布

豆几宝

最新推荐文章于 2021-05-21 09:59:39 发布

阅读量572

点赞数 2

文章标签： C语言old阶乘

阶乘的整除性

使m!能被10整除的最小的m为m=5。

使m!能被25整除的最小的m为m=10。

记使m!能被n整除的最小的m为s(n)。

因此s(10)=5，s(25)=10。

记S(n)为∑s(i)，其中2 ≤ i ≤ n。

已知S(100)=2012。

求S(108)。

用SQL先看一下s(n)的组成，显然s(质数)=质数本身，因数越多的数，s(n)越小，如果一个数的最大因数是质数k,s(n)=k，如15=3*5，则s(15)=5，如果一个数的最大因数是合数，暂时没发现规律。

with t as(select level l from dual connect by level<=25)

,fact(m,a) as(select 1,1 from dual

union all

select l,a*l from t,fact where l=m+1)

select l n,min(m) sn from fact,t where mod(a,l)=0 group by l order by 1;

N SN

----- ----------

1 1

2 2

3 3

4 4

5 5

6 3

7 7

8 4

9 6

10 5

11 11

N SN

----- ----------

12 4

13 13

14 7

15 5

16 6

17 17

18 6

19 19

20 5

21 7

22 11

N SN

----- ----------

23 23

24 4

25 10

因为阶乘很容易就变得很大，肯定不能直接计算，而是把它当作一个包含所有1~m因数的数。

还是求n的质因数积，比如48=24*3，再求m!的质因数积，根据幂次去找最小的，而m!的质因数积可以通过n的质因数积推导。

求质因数积还是类似筛法比较好，否则有很多除法和取模操作。

先提供一个简易的答案,p549.h中包含了一个质数表。这个实现显然不可能解决S(108).

#include

typedef std::map mapii;

#include "p549.h"

mapii mp[10000]; //save n

mapii mp2[10000]; //save m!

int get__prime_fact_prod(int n)

{

int oldn=n;

for(int j=0; n!=1 ; j++)

{

int t=prime[j];

if(n%t==0)

{

mp[oldn][t]++;

n/=t;

j--;

}

int get_m_fact_prod(int m) //2!->m!'s prime_fact_prod

{

mp2[2]=mp[2];

for(int i=3; i<=m; i++)

{

mp2[i]=mp2[i-1];

mapii:: iterator it;

for(it=mp[i].begin(); it!=mp[i].end(); it++)

mp2[i][it->first]+=it->second;

}

int print_map(mapii mpa[],int b,int e)

{

for(int i=b; i<=e; i++)

{

printf("%d: ",i);

mapii:: iterator it;

for(it=mpa[i].begin(); it!=mpa[i].end(); it++)

printf("%d^%d ",it->first,it->second);

puts("");

}

int match_fact_prod(mapii mpa[],mapii mpa2[],int n)

{

int match=0;

for(int i=2; match==0; i++)

{

mapii:: iterator it;

for(it=mpa[n].begin(); it!=mpa[n].end(); it++)

{

if(mpa2[i][it->first]second)

break;

}

if (it==mpa[n].end())

{

match=i;

break;

}

return match;

}

int main()

{

for(int i=2; i<=100; i++)

get__prime_fact_prod(i);

//print_map(mp,2,100);

get_m_fact_prod(97); //max prime in 100

//print_map(mp2,2,97);

int sum=0;

for(int i=2; i<=100; i++)

{

sum+=match_fact_prod(mp,mp2,i);

//printf("%d %d\n",i,x);

}

printf("%d\n",sum);

return 0;

}

上述程序很慢，现在逐步来修改它，看在合理的时间内，能算到多少。

前面指出质数的s(n)即本身，因此可以直接加总，而不必调用函数。

我们用一个char类型数组来保存质数标记。并在主函数作了修改，其他不动。

#include

typedef std::map mapii;

#include "p549.h"

#define NN 10000

mapii mp[NN+1]; //save n

mapii mp2[NN+1]; //save m!

char a[NN+1]={0}; //prime flag

int main()

{

for(int i=2; i<1225; i++) //set prime pos ,total 1228 prime number in p594.h

a[prime[i]]=1;

for(int i=2; i<=NN; i++)

get__prime_fact_prod(i);

//print_map(mp,2,100);

get_m_fact_prod(NN);

//print_map(mp2,2,97);

int sum=0;

for(int i=2; i<=NN; i++)

{

if(a[i]==1)

sum+=i;

else

sum+=match_fact_prod(mp,mp2,i);

//printf("%d %d\n",i,x);

}

printf("%d\n",sum);

return 0;

}

--修改前

10125843

Kernel Time = 0.156 = 00:00:00.156 = 3%

User Time = 4.820 = 00:00:04.820 = 96%

Process Time = 4.976 = 00:00:04.976 = 99%

Global Time = 4.977 = 00:00:04.977 = 100%

--修改后

D:\>a\timer a

10125843

Kernel Time = 0.187 = 00:00:00.187 = 10%

User Time = 1.638 = 00:00:01.638 = 89%

Process Time = 1.825 = 00:00:01.825 = 100%

Global Time = 1.825 = 00:00:01.825 = 100%

这个程序，clang++编译的结果和g++差不多。

时间都花在哪里了，下面用gnomon测试一下,可见，主要是后2步比较耗时。

D:\>a |gnomon -i

0.0112s Start

0.0010s get_n_fact

0.2938s get_m!_fact

0.0003s 10125843

0.3339s all done

Total 0.6444s

实际上，没必要计算出所有m!的质因数积，因为质数s(n)就等于本身，用不着它。而合数，最大因数不过n/2，所以s(n)不会超过(n/2)!。所以，把第2步计算的范围减少1半。get_m_fact_prod(NN/2+2);

D:\>a\timer a

Start

get_n_fact

get_m!_fact

10125843

all done

Kernel Time = 0.031 = 00:00:00.031 = 2%

User Time = 1.388 = 00:00:01.388 = 97%

Process Time = 1.419 = 00:00:01.419 = 99%

Global Time = 1.420 = 00:00:01.420 = 100%

考虑质数m的倍数j，如果这个倍数j小于质数本身，那么质数m的阶乘中必然已经包含此倍数j，因此也不用计算。直接就知道所有的s(m*j)(其中j

int main()

{

int t=clock();

P("Start")

for(int i=2; i<1230; i++) //set prime pos ,total 1228 prime number in p594.h

a[prime[i]]=1;

for(int i=2; i<=NN; i++)

get__prime_fact_prod(i);

P("get_n_fact")

//print_map(mp,2,100);

get_m_fact_prod(NN/2+2); //opt

P("get_m!_fact")

//print_map(mp2,2,97);

int sum=0;

//fill prime *2

for(int i=2; i<=NN/2; i++)

{

if(a[i]==1)

{

sum+=(i+(i==2)); // only 2*2 need +1 , other prime >2 , prime *2 ->m! too

a[i*2]=2;

}

//fill prime *3

for(int i=2; i<=NN/3; i++)

{

if(a[i]==1)

{

sum+=(i+(i==3)); // only 3*3 need +1 , other prime >3 , prime *3 ->m! too ,2*3 already filled

a[i*3]=3;

}

for(int i=2; i<=NN; i++)

{

if(a[i]==1)

sum+=i;

else if(a[i]==0) //other flag already add to sum

sum+=match_fact_prod(mp,mp2,i);

}

printf("%d\n",sum);

P("all done")

return 0;

}

--原始

1591ms

--添加处理2

826ms

--添加2，3处理

577ms

如果把上述2,3的处理再继续扩大到其他倍数,速度提高很多，可惜算错了。

for(int j=2;j<=10;j++) //if NN=1e4, j less than sqrt(NN)

for(int i=j; i<=NN/j; i++)

{

if(a[i]==1)

{

sum+=(i+(i==j)); // only 2*2 need +1 , other prime >2 , prime *2 ->m! too

a[i*j]=j;

}

--j to 100

D:\p549>g++ p549f.cpp -O2

D:\p549>a

Start

0ms

get_n_fact

15ms

get_m!_fact

200ms

10124811 2660

all done

210ms

--j to 10

D:\p549>g++ p549f.cpp -O2

D:\p549>a

Start

0ms

get_n_fact

13ms

get_m!_fact

198ms

10125833 5945

all done

268ms

能正确计算S(100)的改进程序，思路见回帖。

#include

#define NN 1000000

#define P(X) printf("%s %ldms\n",X,clock()-t);fflush(stdout);

char a[NN+1]= {0}; //prime flag,1 is primer number

int prime[NN/4+1]= {0}; //min(NN)=100,NN is 100^n (n=1->4)

int dat[8][32]= {0}; //save f(t,e)'s return value

int filldat()

{

for(int i=2; i<=7; i++) // i is base

{

//puts("");

int sum_e[32]= {0};

int sumover=0;

int j=0;

if(a[i]==1) //only prime

{

for(j=1; pow(i,sumover)<=NN; j++) //i*j is i times 1 2 3..

{

int k=i*j;

int n=0;

for (; k%i==0; n++)

k/=i;

sumover+=n;

sum_e[j]=sumover;

}

for(j--; j>0; j--)

for(int k=sum_e[j]; k>sum_e[j-1]; k--)

dat[i][k]=j;

}

int f(int t,int i)

{

if(i==1)

return 1;

if(a[t]==1 && a[t]<=7)

return dat[t][i];

else

return i;

}

void sieve()

{

a[2]=1;

for(int i=3; i<=NN; i++)

{

a[i++]=1; //nod

a[i]=0; //even

}

for(int i=3; i<=sqrt(NN); i++)

if(a[i]==1)

for(int j=i*i; j<=NN; j+=(2*i))

a[j]=0;

int d=0;

for(int i=2; i<=NN; i++)

if(a[i]==1)

{

prime[d++]=i;

}//if(i<400)printf("%d ",i);

printf("size of sieve is %d\n",d);

}

int get_prime_fact_prod(int n)

{

int tmpn=n;

int e[10000+1]= {0}; //max prime =sqrt(NN)

for(int j=0; tmpn!=1 ; j++)

{

int t=prime[j];

if(tmpn%t==0)

{

e[j]++; //exp

tmpn/=t;

j--;

}

//all data got, begin to match

int m=0;

for(int j=0; j<10 ; j++)

{

int x=0;

if(e[j]>0)

x=f(prime[j],e[j])*prime[j];

if (m

m=x;

}

//printf("%d=%d\n",n,m);

return m;

}

int main()

{

int t=clock();

P("Start")

sieve();

filldat();

P("Sieve")

long long sum=0;

//fill prime *<=prime

for(int j=2; j<=sqrt(NN); j++) //if NN=1e4, j less than sqrt(NN)

{

if(a[j]==1)

{

a[j*j]=j*2;

sum+=j*2;

}

for(int i=j+1; i<=NN/j; i++)

{

if(a[i]==1)

{

sum+=i;

a[i*j]=i;

}

P("Prime fact biggest")

int c=0;

for(int i=2; i<=NN; i++)

if(a[i]==0)

{

int temp=get_prime_fact_prod(i);

sum+=temp;

a[i]=temp;

c++;

}

P("get_n_fact")

for(int i=2; i<=NN; i++)

{

if(a[i]==1)

{

sum+=i;

a[i]=i;

}

printf("S(%d)=%lld call %d\n",NN,sum,c);

P("all done")

return 0;

}

豆几宝

关注

2
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
C语言old阶乘,欧拉计划549题:阶乘的整除性

阶乘的整除性使m!能被10整除的最小的m为m=5。使m!能被25整除的最小的m为m=10。记使m!能被n整除的最小的m为s(n)。因此s(10)=5，s(25)=10。记S(n)为∑s(i)，其中2 ≤ i ≤ n。已知S(100)=2012。求S(108)。用SQL先看一下s(n)的组成，显然s(质数)=质数本身，因数越多的数，s(n)越小，如果一个数的最大因数是质数k,s(n)=k，如15=3...
复制链接

扫一扫