C, 0,026119s (12 marca 2016 r.)
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define cache_size 16384
#define Phi_prec_max (47 * a)
#define bit(k) (1ULL << ((k) & 63))
#define word(k) sieve[(k) >> 6]
#define sbit(k) ((word(k >> 1) >> (k >> 1)) & 1)
#define ones(k) (~0ULL >> (64 - (k)))
#define m2(k) ((k + 1) / 2)
#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))
#define ns(t) (1000000000 * t.tv_sec + t.tv_nsec)
#define popcnt __builtin_popcountll
#define mask_build(i, p, o, m) mask |= m << i, i += o, i -= p * (i >= p)
#define Phi_prec_bytes ((m2(Phi_prec_max) + 1) * sizeof(int16_t))
#define Phi_prec(i, j) Phi_prec_pointer[(j) * (m2(Phi_prec_max) + 1) + (i)]
#define Phi_6_next ((i / 1155) * 480 + Phi_5[i % 1155] - Phi_5[(i + 6) / 13])
#define Phi_6_upd_1() t = Phi_6_next, i += 1, *(l++) = t
#define Phi_6_upd_2() t = Phi_6_next, i += 2, *(l++) = t, *(l++) = t
#define Phi_6_upd_3() t = Phi_6_next, i += 3, *(l++) = t, *(l++) = t, *(l++) = t
typedef unsigned __int128 uint128_t;
struct timespec then, now;
uint64_t a, primes[4648] = { 2, 3, 5, 7, 11, 13, 17, 19 }, *primes_fastdiv;
uint16_t *Phi_6, *Phi_prec_pointer;
inline uint64_t Phi_6_mod(uint64_t y)
{
if (y < 30030)
return Phi_6[m2(y)];
else
return (y / 30030) * 5760 + Phi_6[m2(y % 30030)];
}
inline uint64_t fastdiv(uint64_t dividend, uint64_t fast_divisor)
{
return ((uint128_t) dividend * fast_divisor) >> 64;
}
uint64_t Phi(uint64_t y, uint64_t c)
{
uint64_t *d = primes_fastdiv, i = 0, r = Phi_6_mod(y), t = y / 17;
r -= Phi_6_mod(t), t = y / 19;
while (i < c && t > Phi_prec_max) r -= Phi(t, i++), t = fastdiv(y, *(d++));
while (i < c && t) r -= Phi_prec(m2(t), i++), t = fastdiv(y, *(d++));
return r;
}
uint64_t Phi_small(uint64_t y, uint64_t c)
{
if (!c--) return y;
return Phi_small(y, c) - Phi_small(y / primes[c], c);
}
uint64_t pi_small(uint64_t y)
{
uint64_t i, r = 0;
for (i = 0; i < 8; i++) r += (primes[i] <= y);
for (i = 21; i <= y; i += 2)
r += i % 3 && i % 5 && i % 7 && i % 11 && i % 13 && i % 17 && i % 19;
return r;
}
int output(int result)
{
clock_gettime(CLOCK_REALTIME, &now);
printf("pi(x) = %9d real time:%9ld ns\n", result , ns(now) - ns(then));
return 0;
}
int main(int argc, char *argv[])
{
uint64_t b, i, j, k, limit, mask, P2, *p, start, t = 8, x = atoi(argv[1]);
uint64_t root2 = sqrt(x), root3 = pow(x, 1./3), top = x / root3 + 1;
uint64_t halftop = m2(top), *sieve, sieve_length = (halftop + 63) / 64;
uint64_t i3 = 1, i5 = 2, i7 = 3, i11 = 5, i13 = 6, i17 = 8, i19 = 9;
uint16_t Phi_3[] = { 0, 1, 1, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 7, 7, 8 };
uint16_t *l, *m, Phi_4[106], Phi_5[1156];
clock_gettime(CLOCK_REALTIME, &then);
sieve = malloc(sieve_length * sizeof(int64_t));
if (x < 529) return output(pi_small(x));
for (i = 0; i < sieve_length; i++)
{
mask = 0;
mask_build( i3, 3, 2, 0x9249249249249249ULL);
mask_build( i5, 5, 1, 0x1084210842108421ULL);
mask_build( i7, 7, 6, 0x8102040810204081ULL);
mask_build(i11, 11, 2, 0x0080100200400801ULL);
mask_build(i13, 13, 1, 0x0010008004002001ULL);
mask_build(i17, 17, 4, 0x0008000400020001ULL);
mask_build(i19, 19, 12, 0x0200004000080001ULL);
sieve[i] = ~mask;
}
limit = min(halftop, 8 * cache_size);
for (i = 21; i < root3; i += 2)
if (sbit(i))
for (primes[t++] = i, j = i * i / 2; j < limit; j += i)
word(j) &= ~bit(j);
a = t;
for (i = root3 | 1; i < root2 + 1; i += 2)
if (sbit(i)) primes[t++] = i;
b = t;
while (limit < halftop)
{
start = 2 * limit + 1, limit = min(halftop, limit + 8 * cache_size);
for (p = &primes[8]; p < &primes[a]; p++)
for (j = max(start / *p | 1, *p) * *p / 2; j < limit; j += *p)
word(j) &= ~bit(j);
}
P2 = (a - b) * (a + b - 1) / 2;
for (i = m2(root2); b --> a; P2 += t, i = limit)
{
limit = m2(x / primes[b]), j = limit & ~63;
if (i < j)
{
t += popcnt((word(i)) >> (i & 63)), i = (i | 63) + 1;
while (i < j) t += popcnt(word(i)), i += 64;
if (i < limit) t += popcnt(word(i) & ones(limit - i));
}
else if (i < limit) t += popcnt((word(i) >> (i & 63)) & ones(limit - i));
}
if (a < 7) return output(Phi_small(x, a) + a - 1 - P2);
a -= 7, Phi_6 = malloc(a * Phi_prec_bytes + 15016 * sizeof(int16_t));
Phi_prec_pointer = &Phi_6[15016];
for (i = 0; i <= 105; i++)
Phi_4[i] = (i / 15) * 8 + Phi_3[i % 15] - Phi_3[(i + 3) / 7];
for (i = 0; i <= 1155; i++)
Phi_5[i] = (i / 105) * 48 + Phi_4[i % 105] - Phi_4[(i + 5) / 11];
for (i = 1, l = Phi_6, *l++ = 0; i <= 15015; )
{
Phi_6_upd_3(); Phi_6_upd_2(); Phi_6_upd_1(); Phi_6_upd_2();
Phi_6_upd_1(); Phi_6_upd_2(); Phi_6_upd_3(); Phi_6_upd_1();
}
for (i = 0; i <= m2(Phi_prec_max); i++)
Phi_prec(i, 0) = Phi_6[i] - Phi_6[(i + 8) / 17];
for (j = 1, p = &primes[7]; j < a; j++, p++)
{
i = 1, memcpy(&Phi_prec(0, j), &Phi_prec(0, j - 1), Phi_prec_bytes);
l = &Phi_prec(*p / 2 + 1, j), m = &Phi_prec(m2(Phi_prec_max), j) - *p;
while (l <= m)
for (k = 0, t = Phi_prec(i++, j - 1); k < *p; k++) *(l++) -= t;
t = Phi_prec(i++, j - 1);
while (l <= m + *p) *(l++) -= t;
}
primes_fastdiv = malloc(a * sizeof(int64_t));
for (i = 0, p = &primes[8]; i < a; i++, p++)
{
t = 96 - __builtin_clzll(*p);
primes_fastdiv[i] = (bit(t) / *p + 1) << (64 - t);
}
return output(Phi(x, a) + a + 6 - P2);
}
Wykorzystuje to metodę Meissel-Lehmera .
Czasy
Na mojej maszynie otrzymuję około 5,7 milisekund dla połączonych przypadków testowych. Jest to procesor Intel Core i7-3770 z pamięcią DDR3 RAM o częstotliwości 1867 MHz i systemem openSUSE 13.2.
$ ./timepi '-march=native -O3' pi 1000
pi(x) = 93875448 real time: 2774958 ns
pi(x) = 66990613 real time: 2158491 ns
pi(x) = 62366021 real time: 2023441 ns
pi(x) = 34286170 real time: 1233158 ns
pi(x) = 5751639 real time: 384284 ns
pi(x) = 2465109 real time: 239783 ns
pi(x) = 1557132 real time: 196248 ns
pi(x) = 4339 real time: 60597 ns
0.00572879 s
Ponieważ wariancja stała się zbyt wysoka , używam taktowania z programu do nieoficjalnych czasów uruchamiania. To jest skrypt, który obliczył średnią połączonych czasów pracy.
#!/bin/bash
all() { for j in ${a[@]}; do ./$1 $j; done; }
gcc -Wall $1 -lm -o $2 $2.c
a=(1907000000 1337000000 1240000000 660000000 99820000 40550000 24850000 41500)
all $2
r=$(seq 1 $3)
for i in $r; do all $2; done > times
awk -v it=$3 '{ sum += $6 } END { print "\n" sum / (1e9 * it) " s" }' times
rm times
Oficjalne czasy
Ten czas jest na zrobienie przypadków punktów 1000 razy.
real 0m28.006s
user 0m15.703s
sys 0m14.319s
Jak to działa
Formuła
Niech będzie dodatnią liczbą całkowitą.x
Każda dodatnia liczba całkowita spełnia dokładnie jeden z następujących warunków.n≤x
n=1
jest podzielne przez liczbę pierwszą p w [ 1 , 3 √np.[1,x−−√3]
, gdzie p i q to (niekoniecznie różne) liczby pierwsze w ( 3 √n=pqpq.(x−−√3,x2−−√3)
jest liczbą pierwszą, a n > 3 √nn>x−−√3
Niech oznacza liczbę liczb pierwszych p, tak że p ≤ y . Istnieją π ( x ) - π ( 3 √π(y)pp≤yliczby należące do czwartej kategorii.π(x)−π(x−−√3)
Niech oznacza liczbę dodatnich liczb całkowitych m ≤ y, które są iloczynem dokładnie k liczb pierwszych, a nie pierwszych c liczb pierwszych. Istnieją P 2 ( x , π ( 3 √Pk(y,c)m≤ykcliczby należące do trzeciej kategorii.P2(x,π(x−−√3))
Na koniec, niech oznacza liczbę dodatnich liczb całkowitych k ≤ y, które są pierwszymi liczbami pierwszymi c liczb pierwszych. Istnieją x - ϕ ( x , π ( 3 √ϕ(y,c)k≤ycliczby należące do drugiej kategorii.x−ϕ(x,π(x−−√3))
Ponieważ we wszystkich kategoriach występuje liczba ,x
1+x−ϕ(x,π(x−−√3))+P2(x,π(x−−√3))+π(x)−π(x−−√3)=x
i dlatego,
π(x)=ϕ(x,π(x−−√3))+π(x−−√3)−1−P2(x,π(x−−√3))
Liczby w trzeciej kategorii mają unikalną reprezentację, jeśli wymagamy, aby a zatem p ≤ √p≤q . W ten sposób iloczyn liczb pierwszychpiqnależy do trzeciej kategorii wtedy i tylko wtedy, gdy 3 √p≤x−−√pq , więc sąπ(xx−−√3<p≤q≤xpmożliwe wartości dlaqdla stałej wartościp, iP2(x,π(3√π(xp)−π(p)+1qp, gdziepkoznaczak-tąliczbę pierwszą.P2(x,π(x−−√3))=∑π(x√3)<k≤π(x√)(π(xpk)−π(pk)+1)pkkth
Na koniec, każda dodatnia , który jest nie względnie pierwsze dla pierwszych c liczb może być wyrażona w sposób unikalny, co n = p k f , gdzie t k jest najniższym pierwsza współczynnik n . W ten sposób k ≤ c i f są pierwszymi liczbami pierwszymi k - 1 liczb pierwszych.n≤ycn=pkfpknk≤cfk−1
Prowadzi to do wzoru rekurencyjnego . W szczególności suma jest pusta, jeślic=0, więcϕ(y,0)=y.ϕ(y,c)=y−∑1≤k≤cϕ(ypk,k−1)c=0ϕ(y,0)=y
Mamy teraz wzór, który pozwala nam obliczyć , generując tylko pierwszy π ( 3 √)π(x)liczby pierwsze (miliony vs miliardy).π(x2−−√3)
Algorytm
Będziemy musieli obliczyć , gdziepmoże wynosić nawet3√π(xp)p . Chociaż są na to inne sposoby (np. Rekurencyjne stosowanie naszej formuły), najszybszym sposobem wydaje się wyliczyć wszystkie liczby pierwsze do3 √x−−√3 , co można zrobić za pomocą sita Eratostenesa.x2−−√3
Po pierwsze, identyfikujemy i przechowujemy wszystkie liczby pierwsze w i obliczπ( 3 √[1,x−−√]iπ(√π(x−−√3)w tym samym czasie. Następnie obliczamy xπ(x−−√) dla wszystkichkin(π(3√xpkki policz liczby pierwsze do każdego kolejnego ilorazu.(π(x−−√3),π(x−−√)]
Również ma postać zamkniętąπ( 3 √∑π(x√3)<k≤π(x√)(−π(pk)+1) , co pozwala nam zakończyć obliczeniaP2(x,π(3√π(x√3)−π(x√))(π(x√3)+π(x√)−12.P2(x,π(x−−√3))
Pozostawia to obliczenie , które jest najdroższą częścią algorytmu. Samo użycie rekurencyjnej formuły wymagałoby 2 wywołań funkcji c do obliczenia ϕ ( y , c ) .ϕ2cϕ(y,c)
Przede wszystkim dla wszystkich wartości c , więc ϕ ( y , c ) = y - ∑ 1 ≤ k ≤ c , p k ≤ y ϕ ( yϕ(0,c)=0c. Ta obserwacja sama w sobie wystarcza, aby obliczenia były wykonalne. Jest tak, ponieważ dowolna liczba poniżej2⋅109jest mniejsza niż iloczyn dziesięciu różnych liczb pierwszych, więc przeważająca większość znika znika.ϕ(y,c)=y−∑1≤k≤c,pk≤yϕ(ypk,k−1)2⋅109
Ponadto, grupując i pierwsze szczyty c ′ definicji ϕ , otrzymujemy alternatywny wzór ϕ ( y , c ) = ϕ ( y , c ′ ) - ∑ c ′ < k ≤ c , p k ≤ y ϕ ( yyc′ϕ. Zatem wstępne obliczenieϕ(y,c')dla ustalonegoc'i odpowiednich wartościyzapisuje większość pozostałych wywołań funkcji i związanych z nimi obliczeń.ϕ(y,c)=ϕ(y,c′)−∑c′<k≤c,pk≤yϕ(ypk,k−1)ϕ(y,c′)c′y
Jeśli , to ϕ ( m c , c ) = φ ( m c ) , ponieważ liczby całkowite w [ 1 , m c ], które nie są podzielne przez żadną z p 1 , ⋯ , p c są dokładnie tymi, które są chronione prawem autorskim do m c . Ponadto, ponieważ gcd ( z + m c , mmc=∏1≤k≤cpkϕ(mc,c)=φ(mc)[1,mc]p1,⋯,pcmc , mamy że ϕ ( y , c ) = ϕ ( ⌊ ygcd(z+mc,mc)=gcd(z,mc).ϕ(y,c)=ϕ(⌊ymc⌋mc,c)+ϕ(y
Ponieważ funkcja Eulera jest multiplikatywna, i mamy łatwy sposób na uzyskanie ϕ ( y , c ) dla wszystkich y przez precomputing wartości tylko tych y w [ 0 , m c ) .φ(mc)=∏1≤k≤cφ(pk)=∏1≤k≤c(pk−1)ϕ(y,c)yy[0,mc)
Ponadto, jeśli ustawimy , otrzymamy ϕ ( y , c ) = ϕ ( y , c - 1 ) - ϕ ( yc′=c−1, oryginalna definicja z pracy Lehmera. To daje nam prosty sposób na wstępne obliczenieϕ(y,c) wcelu zwiększenia wartościc.ϕ(y, c ) = ϕ ( y, c - 1 ) - ϕ ( ypdo, c - 1 )ϕ ( y, c )do
Oprócz wstępnego obliczania dla pewnej, niskiej wartości c , będziemy również obliczać ją wstępnie dla niskich wartości y , zmniejszając rekurencję krótko po spadku poniżej pewnego progu.ϕ ( y, c )doy
Realizacja
Poprzednia sekcja obejmuje większość części kodu. Pozostałym ważnym szczegółem jest sposób wykonywania podziałów w funkcji Phi.
Ponieważ obliczenia wymagają tylko podzielenia przez pierwszy π ( 3 √ϕliczby pierwsze,zamiast tegomożemy użyćfunkcji. Zamiast zwykłego dzieleniayprzezliczbępierwsząp, mnożymyyprzezdp≈ 2 64π( x--√3))fastdivypy zamiast tego i odzyskajyrep≈ 264p jakdpyyp . Ze względu na sposób mnożenia liczb całkowitych nax64, dzielenie przez264nie jest wymagane; wyższe 64 bitydpysą przechowywane we własnym rejestrze.repy2)642)64repy
Zauważ, że ta metoda wymaga wstępnego obliczenia , które nie jest szybsze niż obliczenie yrep bezpośrednio. Ponieważ jednak musimy dzielić te same liczby pierwsze w kółko, a dzielenie jesto wiele wolniejszeniż mnożenie, skutkuje to znacznym przyspieszeniem. Więcej szczegółów na temat tego algorytmu, a także formalny dowód, można znaleźć wDziale na liczby całkowite niezmienne przy użyciu mnożenia.yp