Euler Projesi 272. Soru

Modüler Küpler, kısım 2

Bir pozitif n sayısı için $1<x<n$ ve $x^3\equiv 1$ mod n olmak üzere x tam sayılarının sayısı C(n) ile tanımlansın.

n=91 iken x için 8 olası değer mevcut:  9, 16, 22, 29, 53, 74, 79, 81.
Yani, C(91)=8.

C(n)=242 olan tüm $n\leq 10^{11}$ pozitif sayılarının toplamını bulunuz.
Cevap: 8495585919506151122
C/C++: #include <stdio.h> #include <stdlib.h> #include <time.h> int main() { int *isprime,sec; long long int sum=0,*pr,interval,i,j,si,B,N,pos[32],ct,depth,backtrack; long long int mul[]={0,0,0,0,0,1,-6,21,-56,126,-252,462,-792,1287,-2002,3003}; scanf("%lld",&B); interval=B/15561+1; pr=(long long int*)(malloc)(interval*sizeof(long long int)); isprime=(int*)(malloc)(interval*sizeof(int)); for(i=0;i<interval;i++) isprime[i]=(i>1); for(i=0;i*i<=interval;i++) { if(isprime[i]) { for(j=i*i;j<interval;j+=i) isprime[j]=0; } } si=4; pr[1]=7; pr[2]=9; pr[3]=13; pr[4]=19; for(i=20;i<interval;i++) if(isprime[i]&&(i%3==1)) si++,pr[si]=i; printf("primes=%lld\n",si); sec=time(NULL); depth=1; pos[1]=0; while(depth>0) { backtrack=0; pos[depth]++; if(pos[depth]>si) backtrack=1; else { N=1; for(i=1;i<=depth;i++) N*=pr[pos[i]]; if(N>B) backtrack=1; else { ct=B/N; sum+=mul[depth]*ct*(ct+1)/2*N; depth++; pos[depth]=pos[depth-1]; } } if(backtrack) depth--; } printf("sum=%lld,time=%ld sec.\n",sum,time(NULL)-sec); } Haskell: import ONeillPrimes (primes) gen1 :: Integral a => a -> [a] -> Int -> [a] gen1 _ _ 0 = [1] gen1 u (p:ps) (k+1) | u < p^(k+1) = [] | otherwise = concat [map (q*) $ gen1 (div u q) ps k| q <- takeWhile (<= u).iterate (p*) $ p] ++ gen1 u ps (k+1) gen2 :: Integral a => a -> [a] -> [a] gen2 u (p:ps) | u < p = [1] | otherwise = concat [map (q*) $ gen2 (div u q) ps| q <- takeWhile (<= u).iterate (p*) $ 1] p270 :: Integral a => a -> Int -> a p270 u k = sum.map count $ gs1 ++ gs2 where ps1 = filter ((1==).(`mod` 3)) primes ps2 = filter ((2==).(`mod` 3)) primes gs1 = gen1 u ps1 k gs2 = concat [map (q*) $ gen1 (div u q) ps1 (k-1)| q <- takeWhile (<= div u q).iterate (3*) $ 9] q = product.take (k-2) $ ps1 count g | mod g 9 == 0 = g * sum hs | otherwise = g * (sum hs + 3 * sum (filter (<= (div u $ 3*g)) hs)) where hs = gen2 (div u g) ps2 main :: IO () main = print.p270 (10^11) $ 5 Python: import math n = 10 ** 11 u = [0] * 5 + [1, -6, 21, -56] c = 0 nx = (n + 1) // (9 * 7 * 13 * 19 - 1) def v(t, i, a, r): if a == r: global c m = n // t c += u[r] * t * (m * (m + 1)) // 2 return elif a < r: q = xx[i] while t * q ** (r - a) <= n and i < len(xx) - 1: v(t * q, i + 1, a + 1, r) i += 1 q = xx[i] d = [True] * (nx + 1) for x in range(2, int(math.sqrt(len(d))) + 3): if d[x]: for y in range(2 * x, len(d), x): d[y] = False xx = [x for x in range(3, len(d)) if x == 9 or d[x] and x % 3 == 1] for i in range(5, 9): v(1, 0, 0, i) print(c) Delphi: program Euler272; uses sysutils, windows, funcs1; var primeTrivial, PrimeNonTrivial: array[1..4000000] of integer; noTrivial, noNonTrivial: integer; const maxProduct: int64 = 100000000000; const maxi = 440087; var sum: extended; //7*9*13*19 15561 prod1 15561 // index 440087 nthPrime(index) 6426349 // prod1*nthPrime(index) 100000416789 const index=440087; procedure doProdsTrivial(start: integer; Prod: extended; has3: boolean); var i, p: integer; newProd, newProd3: extended; begin if not has3 then begin p:= 3; newProd3:= Prod * p; if newProd3<=maxProduct then begin sum:= sum+newProd3; doProdsTrivial(start, newProd3, true); // true: only one single 3 doProdsTrivial(start, Prod, true); // true: only one single 3 exit; end; end (* exact one 3 *); for i:= start to noTrivial do begin p:= PrimeTrivial[i]; if p>207300 then exit; newProd:= Prod * p; if newProd>maxProduct then exit; // a new number: sum:= sum+newProd; doProdsTrivial(i, newProd, has3); end (* for i *); end; procedure doProdsNonTrivial(start: integer; aktNumber: integer; Prod: extended); var i, p: integer; newProd: extended; has3: boolean; // doProdsNonTrivial: multiply 5 non-triviale primes // - even in higher power - and then with as much as possible // trivial primes. - beware of the 3 - begin if aktNumber=5 then begin sum:= sum+Prod; has3:= round(Prod) mod 3 = 0; doProdsTrivial(1, Prod, has3); p:= PrimeNonTrivial[start]; newProd:= Prod * p; if newProd>maxProduct then exit; doProdsNonTrivial(start, aktNumber, newProd); exit; end; for i:= start to noNonTrivial do begin p:= PrimeNonTrivial[i]; if aktNumber=0 then if p>160 then exit; if aktNumber=1 then if p>563 then exit; if aktNumber=2 then if p>4642 then exit; if aktNumber=3 then if p>316228 then exit; newProd:= Prod * p; if aktNumber = 0 then if p=3 then newProd:= newProd*3; if newProd>maxProduct then exit; if (aktNumber=0) and (p=3) then begin doProdsNonTrivial(i, aktNumber+1, newProd); end else if i=start then doProdsNonTrivial(i, aktNumber, newProd) else doProdsNonTrivial(i, aktNumber+1, newProd); end (* for i *); end; var theTime: longword; var idx, i, p: integer; s: string; begin theTime:= GetTickCount; fillPrimesNumbered(7000000, idx); noTrivial:= 0; noNonTrivial:= 0; for i:= 1 to index do begin p:= nthPrime(i); if p mod 3 = 2 then begin inc(noTrivial); PrimeTrivial[noTrivial]:= p; end else begin inc(noNonTrivial); PrimeNonTrivial[noNonTrivial]:= p; end; end; sum:= 0; doProdsNonTrivial(1, 0, 1); str(round(sum):18, s); writeln('sum ',s); theTime:= GetTickCount-theTime; writeln('time ',theTime); end.