阶乘的整除性
使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
#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
#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;
}
--修改前
Timer 3.01 Copyright (c) 2002-2003 Igor Pavlov 2003-07-10
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
Timer 3.01 Copyright (c) 2002-2003 Igor Pavlov 2003-07-10
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
Timer 3.01 Copyright (c) 2002-2003 Igor Pavlov 2003-07-10
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
#include
#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;
}