Euler Projesi 214. Soru

Totient Zincirleri

Euler totient fonksiyonu φ olsun. Yani herhangi bir n doğal sayısı için, φ(n) değeri gcd(k,n) = 1 eşitliğini sağlayan k, 1 ≤ k ≤ n, değerlerinin sayısıdır.
φ fonksiyonunu öteleyerek, her pozitif tamsayı 1 ile biten bir azalan sayı zinciri üretir.
Örneğin, 5 ile başlarsak 5,4,2,1 zinciri üretilir.
Aşağıda 4 uzunluktaki tüm zincirlerin listesi veriliyor:

5,4,2,1
7,6,2,1
8,4,2,1
9,6,2,1
10,4,2,1
12,4,2,1
14,6,2,1
18,6,2,1

Bunlardan sadece 2 tanesi bir asal sayı ile başlar, ki toplamları 12 dir.
25 uzunlukta zincirler üreten  40000000den küçük asalların toplamı kaçtır?

Cevap: 1677366278943

Java:

public class Euler214 {
   int size = 40000000;
   int totient[] = new int[size];
   Primes pr = new Primes(size);
   Euler214(){
      for( int i=0; i<size; i++ ) totient[i] = 1;
      int p = 2;
      do{
         for( int n = p; n<size; n+=p ){
            // get power of p dividing n
            int k = n/p;
            int e = 1;
            int f = p-1;
            while( k%p==0 ){
               k/=p;
               f*=p;
               e++;
            }
            
            totient[n] *= f;
         }
         p = (int)pr.getNextPrime(p);
      }while( p<size );
      
      // replace each number by its chain length
      totient[0] = totient[1] = 1; 
      for( int i=2; i<size; i++ ){
         totient[i] = totient[totient[i]]+1;
      }
      
      // sum all primes with chain length 25
      long total = 0;
      for( int i=1; i<size; i++ ){
         if( totient[i]==25 && pr.isPrime(i) )
            total += i;
      }
      System.out.println(total);
   }
}

C++:

#include <stdio.h>
#include <algorithm>
#define N 40000000 
int fac[N+1], phi[N+1], cyc[N+1];
int np, prime[N+1];
int main()
{
	int i, j, c = 0;
	__int64 s = 0;
	for(i = 0; i <= N; i++)
		fac[i] = i;
	np = 0; phi[1] = 1; cyc[1] = 1;
	for(i = 2; i <= N; i++) {
		j = i/fac[i];
		if(j % fac[i] == 0)
			phi[i] = phi[j] * fac[i];
		else
			phi[i] = phi[j] * (fac[i]-1);
		cyc[i] = cyc[phi[i]] + 1;
		if(fac[i] != i) continue;
		prime[np++] = i;
		if(cyc[i] == 25) {
			//printf("cyc(%d) = %d, phi[%d] = %d\n", i, cyc[i], i, phi[i]);
			++c; s += i;
		}
		for(j = i*i; i < 32768 && j <= N; j+=i) {
			if(fac[j] < j) continue;
			fac[j] = i;
		}
	}
	printf("\n%d %d %I64d\n",np, c, s);
	return 0;
}

Pascal:

program rox;
uses minimaple;
var phi:array[1..40000000] of longint;
    ii:longint;

procedure phig(n:longint);
var i,j,k:longint;
begin
  writeln('phi jusqu a 40 millions en cours de generation');
  phi[1]:=0;
  for i:=2 to n do
  phi[i]:=i;
  for i:=2 to n div 2 do
  if phi[i]=i then
  begin
    phi[i]:=i-1;
	k:=2;
	j:=2*i;
    while j<=n do
    begin
      phi[j]:=(phi[j] div i)* (i-1);
      k:=k+1;
	  j:=i*k;
	end;
  end;
  for i:=2 to n do
  if phi[i]=i then phi[i]:=i-1;
  writeln('phi generes en: ');
  returntime(time);
  writeln;
end;

function chaine(n:longint):longint;
var compteur:longint;
begin
  compteur:=1;
  while n<>1 do
  begin
    compteur:=compteur+1;
    n:=phi[n];
  end;
  chaine:=compteur;
end;

function nextpr(p:longint):longint;
var n:longint;
begin
  if p>=39999983 then nextpr:=40000003 else  begin
  n:=p;
  if n mod 2=0 then n:=n+1 else n:=n+2;
  while phi[n]<>n-1 do
  n:=n+2;
  if p<=1 then nextpr:=2 else
  nextpr:=n;
end;
end;

function resultat(n,taille:longint):int64;
var ch,p:longint;
    som:int64;
begin
  phig(n);
  som:=0;
  p:=15;
  while p<=n do
  begin
    ch:=chaine(p);
    if ch=taille then som:=som+p;
    p:=nextpr(p);
  end;
  resultat:=som;
end;

begin
  writeln(resultat(40000000,25));
end.

Mathematica:

m = 0;
For[t = 1, t <= PrimePi[40000000], ++t, {
  n = Prime[t];
  k = 0;
  While[EulerPhi[n] != 1, {x = EulerPhi[n], n = x, k = k + 1}]
   b = k + 2;
  If[b == 25, m = m + Prime[t]]}]
m

Maple:

> with(numtheory):

> nbiter:=proc(p)
> local k,compteur;
> k:=p;
> compteur:=1;
> while k<>1 do 
>   k:=phi(k);
>   compteur:=compteur+1;
> od;
> compteur;
> end:
> 
> cherche:=proc(n)
> local p,som;
> p:=5;
> som:=0;
> while p<n do
>   if nbiter(p)=25 then som:=som+p fi;
>   p:=nextprime(p)
> od;
> som
> end:
> cherche(40000000);

Matlab:

% SL: sequence length for all numbers
% SL_p sequence length just for the primes
% p: vector of primes till 4e7
for np = 1:length(p)
    curv = p(np)-1;
    hist = p(np);
    while SL(curv) == 0
        hist = [hist,curv];
        curv = totient(curv);
    end
    SL(hist) = SL(curv) + (length(hist):-1:1);
    SL_p(np) = SL(p(np));
end
TheSum = sum(p(SL_p==25));