Come posso migliorare le prestazioni del mio programma Julia per numeri eccellenti?

Question

Ho giocato con alcuni Perl programs to calculate excellent numbers. Sebbene i tempi di esecuzione delle mie soluzioni fossero accettabili, ho pensato che un linguaggio diverso, in particolare quello progettato per materiale numerico, potesse essere più veloce. Un amico ha suggerito Julia, ma la prestazione che sto vedendo è così grave che devo fare qualcosa di sbagliato. Ho guardato attraverso il Performance Tips e non vedo che cosa devo migliorare:

digits = int(ARGS[1]) 

const k = div(digits, 2) 

for a = (10^(k - 1)) : (10^(k) - 1) 
    front = a * (10^k + a) 

    root = floor(front^0.5) 

    for b = (root - 1): (root + 1) 
     back = b * (b - 1); 
     if back > front 
      break 
     end 
     if log(10,b) > k 
      continue 
     end 
     if front == back 
      @printf "%d%d\n" a b 
     end 
    end 
end 

Ho un programma di C equivalente di un ordine di grandezza più veloce al posto del fattore 2 annotati sulla Julia page (anche se la maggior parte di le domande StackOverflow circa la velocità di Julia sembrano indicare benchmark viziate da quella pagina):

E il non ottimizzato pura Perl ho scritto prende la metà del tempo:

use v5.20; 

my $digits = $ARGV[0] // 2; 
die "Number of digits must be even and non-zero! You said [$digits]\n" 
    unless($digits > 0 and $digits % 2 == 0 and int($digits) eq $digits); 

my $k = ($digits/2); 

foreach my $n (10**($k-1) .. 10**($k) - 1) { 
    my $front = $n*(10**$k + $n); 
    my $root = int(sqrt($front)); 

    foreach my $try ($root - 2 .. $root + 2) { 
     my $back = $try * ($try - 1); 
     last if length($try) > $k; 
     last if $back > $front; 
     # say "\tn: $n back: $back try: $try front: $front"; 
     if($back == $front) { 
      say "$n$try"; 
      last; 
      } 
     } 
    } 

sto usando il pre-compilato Julia per Mac OS X da allora Non ho potuto ottenere la fonte da compilare (ma non ho provato oltre a farlo esplodere la prima volta). Immagino che ne faccia parte.

Inoltre, vedo circa 0,7 secondi di tempo di avvio per qualsiasi programma Julia (vedere), il che significa che il programma C compilato in modo equivalente può essere eseguito circa 200 volte prima che Julia termini una volta. Con l'aumentare del tempo di esecuzione (maggiori valori di digits) e il tempo di avvio significa meno, il mio programma Julia è ancora molto lento.

Non ho ottenuto la parte per numeri molto grandi (numeri eccellenti di 20 cifre) che non ho capito che Julia non gestisce quelli meglio della maggior parte delle altre lingue.

---

Ecco il mio codice C, che è un po 'diverso da quando ho iniziato questo. Le mie abilità C arrugginite e ineleganti sono essenzialmente la stessa cosa del mio Perl.

#include <math.h> 
#include <stdio.h> 
#include <stdlib.h> 

int main(int argc, char *argv[]) { 
    long 
     k, digits, 
     start, end, 
     a, b, 
     front, back, 
     root 
     ; 

    digits = atoi(argv[1]); 

    k = digits/2; 

    start = (long) pow(10, k - 1); 
    end = (long) pow(10, k); 

    for(a = start; a < end; a++) { 
     front = (long) a * (pow(10,k) + a); 
     root = (long) floor(sqrt(front)); 

     for(b = root - 1; b <= root + 1; b++) { 
      back = (long) b * (b - 1); 
      if(back > front) { break; } 
      if(log10(b) > k) { continue; } 
      if(front == back) { 
       printf("%ld%ld\n", a, b); 
       } 
      } 
     } 

    return 0; 
    } 

Answer 1

È possibile provare a inserire il codice in una funzione.

function excellent(k) 
    for a = (10^(k - 1)) : (10^(k) - 1) 
    front = a * (10^k + a) 
    root = ifloor(sqrt(front)) # floor() returns a double 

    for b = (root - 1): (root + 1) 
     back = b * (b - 1); 
     if back > front 
     break 
     end 
     if log(10,b) > k 
     continue 
     end 
     if front == back 
     @printf "%d%d\n" a b 
     end 
    end 
    end 
end 

@time excellent(7) 
## 33333346666668 
## 48484848484848 
## elapsed time: 1.451842881 seconds (14680 bytes allocated) 

Per i numeri precisione arbitraria, è possibile utilizzare BigInt (ad esempio, ten = BigInt("10"), ma le prestazioni gocce ...

Answer 2

ottengo un grande aumento di velocità quando uso ifloor (che non è elencato in Mathematical Operations and Elementary Functions) invece di floor. Come liberarsi di entrambi floor a favore della isqrt mostra lo stesso aumento di velocità. non vedo dove è documentato che sia.

Ora vedo le prestazioni che ci si aspetta, anche se sul mio Mac, sembra come Julia non poteva iniziare k = 10. BigInt può aiutare, ma poi le prestazioni sono marce. Parte del motivo per cui ho guardato Julia era la mia speranza che potesse facilmente gestire numeri più grandi, quindi dovrò continuare a guardarlo.

Il resto della velocità prevista potrebbe essere nascosto nell'implementazione degli algoritmi del logaritmo, come notato da Colin nei commenti.

Mi sono divertito a controllare questa lingua, e forse ci riproverò quando sarà maturo.

Answer 3

ho benchmark codice (brian.jl) contro il codice seguente, che tenta di apportare modifiche minime al codice e segue lo stile Julian:

function excellent(digits) 
    k = div(digits, 2) 
    l = 10^(k - 1) 
    u = (10^k) - 1 

    for a in l:u 
     front = a * (10^k + a) 
     root = isqrt(front) 
     for b = (root - 1):(root + 1) 
      back = b * (b - 1) 
      back > front && break 
      log10(b) > k && continue 
      front == back && println(a,b) 
     end 
    end 
end 
excellent(int(ARGS[1])) 

separare u e l era una preferenza personale per migliorare la leggibilità . Come una linea di base, il tempo di avvio Julia sulla mia macchina è:

$ time julia -e '' 
real 0m0.248s 
user 0m0.306s 
sys 0m0.091s 

Quindi, se il calcolo si esegue per ogni esecuzione di Julia da un avviamento a freddo è dell'ordine di 0,3 secondi, quindi Julia potrebbe non essere una grande scelta per te in questa fase. Passai in 16 agli script, ed ha ottenuto:

$ time julia brian.jl 16 
1045751633986928 
1140820035650625 
3333333466666668 

real 0m15.973s 
user 0m15.691s 
sys  0m0.586s 

e

$ time julia iain.jl 16 
1045751633986928 
1140820035650625 
3333333466666668 

real 0m9.691s 
user 0m9.839s 
sys  0m0.155s 

Un limite di questo codice come scritta è che se digits>=20 supereremo la memorizzazione di un Int64. Julia, per ragioni di prestazioni, non promuove automaticamente i tipi interi da numeri interi arbitrari di precisione. Possiamo usare la nostra conoscenza del problema da affrontare questo modificando l'ultima riga a:

digits = int(ARGS[1]) 
excellent(digits >= 20 ? BigInt(digits) : digits) 

otteniamo la versione BigInt di ottima per libero, che è bello. Ignorandolo per ora, sulla mia versione del profilo ho scoperto che circa il 74% del tempo viene impiegato per calcolare il log10, seguito da ~ 19% su isqrt. Ho fatto questo profiling sostituendo l'ultima linea con

excellent(4) # Warm up to avoid effects of JIT 
@profile excellent(int(ARGS[1])) 
Profile.print() 

Ora, se volessimo a dilettarsi con le modifiche algoritmiche minori, dato quello che sappiamo ora dal profiler, siamo in grado di sostituire la linea log10 (che è solo controllando il numero di cifre in vigore) con ndigits(b) > k && continue, che ci

$ time julia iain.jl 16 
1045751633986928 
1140820035650625 
3333333466666668 

real 0m3.634s 
user 0m3.785s 
sys  0m0.153s 

dà Ciò cambia l'equilibrio di circa ~ 56% dal isqrt e ~ 28% da ndigits. Scavando ulteriormente in quel 56%, circa la metà viene spesa eseguendo lo this line che sembra un algoritmo abbastanza sensibile, quindi qualsiasi miglioramento probabilmente cambierebbe lo spirito del confronto poiché sarebbe davvero un approccio completamente diverso. Indagare il codice macchina con @code_native tende a suggerire che non c'è nient'altro di strano in corso, anche se non ho approfondito questo argomento.

Se mi permetto di impegnarsi in alcuni miglioramenti algoritmici più piccoli, posso iniziare da root+1 e solo facendo il controllo ndigits una volta, vale a dire

for a in l:u 
    front = a * (10^k + a) 
    root = isqrt(front) 
    b = root + 1 
    ndigits(b) > k && continue 
    front == b*(b-1) && println(a,b) 
    b = root 
    front == b*(b-1) && println(a,b) 
    b = root - 1 
    front == b*(b-1) && println(a,b) 
end 

che mi si mette al

real 0m2.901s 
user 0m3.050s 
sys  0m0.154s 

(Non sono convinto che siano necessari i secondi due controlli di uguaglianza, ma sto cercando di minimizzare le differenze!). Alla fine ho pensato di aumentare la velocità anticipando lo 10^k, ovvero k10 = 10^k, che sembra essere calcolato ogni nuova iterazione.Con quello, ottengo a

real 0m2.518s 
user 0m2.670s 
sys  0m0.153s 

Che è un miglioramento piuttosto 20x rispetto al codice originale.

Answer 4

Ero curioso di sapere come Perl ottenesse prestazioni così buone da questo codice, quindi ho sentito che dovevo fare un confronto. Dal momento che ci sono alcune differenze apparentemente inutili sia nel flusso di controllo che nelle operazioni tra le versioni Perl e Julia del codice nella domanda, ho trasferito ciascuna versione nell'altra lingua e ho confrontato tutte e quattro le versioni. Ho anche scritto una quinta versione di Julia usando più funzioni numeriche idiomatiche ma con la stessa struttura di controllo del flusso della versione Perl della domanda.

La prima variante è essenzialmente il codice Perl dalla domanda, ma avvolto in una funzione:

sub perl1 { 
    my $k = $_[0]; 
    foreach my $n (10**($k-1) .. 10**($k)-1) { 
     my $front = $n * (10**$k + $n); 
     my $root = int(sqrt($front)); 

     foreach my $t ($root-2 .. $root+2) { 
      my $back = $t * ($t - 1); 
      last if length($t) > $k; 
      last if $back > $front; 
      if ($back == $front) { 
       print STDERR "$n$t\n"; 
       last; 
      } 
     } 
    } 
} 

successiva, ho tradotto questo per Julia, mantenendo lo stesso flusso di controllo e con le stesse operazioni - it prende il pavimento intero della radice quadrata del front nel circuito esterno e prende la lunghezza della "in stringa" della t nel ciclo interno:

function julia1(k) 
    for n = 10^(k-1):10^k-1 
     front = n*(10^k + n) 
     root = floor(Int,sqrt(front)) 

     for t = root-2:root+2 
      back = t * (t - 1) 
      length(string(t)) > k && break 
      back > front && break 
      if back == front 
       println(STDERR,n,t) 
       break 
      end 
     end 
    end 
end 

Ecco codice Julia della domanda con qualche piccola modifica formattazione s, avvolto in una funzione:

function julia2(k) 
    for a = 10^(k-1):10^k-1 
     front = a * (10^k + a) 
     root = floor(front^0.5) 

     for b = root-1:root+1 
      back = b * (b - 1); 
      back > front && break 
      log(10,b) > k && continue 
      if front == back 
       @printf STDERR "%d%d\n" a b 
       # missing break? 
      end 
     end 
    end 
end 

Tradussi questo torna a Perl, mantenendo la stessa struttura flusso di controllo e con le stesse operazioni del codice Perl - prendendo piano di root elevato a potenza 0,5 nella esterno loop, e prendendo il logaritmo in base 10 nel ciclo interno:

sub perl2 { 
    my $k = $_[0]; 
    foreach my $a (10**($k-1) .. 10**($k)-1) { 
     my $front = $a * (10**$k + $a); 
     my $root = int($front**0.5); 

     foreach my $b ($root-1 .. $root+1) { 
      my $back = $b * ($b - 1); 
      last if $back > $front; 
      next if log($b)/log(10) > $k; 
      if ($front == $back) { 
       print STDERR "$a$b\n" 
      } 
     } 
    } 
} 

Infine, ho scritto una versione Julia che ha la stessa struttura del flusso di controllo come Perl versione della domanda, ma utilizza operazioni numeriche più idiomatiche - il isqrt e ndigits funzioni:

function julia3(k) 
    for n = 10^(k-1):10^k-1 
     front = n*(10^k + n) 
     root = isqrt(front) 

     for t = root-2:root+2 
      back = t * (t - 1) 
      ndigits(t) > k && break 
      back > front && break 
      if back == front 
       println(STDERR,n,t) 
       break 
      end 
     end 
    end 
end 

Per quanto ne so (ho usato per fare un sacco di programmazione Perl, ma è stato un po '), non ci sono versioni Perl di una di queste operazioni, quindi non c'è alcun corrispondente perl3 variante.

Ho eseguito tutte e cinque le variazioni con Perl 5.18.2 e Julia 0.3.9, rispettivamente, dieci volte ciascuna per 2, 4, 6, 8, 10, 12 e 14 cifre. Ecco i risultati di temporizzazione:

L'asse x è il numero di cifre richieste. L'asse y è il tempo mediano in secondi necessario per calcolare ciascuna funzione. L'asse y è tracciato su una scala di log (c'è qualche bug di rendering nel backend Cairo del Gadfly, quindi gli apici non appaiono molto sollevati). Possiamo vedere che, ad eccezione del numero molto piccolo di cifre (2), tutte e tre le varianti di Julia sono più veloci di entrambe le varianti Perl e che lo julia3 è sostanzialmente più veloce di qualsiasi altra cosa. Quanto più veloce? Ecco un confronto delle altre quattro varianti rispetto alla julia3 (non scala logaritmica):

L'asse x è il numero di cifre richiesto nuovamente, mentre l'asse y è quante volte più lento ogni variante era di julia3.Come potete vedere qui, non ero in grado di riprodurre le prestazioni Perl rivendicate nella domanda - il codice Perl non era 2x più veloce di Julia - era da 7 a 40 volte più lento di julia3 e almeno 2x più lento della variante Julia più lenta per qualsiasi numero non banale di cifre. Non ho provato con Perl 5.20 - forse qualcuno potrebbe seguire eseguendo questi benchmark con un nuovo Perl e vedere se questo spiega i diversi risultati? Il codice per eseguire i benchmark può essere trovato qui: excellent.pl, excellent.jl. Ho eseguito il loro in questo modo:

cat /dev/null >excellent.csv 
for d in 2 4 6 8 10 12 14; do 
    perl excellent.pl $d >>excellent.csv 
    julia excellent.jl $d >>excellent.csv 
done 

ho analizzato il excellent.csv file risultante utilizzando this Julia script.

Infine, come è stato menzionato nei commenti, utilizzando BigInt o Int128 è un'opzione per esplorare grandi numeri eccellenti in Julia. Tuttavia, ciò richiede un po 'di attenzione per la scrittura generica dell'algoritmo. Ecco una quarta variante che funziona genericamente:

function julia4(k) 
    ten = oftype(k,10) 
    for n = ten^(k-1):ten^k-1 
     front = n*(ten^k + n) 
     root = isqrt(front) 

     for t = root-2:root+2 
      back = t * (t - 1) 
      ndigits(t) > k && break 
      back > front && break 
      if back == front 
       println(STDERR,n,t) 
       break 
      end 
     end 
    end 
end 

Questa è la stessa julia3 ma funziona per tipi interi generici convertendo 10 al tipo del suo argomento. Dal momento che l'algoritmo di scale in modo esponenziale, però, ci vuole ancora molto tempo per calcolare per qualsiasi numero di cifre molto più grandi di 14:

julia> @time julia4(int128(10)) # digits = 20 
21733880705143685100 
22847252005297850625 
23037747345324014028 
23921499005444619376 
24981063345587629068 
26396551105776186476 
31698125906461101900 
33333333346666666668 
34683468346834683468 
35020266906876369525 
36160444847016852753 
36412684107047802476 
46399675808241903600 
46401324208242096401 
48179452108449381525 
elapsed time: 2260.27479767 seconds (5144 bytes allocated) 

Funziona, ma al 37 ° è una specie di lungo tempo di attesa. L'uso di un linguaggio di programmazione più veloce ti dà solo un aumento costante del fattore - 40 volte in questo caso - ma acquista solo un paio di cifre aggiuntive. Per esplorare davvero numeri più grandi, è necessario cercare algoritmi migliori.