19 Temmuz 2018

Euler Projesi 245. Soru

Eşdirenç

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;