Sadeleştirilemeyen bir kesre 'dirençli' kesir adı verilir. Bir d doğal sayısı için, paydası d olan tüm dirençli basit kesirlerin sayısının tamamına oranı R(d) olsun; örneğin R(12)=4/11.
Bu durumda d>1 sayısının direnci $φ(d)/(d-1)$ ile tanımlanır. Burada φ, Euler totient fonksiyonudur.
Ayrıca bir n>1 sayısının eşdirencini $C(n)=[n-φ(n)]/(n-1)$ olarak tanımlanır. Bir p asal sayısının eşdirenci $C(p)=1/(p-1)$ olur.
$1<n \le 2\times 10^{11}$ için C(n) bir birim kesir olacak şekilde tüm bileşik n tam sayılarının toplamını bulunuz.
Cevap: 288084712410001
Fortran:
program p245
use MillerRabin
implicit none
integer(LONG),parameter:: LIMIT=10_LONG**11
integer(LONG),parameter:: MAXPRIME=4*10**5
integer(LONG):: list(10)
logical:: prime(MAXPRIME)
integer(LONG):: p,s=0
call init_primes
do p=3,MAXPRIME,2
if (.not. prime(p)) cycle
list(1)=p
call factor(2,p,p,p-1)
end do
write(*,*) 's=',s
contains
subroutine init_primes
implicit none
integer:: p
prime(:)=.true.
prime(1)=.false.
prime(2:MAXPRIME:2)=.false. ! prime factor 2 does not be considered
do p=3,MAXPRIME,2
if (p*p>MAXPRIME) exit
if (.not. prime(p)) cycle
prime(p*p:MAXPRIME:p)=.false.
end do
end subroutine init_primes
logical function is_prime2(p)
implicit none
integer(LONG):: p
if (p<=MAXPRIME) then
is_prime2=prime(p) ! use sieved list
else
is_prime2=is_prime(p) ! use Miller-Rabin
end if
end function is_prime2
recursive subroutine factor(i,p0,PP,QQ)
implicit none
integer,intent(in):: i
integer(LONG),intent(in):: p0,PP,QQ ! PP=prod(p_i), QQ=prod(p_i-1)
integer:: ii
integer(LONG):: PP1,QQ1
integer(LONG):: p,pmin,pmax,r,rmin,rmax
! one more factor
pmin=p0+2
pmax=LIMIT/PP
rmin=(PP*pmin-2)/(PP*pmin-QQ*(pmin-1))+1
rmax=(PP*pmax-1)/(PP*pmax-QQ*(pmax-1))
do r=rmin,rmax
if (mod(QQ*r+1,PP - (PP-QQ)*r)/=0) cycle
p=(QQ*r+1)/(PP - (PP-QQ)*r)
if (.not. is_prime2(p)) cycle
PP1=PP*p
list(i)=p
s=s+PP1
write(*,*) PP1,'=',(list(ii),ii=1,i)
end do
! more than one factor
do p=p0+2,MAXPRIME,2
if (.not. prime(p)) cycle
list(i)=p
PP1=PP*p
if (PP1*(p+2)>LIMIT) exit
QQ1=QQ*(p-1)
call factor(i+1,p,PP1,QQ1)
end do
end subroutine factor
end program p245
C/C++:
#define LIMIT 100000000000ULL
#define PRIME_LIMIT 22000000
typedef unsigned long long int int_64;
int_64 p245(void) {
int i = 0;
int p;
int sqrt_limit = (int)sqrt(LIMIT);
int_64 result = 0;
primes = sieve(PRIME_LIMIT);
for ( i = 0; (p = primes[i]) < sqrt_limit ; i++) {
result += search(p, p-1, i+1, p-1);
}
return result;
}
int_64 search(int_64 n, int_64 phi, int i, int_64 m) {
int_64 m1, r;
int_64 n1, phi1;
int_64 result;
int p, q;
result = 0;
// check the numbers at this level, (and go deeper as required)
while (1) {
p = primes[i++];
n1 = n*p;
if ( n1 > LIMIT ) {
break;
}
phi1 = phi*(p-1);
r = (n1-1) % (n1-phi1);
if ( r == 0 ) {
result += n1;
}
m1 = (n1-1) / (n1-phi1);
if ( m1 >= m ) {
i--;
break;
}
result += search(n1, phi1, i, m1);
}
// check if we need to go a level deeper, although this level is exhuasted
while (1) {
p = primes[i++];
q = primes[i];
n1 = n*p*q;
if ( n1 > LIMIT ) {
break;
}
phi1 = phi*(p-1)*(q-1);
m1 = (n1-1)/ (n1-phi1);
if (m1 >= m) {
break;
}
n1 = n*p;
phi1 = phi*(p-1);
result += search(n1, phi1, i, (n1-1)/(n1-phi1));
}
return result;
}
Python:
def rec2(pi,X,Y,st=[]):
global total;
if(len(st)>1 and (X-1) % (X-Y)==0): total += X;
for i in xrange(pi,len(primes)):
p=primes[i];
if(X*p>mx): break;
rec2(i+1,X*p,Y*(p-1),st+[p]);
#================================
from globals import getprimes;
from time import time;
import psyco;
psyco.full();
mx=10**11;
total=70945825053;
start=time();
primes=getprimes(22*10**6);
for pi in xrange(11,len(primes)):
p=primes[pi];
rec2(pi+1,p,p-1,[p]);
print "Total = %d"%total;
print "Total time taken: %fs"%(time()-start);
Delphi:
procedure TForm1.Button11Click(Sender: TObject);
const limit=trunc(1e11);
var s:array of int64;
var pf:array of array of integer;
i,j:cardinal;k,m,p:int64;
slim:integer; count:integer;sum:int64;
divisors:array[0..1000] of cardinal;
start,stop,freq:int64;
procedure process(p:integer);
var top,curtop:integer;q,q2,f2,h:int64;i,j,k,mul,f:integer;
begin
q:=int64(p)*p-p+1;
q2:=0;
if high(pf[p])=-1 then
begin
q:=q-(p-1);
if (p*q<=limit) and primesieve.isprime(q) then
begin
inc(count);inc(sum,p*q);
end;
exit;
end;
divisors[0]:=1;
top:=0;curtop:=0;
q2:=q;
for i:=0 to high(pf[p]) do
begin
f:=pf[p,i];
q2:=trunc(q2/f);
mul:=1;
while q2 mod f=0 do begin q2:=q2 div f; inc(mul) end;
f2:=f;
for j:=1 to mul do
begin
for k:=0 to top do
begin
h:=divisors[k]*f2;if h
p) and (p*q2<=limit) and primesieve.isprime(q2) then
begin
inc(count);
inc(sum,p*q2);
end;
end;
end;
begin
queryperformancecounter(start);
queryperformancefrequency(freq);
primesieve.init_sieve(trunc(1.1*power(limit,2/3)));
slim:=trunc(sqrt(limit+0.0));
setlength(s,slim+100);
setlength(pf,slim+100);
for i:=3 to slim do s[i]:=int64(i)*i-i+1;
i:=5;
while i<=slim do
begin
if i shr 1>3 then
begin
setlength(pf[i],1);
pf[i,0]:=3;
end;
s[i]:=trunc(s[i]/3);
while s[i]-trunc(s[i]/3)*3 =0 do s[i]:=trunc(s[i]/3);
inc(i,3);
end;
for i:=3 to slim do
begin
if s[i]>1 then
begin
p:=s[i];
if (i shr 1>p) and primesieve.isprime(i) then
begin
setlength(pf[i],length(pf[i])+1);
pf[i,high(pf[i])]:=p;
end;
//s[i]:=1;
k:=p-i+1;
while k<=slim do
begin
if (k shr 1>p) and primesieve.isprime(k) then
begin
setlength(pf[k],length(pf[k])+1);
pf[k,high(pf[k])]:=p;
end;
s[k]:=trunc(s[k]/p);
while s[k]=trunc(s[k]/p)*p do s[k]:=trunc(s[k]/p);
k:=k+p
end;
k:=p+i;
while k<=slim do
begin
if (k shr 1>p) and primesieve.isprime(k) then
begin
setlength(pf[k],length(pf[k])+1);
pf[k,high(pf[k])]:=p;
end;
s[k]:=trunc(s[k]/p);
while s[k]=trunc(s[k]/p)*p do s[k]:=trunc(s[k]/p);
k:=k+p
end;
end;
end;
i:=primesieve.nextprime(2);
sum:=0;count:=0;
while i<=slim do
begin
process(i);
i:=primesieve.nextprime;
end;
queryperformancecounter(stop);
memo1.lines.add(inttostr(count)+' '+inttostr(sum));
memo1.lines.add(floattostr((stop-start)/freq*1000));
end;