Superszybka funkcja totient

22

Cel jest prosty: oblicz funkcję sumaryczną dla jak największej liczby liczb w 10 sekund i zsumuj liczby.

Musisz wydrukować wynik na końcu i faktycznie go obliczyć. Żadna automatyczna funkcja totient nie jest dozwolona, ​​ale biblioteki bignum są. Musisz zacząć od 1 i policzyć kolejno liczby całkowite. Teraz nie wolno pominąć numery.

Twój wynik to ile liczb twój program może obliczyć na twoim komputerze / ile mój program może obliczyć na twoim komputerze . Mój kod jest prostym programem w C ++ (wyłączone optymalizacje), mam nadzieję, że możesz go uruchomić.

Ważne właściwości, których możesz użyć!

  • Jeśli gcd(m,n) = 1, phi(mn) = phi(m) * phi(n)
  • jeśli pjest liczbą pierwszą phi(p) = p - 1(dla p < 10^20)
  • jeśli njest nawetphi(2n) = 2 phi(n)
  • inne wymienione w pierwszym linku

Mój kod

#include <iostream>
using namespace std;

int gcd(int a, int b)
{
    while (b != 0)
    {
        int c = a % b;
        a = b;
        b = c;
    }
    return a;
}

int phi(int n)
{
    int x = 0;
    for (int i=1; i<=n; i++)
    {
        if (gcd(n, i) == 1)
            x++;
    }
    return x;
}

int main()
{
    unsigned int sum = 0;
    for (int i=1; i<19000; i++) // Change this so it runs in 10 seconds
    {
        sum += phi(i);
    }
        cout << sum << endl;
        return 0;
}
qwr
źródło
2
Może warto dodać, że liczby wejściowe powinny być kolejnymi liczbami całkowitymi. W przeciwnym razie może mnie pokusić się o obliczenie funkcji sumarycznej tylko dla potęg 2.
Howard
Czy mogę to zrobić 1, 3, 5, 2, 4?
Leaky Nun

Odpowiedzi:

14

Nimrod: ~ 38 667 (580 000 000/15 000)

W tej odpowiedzi zastosowano dość proste podejście. W kodzie zastosowano proste sito liczby pierwszej, które przechowuje liczbę pierwszą najmniejszej mocy liczby pierwszej w każdym gnieździe dla liczb zespolonych (zero dla liczb pierwszych), a następnie używa programowania dynamicznego do skonstruowania funkcji totient w tym samym zakresie, a następnie sumuje wyniki. Program spędza praktycznie cały czas na budowie sita, a następnie oblicza funkcję totient w ułamku czasu. Wygląda na to, że sprowadza się to do zbudowania wydajnego sita (z lekkim zwrotem, że trzeba również wyodrębnić główny czynnik dla liczb zespolonych z wyniku i utrzymać zużycie pamięci na rozsądnym poziomie).

Aktualizacja: Poprawiona wydajność poprzez zmniejszenie zajmowanego miejsca w pamięci i poprawę zachowania pamięci podręcznej. Można wycisnąć o 5% -10% większą wydajność, ale wzrost złożoności kodu nie jest tego wart. Ostatecznie algorytm ten przede wszystkim wykorzystuje wąskie gardło procesora von Neumanna i istnieje bardzo niewiele poprawek algorytmicznych, które mogą to obejść.

Zaktualizowałem również dzielnik wydajności, ponieważ kod C ++ nie miał być kompilowany przy wszystkich optymalizacjach i nikt tego nie zrobił. :)

Aktualizacja 2: Zoptymalizowane działanie sita dla lepszego dostępu do pamięci. Teraz obsługa dużych liczb pierwszych za pomocą memcpy () (przyspieszenie ~ 5%) i pomijanie wielokrotności 2, 3 i 5 podczas przesiewania większych liczb pierwszych (przyspieszenie ~ 10%).

Kod C ++: 9,9 sekundy (z g ++ 4.9)

Kod Nimrod: 9,9 sekundy (z -d: release, backend gcc 4.9)

proc handleSmallPrimes(sieve: var openarray[int32], m: int) =
  # Small primes are handled as a special case through what is ideally
  # the system's highly optimized memcpy() routine.
  let k = 2*3*5*7*11*13*17
  var sp = newSeq[int32](k div 2)
  for i in [3,5,7,11,13,17]:
    for j in countup(i, k, 2*i):
      sp[j div 2] = int32(i)
  for i in countup(0, sieve.high, len(sp)):
    if i + len(sp) <= len(sieve):
      copyMem(addr(sieve[i]), addr(sp[0]), sizeof(int32)*len(sp))
    else:
      copyMem(addr(sieve[i]), addr(sp[0]), sizeof(int32)*(len(sieve)-i))
  # Fixing up the numbers for values that are actually prime.
  for i in [3,5,7,11,13,17]:
    sieve[i div 2] = 0

proc constructSieve(m: int): seq[int32] =
  result = newSeq[int32](m div 2 + 1)
  handleSmallPrimes(result, m)
  var i = 19
  # Having handled small primes, we only consider candidates for
  # composite numbers that are relatively prime with 31. This cuts
  # their number almost in half.
  let steps = [ 1, 7, 11, 13, 17, 19, 23, 29, 31 ]
  var isteps: array[8, int]
  while i * i <= m:
    if result[i div 2] == 0:
      for j in 0..7: isteps[j] = i*(steps[j+1]-steps[j])
      var k = 1 # second entry in "steps mod 30" list.
      var j = 7*i
      while j <= m:
        result[j div 2] = int32(i)
        j += isteps[k]
        k = (k + 1) and 7 # "mod 30" list has eight elements.
    i += 2

proc calculateAndSumTotients(sieve: var openarray[int32], n: int): int =
  result = 1
  for i in 2'i32..int32(n):
    var tot: int32
    if (i and 1) == 0:
      var m = i div 2
      var pp: int32 = 2
      while (m and 1) == 0:
        pp *= 2
        m = m div 2
      if m == 1:
        tot = pp div 2
      else:
        tot = (pp div 2) * sieve[m div 2]
    elif sieve[i div 2] == 0: # prime?
      tot = i - 1
      sieve[i div 2] = tot
    else:
      # find and extract the first prime power pp.
      # It's relatively prime with i/pp.
      var p = sieve[i div 2]
      var m = i div p
      var pp = p
      while m mod p == 0 and m != p:
        pp *= p
        m = m div p
      if m == p: # is i a prime power?
        tot = pp*(p-1)
      else:
        tot = sieve[pp div 2] * sieve[m div 2]
      sieve[i div 2] = tot
    result += tot

proc main(n: int) =
  var sieve = constructSieve(n)
  let totSum = calculateAndSumTotients(sieve, n)
  echo totSum

main(580_000_000)
Reimer Behrends
źródło
Epicki! +1. Nimrod zaczyna być coraz bardziej popularny; 3
cjfaure
Czekać. Łał Głosuję za Twoją drugą odpowiedzią. : P
cjfaure
1
Czy Nimrod to skrzyżowanie Pythona i C?
mbomb007
Nimrod został niedawno przemianowany na Nim; chociaż pożycza styl składniowy Pythona, semantyka jest inna, aw przeciwieństwie do C, jest bezpieczny dla pamięci (chyba że używasz niebezpiecznych funkcji) i ma funkcję wyrzucania elementów bezużytecznych.
Reimer Behrends
9

Java, wynik ~ 24 000 (360 000 000/15 000)

Poniższy kod Java oblicza razem funkcję totient i sito główne. Zauważ, że w zależności od twojego komputera musisz zwiększyć początkową / maksymalną wielkość sterty (na moim raczej wolnym laptopie musiałem iść do góry -Xmx3g -Xms3g).

public class Totient {

    final static int size = 360000000;
    final static int[] phi = new int[size];

    public static void main(String[] args) {
        long time = System.currentTimeMillis();
        long sum = 0;

        phi[1] = 1;
        for (int i = 2; i < size; i++) {
            if (phi[i] == 0) {
                phi[i] = i - 1;
                for (int j = 2; i * j < size; j++) {
                    if (phi[j] == 0)
                        continue;

                    int q = j;
                    int f = i - 1;
                    while (q % i == 0) {
                        f *= i;
                        q /= i;
                    }
                    phi[i * j] = f * phi[q];
                }
            }
            sum += phi[i];
        }
        System.out.println(System.currentTimeMillis() - time);
        System.out.println(sum);
    }
}
Howard
źródło
9

Nimrod: ~ 2 333 333 (42 000 000 000/18 000)

To używa zupełnie innego podejścia niż moja poprzednia odpowiedź. Zobacz komentarze, aby uzyskać szczegółowe informacje. longintModuł można znaleźć tutaj .

import longint

const max = 500_000_000

var ts_mem: array[1..max, int]

# ts(n, d) is defined as the number of pairs (a,b)
# such that 1 <= a <= b <= n and gcd(a,b) = d.
#
# The following equations hold:
#
# ts(n, d) = ts(n div d, 1)
# sum for i in 1..n of ts(n, i) = n*(n+1)/2
#
# This leads to the recurrence:
# ts(n, 1) = n*(n+1)/2 - sum for i in 2..n of ts(n, i)
#
# or, where ts(n) = ts(n, 1):
# ts(n) = n*(n+1)/2 - sum for i in 2..n of ts(n div i)
#
# Note that the large numbers that we deal with can
# overflow 64-bit integers.

proc ts(n, gcd: int): int =
  if n == 0:
    result = 0
  elif n == 1 and gcd == 1:
    result = 1
  elif gcd == 1:
    result = n*(n+1) div 2
    for i in 2..n:
      result -= ts(n, i)
  else:
    result = ts(n div gcd, 1)

# Below is the optimized version of the same algorithm.

proc ts(n: int): int =
  if n == 0:
    result = 0
  elif n == 1:
    result = 1
  else:
    if n <= max and ts_mem[n] > 0:
      return ts_mem[n]
    result = n*(n+1) div 2
    var p = n
    var k = 2
    while k < n div k:
      let pold = p
      p = n div k
      k += 1
      let t = ts(n div pold)
      result -= t * (pold-p)
    while p >= 2:
      result -= ts(n div p)
      p -= 1
    if n <= max:
      ts_mem[n] = result

proc ts(n: int128): int128 =
  if n <= 2_000_000_000:
    result = ts(n.toInt)
  else:
    result = n*(n+1) div 2
    var p = n
    var k = 2
    while k < n div k:
      let pold = p
      p = n div k
      k += 1
      let t = ts(n div pold)
      result = result - t * (pold-p)
    while p >= 2:
      result = result - ts(n div p)
      p = p - 1

echo ts(42_000_000_000.toInt128)
Reimer Behrends
źródło
Panie i panowie, to właśnie nazywam magią.
Anna Jokela,
2
Świetne podejście do obliczania sumy bezpośrednio, ale niestety nie oblicza ona funkcji sumarycznej dla tylu liczb, ile możesz, co jest wyzwaniem podanym powyżej. Twój kod faktycznie oblicza wyniki (nawet wynik funkcji sumy) tylko dla kilku tysięcy liczb (w przybliżeniu 2*sqrt(n)), co oznacza znacznie niższy wynik.
Howard
7

C #: 49 000 (980 000 000/20 000)

/codegolf//a/26800 „Kod Howarda”.
Ale zmodyfikowane wartości ph są obliczane dla nieparzystych liczb całkowitych.

using System;
using sw = System.Diagnostics.Stopwatch;
class Program
{
    static void Main()
    {
        sw sw = sw.StartNew();
        Console.Write(sumPhi(980000000) + " " + sw.Elapsed);
        sw.Stop(); Console.Read();
    }

    static long sumPhi(int n)  // sum phi[i] , 1 <= i <= n
    {
        long s = 0; int[] phi;
        if (n < 1) return 0; phi = buildPhi(n + 1);
        for (int i = 1; i <= n; i++) s += getPhi(i, phi);
        return s;
    }

    static int getPhi(int i, int[] phi)
    {
        if ((i & 1) > 0) return phi[i >> 1];
        if ((i & 3) > 0) return phi[i >> 2];
        int z = ntz(i); return phi[i >> z >> 1] << z - 1;
    }

    static int[] buildPhi(int n)  // phi[i >> 1] , i odd , i < n
    {
        int i, j, y, x, q, r, f; int[] phi;
        if (n < 2) return new int[] { 0 };
        phi = new int[n / 2]; phi[0] = 1;
        for (j = 2, i = 3; i < n; i *= 3, j *= 3) phi[i >> 1] = j;
        for (x = 4, i = 5; i <= n >> 1; i += x ^= 6)
        {
            if (phi[i >> 1] > 0) continue; phi[i >> 1] = i ^ 1;
            for (j = 3, y = 3 * i; y < n; y += i << 1, j += 2)
            {
                if (phi[j >> 1] == 0) continue; q = j; f = i ^ 1;
                while ((r = q) == i * (q /= i)) f *= i;
                phi[y >> 1] = f * phi[r >> 1];
            }
        }
        for (; i < n; i += x ^= 6)  // primes > n / 2 
            if (phi[i >> 1] == 0)
                phi[i >> 1] = i ^ 1;
        return phi;
    }

    static int ntz(int i)  // number of trailing zeros
    {
        int z = 1;
        if ((i & 0xffff) == 0) { z += 16; i >>= 16; }
        if ((i & 0x00ff) == 0) { z += 08; i >>= 08; }
        if ((i & 0x000f) == 0) { z += 04; i >>= 04; }
        if ((i & 0x0003) == 0) { z += 02; i >>= 02; }
        return z - (i & 1);
    }
}

Nowy wynik: 61 000 (1 220 000 000/20 000)
W „App.config” musiałem dodać „gcAllowVeryLargeObjects enabled = true”.

    static long sumPhi(int n)
    {
        int i1, i2, i3, i4, z; long s1, s2, s3, s4; int[] phi;
        if (n < 1) return 0; phi = buildPhi(n + 1); n -= 4; z = 2;
        i1 = 1; i2 = 2; i3 = 3; i4 = 4; s1 = s2 = s3 = s4 = 0;
        if (n > 0)
            for (; ; )
            {
                s1 += phi[i1 >> 1];
                s2 += phi[i2 >> 2];
                s3 += phi[i3 >> 1];
                s4 += phi[i4 >> z >> 1] << z - 1;
                i1 += 4; i2 += 4; i3 += 4; i4 += 4;
                n -= 4; if (n < 0) break;
                if (z == 2)
                {
                    z = 3; i4 >>= 3;
                    while ((i4 & 3) == 0) { i4 >>= 2; z += 2; }
                    z += i4 & 1 ^ 1;
                    i4 = i3 + 1;
                }
                else z = 2;
            }
        if (n > -4) s1 += phi[i1 >> 1];
        if (n > -3) s2 += phi[i2 >> 2];
        if (n > -2) s3 += phi[i3 >> 1];
        if (n > -1) s4 += phi[i4 >> z >> 1] << z - 1;
        return s1 + s2 + s3 + s4;
    }

    static int[] buildPhi(int n)
    {
        int i, j, y, x, q0, q1, f; int[] phi;
        if (n < 2) return new int[] { 0 };
        phi = new int[n / 2]; phi[0] = 1;
        for (uint u = 2, v = 3; v < n; v *= 3, u *= 3) phi[v >> 1] = (int)u;
        for (x = 4, i = 5; i <= n >> 1; i += x ^= 6)
        {
            if (phi[i >> 1] > 0) continue; phi[i >> 1] = i ^ 1;
            for (j = 3, y = 3 * i; y < n; y += i << 1, j += 2)
            {
                if (phi[j >> 1] == 0) continue; q0 = j; f = i ^ 1;
                while ((q1 = q0) == i * (q0 /= i)) f *= i;
                phi[y >> 1] = f * phi[q1 >> 1];
            }
        }
        for (; i < n; i += x ^= 6)
            if (phi[i >> 1] == 0)
                phi[i >> 1] = i ^ 1;
        return phi;
    }
P_P
źródło
4

Python 3: ~ 24000 (335 000 000/14 000)

Moja wersja to port Pythona algorytmu Howarda . Moją oryginalną funkcją była modyfikacja algorytmu wprowadzona na tym blogu .

Używam modułów Numpy i Numba, aby przyspieszyć czas wykonywania. Zauważ, że normalnie nie musisz deklarować typów zmiennych lokalnych (przy użyciu Numba), ale w tym przypadku chciałem maksymalnie skrócić czas wykonania.

Edycja: połączone funkcje buildsieve i summarum w jedną funkcję.

C ++: 9,99 s (n = 14 000); Python 3: 9,94s (n = 335 000 000)

import numba as nb
import numpy as np
import time

n = 335000000

@nb.njit("i8(i4[:])", locals=dict(
    n=nb.int32, s=nb.int64, i=nb.int32,
    j=nb.int32, q=nb.int32, f=nb.int32))

def summarum(phi):
    s = 0

    phi[1] = 1

    i = 2
    while i < n:
        if phi[i] == 0:
            phi[i] = i - 1

            j = 2

            while j * i < n:
                if phi[j] != 0:
                    q = j
                    f = i - 1

                    while q % i == 0:
                        f *= i
                        q //= i

                    phi[i * j] = f * phi[q]
                j += 1
        s += phi[i]
        i += 1
    return s

if __name__ == "__main__":
    s1 = time.time()
    a = summarum(np.zeros(n, np.int32))
    s2 = time.time()

    print(a)
    print("{}s".format(s2 - s1))
Anna Jokela
źródło
1
Powinieneś przyznać odpowiedni kredyt, gdy kopiujesz kod od innych użytkowników.
Howard
Aktualizacja z odpowiednimi kredytami!
Anna Jokela,
3

Oto moja implementacja w języku Python, która wydaje się być w stanie wykręcić ~ 60000 liczb w ciągu 10 sekund. Faktoryzuję liczby za pomocą algorytmu rovera i za pomocą testu pierwotności Millera.

from Queue import Queue
import random

def gcd ( a , b ):
    while b != 0: a, b = b, a % b
    return a

def rabin_miller(p):
    if(p<2): return False
    if(p!=2 and p%2==0): return False
    s=p-1
    while(s%2==0): s>>=1
    for _ in xrange(10):
        a=random.randrange(p-1)+1
        temp=s
        mod=pow(a,temp,p)
        while(temp!=p-1 and mod!=1 and mod!=p-1):
            mod=(mod*mod)%p
            temp=temp*2
        if(mod!=p-1 and temp%2==0): return False
    return True

def pollard_rho(n):
    if(n%2==0): return 2;
    x=random.randrange(2,1000000)
    c=random.randrange(2,1000000)
    y=x
    d=1
    while(d==1):
        x=(x*x+c)%n
        y=(y*y+c)%n
        y=(y*y+c)%n
        d=gcd(x-y,n)
        if(d==n): break;
    return d;

def primeFactorization(n):
    if n <= 0: raise ValueError("Fucked up input, n <= 0")
    elif n == 1: return []
    queue = Queue()
    factors=[]
    queue.put(n)
    while(not queue.empty()):
        l=queue.get()
        if(rabin_miller(l)):
            factors.append(l)
            continue
        d=pollard_rho(l)
        if(d==l):queue.put(l)
        else:
            queue.put(d)
            queue.put(l/d)
    return factors

def phi(n):

    if rabin_miller(n): return n-1
    phi = n
    for p in set(primeFactorization(n)):
        phi -= (phi/p)
    return phi

if __name__ == '__main__':

  n = 1
  s = 0

  while n < 60000:
    n += 1
    s += phi(n)
  print(s)
will.fiset
źródło
2

φ (2 n ) = 2 n - 1
Σ φ (2 i ) = 2 i - 1 dla i od 1 do n

Po pierwsze, coś do znalezienia razy:

import os
from time import perf_counter

SEARCH_LOWER = -1
SEARCH_HIGHER = 1

def integer_binary_search(start, lower=None, upper=None, big_jump=1):
    if lower is not None and lower == upper:
        raise StopIteration # ?

    result = yield start

    if result == SEARCH_LOWER:
        if lower is None:
            yield from integer_binary_search(
                start=start - big_jump,
                lower=None,
                upper=start - 1,
                big_jump=big_jump * 2)
        else:
            yield from integer_binary_search(
                start=(lower + start) // 2,
                lower=lower,
                upper=start - 1)
    elif result == SEARCH_HIGHER:
        if upper is None:
            yield from integer_binary_search(
                start=start + big_jump,
                lower=start + 1,
                upper=None,
                big_jump=big_jump * 2)
        else:
            yield from integer_binary_search(
                start=(start + upper) // 2,
                lower=start + 1,
                upper=upper)
    else:
        raise ValueError('Expected SEARCH_LOWER or SEARCH_HIGHER.')

search = integer_binary_search(start=1000, lower=1, upper=None, big_jump=2500)
n = search.send(None)

while True:
    print('Trying with %d iterations.' % (n,))

    os.spawnlp(
        os.P_WAIT,
        'g++', 'g++', '-Wall', '-Wextra', '-pedantic', '-O0', '-o', 'reference',
        '-DITERATIONS=%d' % (n,),
        'reference.cpp')

    start = perf_counter()
    os.spawnl(os.P_WAIT, './reference', './reference')
    end = perf_counter()
    t = end - start

    if t >= 10.1:
        n = search.send(SEARCH_LOWER)
    elif t <= 9.9:
        n = search.send(SEARCH_HIGHER)
    else:
        print('%d iterations in %f seconds!' % (n, t))
        break

Dla mnie kod referencyjny to:


Próbuję z iteracjami 14593.
64724364
14593 iteracji w 9,987747 sekund!

Haskell:

import System.Environment (getArgs)

phiSum :: Integer -> Integer
phiSum n = 2 ^ n - 1

main :: IO ()
main = getArgs >>= print . phiSum . (2^) . read . head

Robi coś z 2525224 cyframi w 0,718 sekundy. A teraz zauważam komentarz @ Howarda.

Ry-
źródło
Czy możesz opublikować wynik z całkowitą liczbą kolejnych liczb od 1, którą udało Ci się zsumować?
qwr
@qwr, to będzie 0. Jeśli chcesz kolejnych liczb, powinieneś podać to w swoim pytaniu =)
Ry-
Zrobiłem. Już go edytowałem, ponownie go edytuję.
qwr
2

Matlab: 1464 = 26355867/18000

Nie mogę przetestować twojego kodu, więc podzieliłem go na 18000, ponieważ jest to najszybszy komputer spośród testujących. Doszedłem do zapisu przy użyciu tej właściwości:

  • jeśli p jest liczbą pierwszą, phi (p) = p - 1 (dla p <10 ^ 20)

Najbardziej podoba mi się to, że jest to jedna wkładka:

sum(primes(500000000)-1)
Dennis Jaheruddin
źródło
1
Co powiesz phi(p)na wszystkich nie-najlepszych p?
Geobits
2
@Geobits Pominąłem je, ponieważ w pytaniu nie wspomniano, jakich liczb należy użyć, o ile są one naprawdę obliczone.
Dennis Jaheruddin
Ach, nie zauważyłem tego w treści. Miły.
Geobity
Nie opublikowałeś nawet wyniku ...
qwr
1
... Jak to możliwe, że nie ma Matlab & C ++ na tym samym komputerze?
Kyle Kanos
1

Python 2.7: 10.999 (165975/15090)

Pypy 2.3.1: 28.496 (430000/15090)

Kilka interesujących metod, których używam:

Silny test Pseudoprime Rabin-Millera - test pierwotności, który zapewnia skuteczny algorytm probabilistyczny do określania, czy dana liczba jest liczbą pierwszą

Formuła produktu Eulera - iloczyn liczbowy jest różny od liczb pierwszych dzielących n

Formuła produktu Eulera

Kod:

import math
import random

#perform a Modular exponentiation
def modular_pow(base, exponent, modulus):
    result=1
    while exponent>0:
        if exponent%2==1:
           result=(result * base)%modulus
        exponent=exponent>>1
        base=(base * base)%modulus
    return result

#Miller-Rabin primality test
def checkMillerRabin(n,k):
    if n==2: return True
    if n==1 or n%2==0: return False

    #find s and d, with d odd
    s=0
    d=n-1
    while(d%2==0):
        d/=2
        s+=1
    assert (2**s*d==n-1)

    #witness loop
    composite=1
    for i in xrange(k):
        a=random.randint(2,n-1)
        x=modular_pow(a,d,n)
        if x==1 or x==n-1: continue
        for j in xrange(s-1):
            composite=1
            x=modular_pow(x,2,n)
            if x==1: return False #is composite
            if x==n-1: 
                composite=0
                break
        if composite==1:
            return False        #is composite
    return True                 #is probably prime

def findPrimes(n):              #generate a list of primes, using the sieve of eratosthenes

    primes=(n+2)*[True]

    for i in range(2,int(math.sqrt(n))+1):
        if primes[i]==True:
            for j in range(i**2,n+1,i):
                primes[j]=False

    primes=[i for i in range(2,len(primes)-1) if primes[i]==True]
    return primes

def primeFactorization(n,primes):   #find the factors of a number

    factors=[]

    i=0
    while(n!=1):
        if(n%primes[i]==0):
            factors.append(primes[i])
            n/=primes[i]
        else:
            i+=1

    return factors

def phi(n,primes):
    #some useful properties
    if (checkMillerRabin(n,10)==True):      #fast prime check
        return n-1

    factors=primeFactorization(n,primes)    #prime factors
    distinctive_prime_factors=set(factors)  

    totient=n
    for f in distinctive_prime_factors:     #phi = n * sum (1 - 1/p), p is a distinctive prime factor
        totient*=(1-1.0/f);

    return totient

if __name__ == '__main__':


    s=0
    N=165975
    # N=430000
    primes=findPrimes(N)    #upper bound for the number of primes
    for i in xrange(1,N):
        s+=phi(i,primes)

    print "Sum =",s 
AlexPnt
źródło
Dzięki za algorytmy! To był jedyny, który mogłem łatwo zrozumieć i nie jest brutalną kontrolą liczenia współ-liczb pierwszych.
użytkownik