Description
Today's math class, crash children learn least common multiple (Least Common multiple). For two positive integers a and B,LCM (A, b) represent the smallest positive integers that can be divisible by A and b at the same time. For example, the LCM (6, 8) = 24. Back home, crash still thinking about class school, in order to study least common multiple, he drew a n*m form. A number is written in each lattice, where the lattice in column J of Row I is written in the LCM (i, J). A 4*5 table is as follows:
1 2 3) 4 5
2 2 6) 4 10
3 6 3) 12 15
4 4 12 4 20 Looking at this form, crash thought of a lot of questions to think about. But the question he most wants to solve is a very simple question: how much is there in this table? When N and M are big, crash is helpless, so he finds the smart you use the program to help him solve the problem. Since the end result can be large, crash only wants to know the values of all the numbers in the table and mod 20101009.
Input
The first line of input contains two positive integers, representing N and M., respectively.
Output
Outputs a positive integer that represents the number of values in the table and mod 20101009.
∑∑LCM (i,j), I=1..N,J=1..M
=∑∑I*J/GCD (i,j), I=1..N,J=1..M
=∑∑∑[GCD (i,j) =d]*i*j/d, D=1..N,I=1..N,J=1..M (enumeration GCD's value D)
=∑d* (∑∑[GCD (i,j) =1]*i*j), d=1..n,i=1..n/d,j=1..m/d
=∑d* (∑a*a*mu[a]* (∑∑i*j)), d=1..n,a=1..n/d,i=1..n/d/a,j=1..m/d/a (from [N=1]=∑mu[d], d|n available)
The value of a is only O (n0.5) and can only be calculated once for each value
The value of each A, (n/d/a,m/d/a) is also O (n0.5), which can be optimized in the same way
∑∑i*j can be calculated using the arithmetic progression summation formula O (1)
#include <cstdio>typedefLong Longi64;Const intn=1e7,p=20101009;intps[n/5],pp=0, mu[n+5];i64 m2[n+5];BOOLisnp[n+5];intn,m;inline i64 F (i64 a,i64 b) {return(A * (A +1) >>1)%p* (b* (b +1) >>1) (%P))%P;} InlineintMinintAintb) { returnA<b?a:b;}intMain () {scanf ("%d%d",&n,&m); if(n>m) n^=m,m^=n,n^=m; mu[1]=1; for(intI=2; i<=n;i++){ if(!isnp[i]) ps[pp++]=i,mu[i]=-1; for(intj=0,k;j<pp&& (K=i*ps[j]) <=n;j++) {Isnp[k]=1; if(I%ps[j]) mu[k]=-Mu[i]; Else Break; } } for(intI=1; i<=n;i++) {M2[i]= (mu[i]*1ll*i*i+m2[i-1])%Q; if(m2[i]<0) m2[i]+=Q; } i64 ans=0; for(intD=1, d2;d<=n;d=d2+1) {D2=min (n/(n/d), m/(m/d)); intx=n/d,y=m/D; i64 s=0; for(intA=1, b;a<=x;a=b+1) {b=min (x/(x/a), y/(y/a)); S+ = (m2[b]-m2[a-1]) *f (x/a,y/a)%P; } s%=P; Ans= (ans+ (1ll* (D+D2) * (d2-d+1) >>1)%p*s%p)%Q; } if(ans<0) ans+=Q; printf ("%lld\n", ans); return 0;}
bzoj2154 Crash's digital table