Função totient super rápida

22

O objetivo é simples: calcule a função totiente para o maior número possível de números em 10 segundos e some os números.

Você deve imprimir o resultado no final e realmente calculá-lo. Nenhuma função automatizada de totiente é permitida, mas as bibliotecas de bignum são. Você precisa começar com 1 e contar todos os números consecutivos. Você não tem permissão para pular números.

Sua pontuação é quantos números seu programa pode calcular em sua máquina / quantos meu programa pode calcular em sua máquina . Meu código é um programa simples em C ++ (otimizações desativadas), espero que você possa executá-lo.

Propriedades importantes que você pode usar!

  • E se gcd(m,n) = 1, phi(mn) = phi(m) * phi(n)
  • se pfor primo, phi(p) = p - 1(para p < 10^20)
  • se né par,phi(2n) = 2 phi(n)
  • outros listados no primeiro link

Meu código

#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
fonte
2
Talvez você queira adicionar que os números de entrada sejam inteiros consecutivos. Caso contrário, posso ficar tentado a calcular a função totiente apenas para potências de 2.
727 Howard
Posso fazer 1, 3, 5, 2, 4ou algo assim?
Nunez #

Respostas:

14

Nimrod: ~ 38.667 (580.000.000 / 15.000)

Esta resposta usa uma abordagem bastante simples. O código emprega uma peneira simples de número primo que armazena o número primo da menor potência primária em cada slot para números compostos (zero para números primos), depois usa a programação dinâmica para construir a função totiente no mesmo intervalo e depois soma os resultados. O programa gasta praticamente todo o tempo construindo a peneira e calcula a função totiente em uma fração do tempo. Parece que tudo se resume à construção de uma peneira eficiente (com a leve torção de que também é preciso extrair um fator primo para números compostos do resultado e manter o uso da memória em um nível razoável).

Atualização: Desempenho aprimorado, reduzindo o consumo de memória e melhorando o comportamento do cache. É possível reduzir de 5% a 10% mais o desempenho, mas o aumento na complexidade do código não vale a pena. Por fim, esse algoritmo exercita principalmente o gargalo de von Neumann de uma CPU, e há muito poucos ajustes algorítmicos que podem contornar isso.

Também atualizou o divisor de desempenho, pois o código C ++ não era para ser compilado com todas as otimizações ativadas e mais ninguém o fazia. :)

Atualização 2: operação de peneira otimizada para melhor acesso à memória. Agora, manipule primos pequenos em massa via memcpy () (~ 5% de aceleração) e pule múltiplos de 2, 3 e 5 ao peneirar primos maiores (~ 10% de aceleração).

Código C ++: 9,9 segundos (com g ++ 4,9)

Código Nimrod: 9,9 segundos (com -d: release, backend do 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
fonte
Épico! +1. Partida de Nimrod ficar mais popular; 3
cjfaure
Esperar. Woah. Estou votando sua outra resposta. : P
cjfaure
1
Nimrod é um cruzamento entre Python e C?
mbomb007
Nimrod foi renomeado recentemente para Nim; Embora empreste o estilo sintático do Python, a semântica é diferente e, ao contrário de C, é segura para a memória (a menos que você use recursos inseguros) e possui coleta de lixo.
Reimer Behrends
9

Java, pontuação ~ 24.000 (360.000.000 / 15.000)

O código java abaixo faz o cálculo da função totiente e da peneira primária juntos. Observe que, dependendo da sua máquina, você precisa aumentar o tamanho inicial / máximo da pilha (no meu laptop bastante lento, eu tive que subir -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
fonte
9

Nimrod: ~ 2.333.333 (42.000.000.000 / 18.000)

Isso usa uma abordagem totalmente diferente da minha resposta anterior. Veja os comentários para obter detalhes. O longintmódulo pode ser encontrado aqui .

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
fonte
Senhoras e senhores, é isso que eu chamo de magia.
Anna Jokela
2
Ótima abordagem para calcular a soma diretamente, mas infelizmente não calcula a função totiente para o maior número possível de números, que é o desafio dado acima. Seu código realmente calcula os resultados (nem mesmo o resultado da função totiente) para apenas vários milhares de números (aprox. 2*sqrt(n)), O que resulta em uma pontuação muito menor.
1011 Howard
7

C #: 49.000 (980.000.000 / 20.000)

/codegolf//a/26800 "Código de Howard".
Mas modificados, os valores phi são calculados para números inteiros ímpares.

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);
    }
}

Nova pontuação: 61.000 (1.220.000.000 / 20.000)
No "App.config", tive que adicionar "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
fonte
4

Python 3: ~ 24.000 (335.000.000 / 14.000)

Minha versão é uma porta Python do algoritmo de Howard . Minha função original era uma modificação de um algoritmo introduzido neste post do blog .

Estou usando os módulos Numpy e Numba para acelerar o tempo de execução. Observe que normalmente você não precisa declarar os tipos de variáveis ​​locais (ao usar o Numba), mas nesse caso eu queria reduzir o tempo de execução o máximo possível.

Editar: funções construtivas e resumidas combinadas em uma única função.

C ++: 9,99s (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
fonte
1
Você deve dar o devido crédito ao copiar o código de outros usuários.
Howard
Atualizado com os devidos créditos!
Anna Jokela
3

Aqui está minha implementação em Python que parece capaz de gerar ~ 60000 números em 10 segundos. Estou fatorando números usando o algoritmo pollard rho e usando o teste de primalidade de Rabin Miller.

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
fonte
2

φ (2 n ) = 2 n - 1
φ 2 (2 i ) = 2 i - 1 para i de 1 a n

Primeiro, algo para encontrar horários:

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

Para o código de referência, para mim, é:


Tentando com 14593 iterações.
64724364
14593 iterações em 9.987747 segundos!

Agora, Haskell:

import System.Environment (getArgs)

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

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

Faz algo com 2525224 dígitos em 0,718 segundos. E agora estou apenas percebendo o comentário de @ Howard.

Ry-
fonte
Você pode postar uma pontuação com o total de números consecutivos a partir de 1 que você conseguiu somar?
Q8
@qwr, que seria 0. Se você quiser números consecutivos, você deve especificá-lo na sua pergunta =)
Ry-
Eu fiz. Eu já editei, vou editá-lo novamente.
QWR
2

Matlab: 1464 = 26355867/18000

Como não posso testar seu código, dividi-o em 18000, pois representa o computador mais rápido daqueles que testaram. Eu vim para a pontuação usando esta propriedade:

  • se p é primo, phi (p) = p - 1 (para p <10 ^ 20)

Eu gosto principalmente que é um forro:

sum(primes(500000000)-1)
Dennis Jaheruddin
fonte
1
E phi(p)para todos os que não são primos p?
Geobits
2
@ Geobits Eu pulei aqueles como a pergunta não menciona quais números você deve usar, desde que sejam realmente calculados.
Dennis Jaheruddin
Ah, não percebi isso no texto. Agradável.
Geobits
Você nem mesmo postar uma pontuação ...
QWR
1
... Como é possível não ter o Matlab & C ++ no mesmo computador?
Kyle Kanos
1

Python 2.7: 10.999 (165975/15090)

Pypy 2.3.1: 28.496 (430000/15090)

Alguns métodos interessantes que eu uso:

Teste de Pseudoprime Forte de Rabin-Miller - Um teste de primalidade que fornece um algoritmo probabilístico eficiente para determinar se um determinado número é primo

Fórmula do produto de Euler - O produto ultrapassa os números primos distintos que dividem n

Fórmula do produto de Euler

Código:

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
fonte
Obrigado pelos algoritmos! Foi o único que eu pude entender facilmente e não é a verificação de força bruta dos primos de contagem.
User