Cvičení č. 5: Prvočísla

Vrátíme se od komplexních k celým číslům a připravíme program pro nalezení všech prvočísel menších než nebo rovných N. Uložena je budeme mít v poli, pole zapíšeme na disk do souboru. Kreslit budeme jen dvě křivky: závislost počtu prvočísel a počtu prvočíselných dvojčat na N. Dáme si však v gnuplotu záležet a řádně popíšeme obrázek i osy. Zocelíme se také v práci s velkými datovými soubory a opět budeme potřebovat Numbu pro zrychlení Pythonu.

Dynamická pole

Pole je strukturovanou proměnnou, která pod jedním jménem obsahuje více prvků zvoleného datového typu. Prvky jsou indexovány obvykle celočíselným indexem, vícerozměrná pole se indexují více indexy. Nazveme-li obyčejnou proměnnou skalárem, můžeme jednorozměrnému poli říkat vektor a dvourozměrnému poli matice. Pole bývají uložena v paměti souvisle, prvek vedle prvku, což prospívá efektivnímu přístupu do nich.

Statická pole mají dánu velikost trvale, popisem ve zdrojovém souboru, nelze ji za chodu programu změnit. Praktičtější jsou dynamická pole, jejichž velikost se může volit až za chodu programu a lze ji dále měnit. Dynamické pole se alokuje (vytváří se v paměti), realokuje (mění se jeho velikost) a dealokuje (likviduje se z paměti) pomocí alokačních příkazů nebo přiřazovacím příkazem.

V Pythonu přináší obvyklá pole a programátorskou výbavu k nim (array language) modul NumPy. Jeho pole jsou dynamická, alokujeme i realokujeme je přiřazením.

# Python a NumPy: dynamické pole
import numpy as np                    # np jako obvyklé synonymum

nmax=5
a=np.empty(nmax,dtype=int)            # alokace numpy int vektoru (pole o jedné dimenzi)
for n in range(nmax): a[n]=n+1        # inicializace cyklem
print(a)
a=np.arange(1,nmax*2+1)**2; print(a)  # realokace přiřazením
a=np.arange(1,nmax+1)**3;   print(a)  # realokace přiřazením
a=None                                # dealokace

Ve Fortranu jsou pole dynamická i statická. Pro alokaci a dealokaci dynamických polí slouží příkazy allocate a deallocate, alokovat a realokovat lze i přiřazením.

! Fortran: dynamické a statické pole
integer,allocatable :: a(:) ! popis integer dynamického vektoru (pole o jedné dimenzi)
integer :: b(5)             ! popis integer statického vektoru

nmax=5
allocate (a(nmax))                                   ! alokace alokovatelného pole
do n=1,nmax; a(n)=n; enddo;    print '(*(i0,1x))',a  ! inicializace cyklem
a=[(n,n=1,nmax*2)]**2;         print '(*(i0,1x))',a  ! realokace přiřazením (zvětšení)
a=[(n,n=1,nmax)]**3;           print '(*(i0,1x))',a  ! realokace přiřazením (zmenšení)
deallocate (a)                                       ! dealokace

do n=1,size(b); b(n)=n; enddo; print '(5(i0,1x))',b  ! inicializace statického pole v cyklu
b=[(n,n=1,size(b))]**2;        print '(5(i0,1x))',b  ! reinicializace přiřazením
end program

V obou jazycích vidíme jak skalární přiřazování do prvků pole v cyklech, tak vektorizované přiřazení pole z pravé strany do pole na levé straně.

Prvočísla v SymPy

Pythonské moduly NumPy a SciPy toho přinášejí mnoho, ale získat pole prvočísel nenabízejí. Nabízí to modul SymPy, pythonský Computer Algebra System (CAS) provádějící výpočty symbolicky, analyticky. Jeho prvním vývojářem je český matfyzák (SymPy na Wikipedii). Pokud SymPy nemáme, obvyklou cestou k němu je příkaz pip install sympy. Vypsat prvočísla získaná v SymPy a zjistit jejich počet můžeme s pomocí funkcí list a np.array vytvářejících seznam a z něho NumPy pole:

# Python a SymPy: prvočísla v seznamu a poli
import sympy
import numpy as np
nmax=1_000
primes=sympy.sieve.primerange(1,nmax)
print(primes)                      # <generator object Sieve.primerange at 0x00000211E4401000>
primes_list=list(primes)
print(primes_list)                 # [2, 3, 5, 7, 11, 13, 17, 19, ...]
print('len: ',len(primes_list))    # len:  168
primes_array=np.array(primes_list)
print(primes_array)                # [  2   3   5   7  11  13  17  19 ...]
print('np.size: ',np.size(primes_array)) # np.size:  168

Jak naznačuje volání sympy.sieve, SymPy pro získání prvočísel používá nějaké síto (sieve). A je to Eratosthenovo síto, a my si ho naprogramujeme po svém. Získáme tím nejen vhled, ale i rychlost: pro nmax rovné 100_000_000 se příslušný seznam prvočísel pomocí SymPy vytváří desítky sekund, a my budeme chtít naše obrázky dotáhnout až k miliardě. Spoiler: výpočet stihneme řádově rychleji.

Eratosthenovo síto

Eratosthenovo síto (zde viz i vizualizaci) je populární algoritmus pro nalezení všech prvočísel po předepsanou mez. Je navíc i rychlý, má lineární časovou složitost. Připravíme si pole logických (bool) hodnot indexované od 0 do nmax (hodnota prvku True nebo False bude nakonec říkat, zda index prvku je či není prvočíslem). Do všech jeho prvků vložíme s důvěrou True, vyškrtneme neprvočísla 0 a 1 (přiřadíme False do prvků s indexy 0 a 1) a uchopíme první prvočíslo n=2 (prvku s indexem 2 ponecháme True). Vyškrtáme všechny násobky tohoto prvočísla (vkládáme False do prvků s indexy n*2, n*3, …). Pak se posuneme k dalšímu prvočíslu (prvku a[n] s hodnotou True) a opakujeme od předchozí věty, než se dostaneme k nmax. Násobky jsou vyškrtány, zůstala – prvočísla.

Algoritmus můžeme urychlit třemi optimalizacemi:

  • vyškrtávání startujeme pro každé n na jeho n-násobku, nižší násobky jsme vyškrtali v předchozích chodech cyklu,

  • když startujeme vyškrtávání až na n*n, stačí zvyšovat n jen k odmocnině sqrt(nmax),

  • vyškrtávat můžeme v podmínce jen dosud nevyškrtnuté prvky, neboť zbytečný zápis do paměti stojí zbytečný čas.

Algoritmus vložíme do funkce sieve_of_eratosthenes, jejímž vstupem bude nmax a výstupem integer numpy-pole prvočísel. To získáme tak, že po dokončení prosívání zjistíme počet prvočísel (prvků a[n] s hodnotou true), alokujeme takto velké integer pole pro návratovou hodnotu funkce a prvočísla do něj z indexů v sítu přesypeme.

# Python: Eratosthenovo síto pro nalezení všech prvočísel menších než nebo rovných nmax
from numba import njit
import numpy as np

# vrací integer numpy pole prvočísel od 2 do nmax včetně
@njit
def sieve_of_eratosthenes(nmax):
  a=np.full(nmax+1,True)       # alokace dynamického pole bool prvků (zde: síto)
  a[0:2]=False                 # 0 a 1 nejsou prvočísla
  n=2                          # první prvočíslo
  while n<=np.sqrt(nmax):      # průchod od 2 do sqrt(nmax)
    nn=n*n                     # nejbližší nevyškrtnutý násobek
    while nn<=nmax:            # vyškrtávání násobků
      a[nn]=False              # - všech
    # if a[nn]: a[nn]=False    # - jen dosud nevyškrtnutých
      nn+=n                    # posun na další násobek
    while True:                # posun na další nevyškrtnutý prvek
      n+=1
      if a[n]: break           # skok z vnitřního cyklu na začátek vnějšího cyklu
  cnt=0
  for n in range(nmax+1):      # zjištění počtu nevyškrtnutých prvků
    if a[n]: cnt+=1
  # cnt=sum(1 for _ in a if _) # alternativa neschůdná pro Numbu
  result=np.zeros(cnt,dtype=np.int32)  # alokace výstupního pole prvočísel
  nprime=0
  for n in range(nmax+1):      # průchod sítem
    if a[n]:
      result[nprime]=n         # přenos prvočísel z indexů v sítu do výstupního pole
      nprime+=1
  return result

nmax=1_000
# nmax=1_000_000_000
prvocisla=sieve_of_eratosthenes(nmax)
for n in range(len(prvocisla)): print(f'{prvocisla[n]:8d}',end='')
print()
print('Pocet prvocisel mensich nez nebo rovnych '+str(nmax)+':',len(prvocisla))

Funkci otestujeme přiřazením jejího výstupu do pole a jeho výpisem v cyklu. Pro větší nmax netrápíme displej lavinou printů a omezíme se jen na souhrnný údaj o počtu nalezených prvočísel. Program na mém notebooku pro nmax rovné miliardě zjišťoval počet prvočísel (bez výpisu prvočísel) s pomocí Numby asi 13 s. Pro nmax=100_000_000 počítal s Numbou 3 s, bez Numby 46 s a se SymPy 28 s. Podmínění vyškrtávacího řádku, if a[nn]: a[nn]=False, pozitivní efekt nepřineslo.

Fortranská verze (gfortran -O3) počítá pro miliardu v nmax 12 s a sto milionů stihne za 1 s, podmíněné vyškrtávání běh zrychluje o necelých 10 %.

! Fortran: Eratosthenovo síto pro nalezení všech prvočísel menších než nebo rovných nmax
module mSieve
implicit none

contains

! vrací integer alokovatelné pole prvočísel od 2 do nmax včetně
function sieve_of_eratosthenes(nmax) result (result)
integer,allocatable :: result(:)
integer,intent(in) :: nmax
logical(1),allocatable :: a(:)  ! alokace dynamického pole logical prvků (zde: síto)
integer n,nn,cnt,nprime
allocate (a(nmax))              ! alokace paměti pro dynamické pole
a=.true.                        ! inicializace celého pole
a(1)=.false.                    ! 1 není prvočíslo
n=2                             ! první prvočíslo
do while (n<=sqrt(real(nmax,8)))  ! průchod od 2 do sqrt(nmax)
  nn=n*n                        ! nejbližší nevyškrtnutý násobek
  do while (nn<=nmax)           ! vyškrtávání násobků
  ! a(nn)=.false.               ! - všech
    if (a(nn)) a(nn)=.false.    ! - jen dosud nevyškrtnutých
    nn=nn+n                     ! posun na další násobek
  enddo
  do                            ! posun na další nevyškrtnutý prvek
    n=n+1
    if (a(n)) exit              ! skok z vnitřního cyklu na začátek vnějšího cyklu
  enddo
enddo
cnt=count(a)                    ! zjištění počtu nevyškrtnutých prvků
allocate (result(cnt))          ! alokace výstupního pole prvočísel
nprime=1
do n=1,nmax
  if (a(n)) then
    result(nprime)=n            ! přenos prvočísel z indexů v sítu do výstupního pole
    nprime=nprime+1
  endif
enddo
end function

end module

program testSieve
use mSieve
integer,allocatable :: prvocisla(:)
nmax=1000000000
prvocisla=sieve_of_eratosthenes(nmax)
! do n=1,size(prvocisla); print '(i8$)',prvocisla(n); enddo; print *
print '(a,i0,a,i0)','Pocet prvocisel mensich nez nebo rovnych ',nmax,': ',size(prvocisla)
end program

Doplňme, že když si člověk neznalý Eratosthenova síta sedne se zadáním získat řadu prvočísel k prázdnému stolu, možná, že ho jako první algoritmus napadne projít všechna (lichá) čísla a u každého zjišťovat zbytek po jeho vydělení prvočísly menšími než ono samo; nulovost některého zbytku kandidáta diskvalifikuje. Asi takhle:

! Fortran: nalezení všech prvočísel menších než nebo rovných nmax pomocí dělení
! ... náhrada cyklu do while (n<=sqrt(real(nmax,8)))
loop1: do n=2,nmax              ! průchod od 2 do nmax
  do nn=2,int(sqrt(real(n,8)))  ! pro prvočísla od 2 do sqrt(n)
    if (a(nn)) then
      if (mod(n,nn)==0) then    ! je-li n dělitelné prvočíslem
        a(n)=.false.            ! vyškrtni
        cycle loop1             ! a jdi na další n
      endif
    endif
  enddo
  ! print '(i8$)',n             ! kontrolní výpis prvočísla
enddo loop1
! ...

(Celý zdrojový soubor je s ostatními v zipu odkázaném dole.) Ano, lze i takto, ovšem časová složitost tohoto algoritmu je kvadratická – je tvořen dvěma vnořenými cykly, přičemž počet iterací v každém cyklu roste lineárně s nmax. To nedovoluje na miliardové nmax ani pomyslet. Zopakujme, že Eratosthenovo síto má časovou složitost lineární.

Textové a binární soubory

Naše ukázky jsou nachystány pro úsporný vícesloupcový výpis na obrazovku. To si upravte pro výpis do textového souboru, jak jsme už zapisovali data pro Lissajousův obrazec i Mandelbrotovu množinu. Bude také vhodné (pro gnuplot) umístit každé prvočíslo na vlastní řádek, pořadové číslo řádku pak bude odpovídat počtu prvočísel menších než nebo rovných prvočíslu na řádku. Tím zabereme až 11 bytů na prvočíslo (9 bytů pro 9ciferná prvočísla plus ve Windows 2 byty, v Linuxu 1 byte, jako symbol konce řádku). Při miliardě v nmax a zhruba 5 procentech prvočísel v ní bude soubor velký přes půl gigabytu.

Pro úspornější ukládání mohou být užitečné binární soubory, v kterých jsou data uložena ve stejné podobě jako v paměti počítače. To pak stojí jen 4 byty na integer (8 bytů na real), tedy v našem případě asi 200 megabytů. Binární soubor si sice sami snadno nepřečteme, ale stačí, umí-li jej přečíst gnuplot, a to on umí. Navíc, programy zapisují a čtou binární soubory zpravidla rychleji než textové. Zápis do textových (neboli ASCII) a binárních souborů v Pythonu a Fortranu:

# Python: zápis do textového a binárního souboru

# Zápis do textového souboru
filename='file-py.dat'                  # jméno souboru
f=open(filename,'w')                    # otevření textového souboru
print(f'{1234567890:0d}',file=f)        # zápis 10+2 bytů (ve Windows, 10+1 v Linuxu)
print(f'{1/3:18.16f}',file=f)           # zápis 18+2 bytů do souboru
print(f'{"ABCD"}',file=f)               # zápis 4+2 bytů do souboru
f.close                                 # zavření souboru (10+2+18+2+4+2=38 bytů, 35 v Linuxu)

# Zápis do binárního souboru
import numpy as np
filename='file-py.bin'                  # jméno souboru
f=open(filename,'wb')                   # otevření binárního souboru
np.array(1234567890,dtype=np.int32).tofile(f) # binární zápis 4 bytů
np.array(1/3).tofile(f)                 # binární zápis 8 bytů
np.array('ABCD',dtype=bytes).tofile(f)  # binární zápis 4 bytů
# f.write((1234567890).to_bytes(4,'little')) # 4bytový integer bez numpy
# f.write('ABCD'.encode())              # string bez numpy
f.close                                 # zavření souboru (4+8+4=16 bytů)
! Fortran: zápis do textového a binárního souboru

integer f                               ! celočíselný deskriptor
character(20) filename                  ! jméno souboru

! Zápis do textového souboru
filename='file-f90.dat'
open (newunit=f,file=filename)          ! otevření textového souboru
write (f,'(i0)') 1234567890             ! zápis 10+2 bytů (ve Windows, 10+1 v Linuxu)
write (f,'(f18.16)') 1._8/3             ! zápis 18+2 bytů do souboru
write (f,'(a)') 'ABCD'                  ! zápis 4+2 bytů do souboru
close (f)                               ! zavření souboru (10+2+18+2+4+2=38 bytů, 35 v Linuxu)

! Zápis do binárního souboru
filename='file-f90.bin'
open (newunit=f,file=filename,form='unformatted',access='stream',status='replace')  ! otevření binárního souboru
write (f) 1234567890                    ! binární zápis 4 bytů
write (f) 1._8/3                        ! binární zápis 8 bytů
write (f) 'ABCD'                        ! binární zápis 4 bytů
close (f)                               ! zavření souboru (4+8+4=16 bytů)

end program

Sto milionů 4bytových integerů zapíše na mém stroji Python do textového, resp. binárního souboru za 90 s, resp. 2 s a Fortran za 63 s, resp. 5 s.

DÚ: Prvočíselná dvojčata

Domácí úkol: upravte funkci s algoritmem Eratosthenova síta a testovací program tak, aby se kromě počtu prvočísel zjistil a vypsal i počet párů prvočíselných dvojčat (prvočísel, jejichž rozdíl je roven 2) menších než nmax. Obě funkce, p(N) a d(N), vykreslete pro kladná N do miliardy nebo (v případě problémů) alespoň do 100 milionů. Místo tohoto domácího úkolu můžete zpracovat úkol sice náročnější, ale za čokoládu – jeho zadání je o trochu níže.

Zbývá nám pravidelná týdenní porce gnuplotu. Chceme vlastně málo, jen kreslit dvě obyčejné čáry, pro které očekáváme data ve dvou jednosloupcových souborech, buď textových, primes.dat a twins.dat, nebo binárních, primes.bin a twins.bin. Tak alespoň ještě doplníme dvouřádkový nadpis s českou diakritikou, popisy os včetně horních indexů, přesuneme legendu do jiného rohu a nastavíme čtvercový formát obrázku:

# gnuplot: počet prvočísel (primes) a prvočíselných dvojčat (twin prime pairs) menších než N
xsize=800             # rozlišení obrázku ve směru x
ysize=800             # rozlišení obrázku ve směru y
xmax=1e9              # horní limit x-osy
ymax=xmax/20          # horní limit y-osy
binary=0              # typ vstupních souborů: 0 text, 1 binary
set title " p(N): počet prvočísel menších než N \n d(N): počet prvočíselných dvojčat menších než N "
set xlabel 'N'
set xtics xmax/5 format "%.0t.10^%1T"
set mxtics 4          # minor xtics
set ytics ymax/5 format "%.0t.10^%1T"
set mytics 4          # minor ytics
set key left top      # legenda vlevo nahoře
set size square       # čtvercový formát obrázku
set terminal pngcairo enhanced size xsize,ysize
set output 'eratosthenes.png'
if (!binary) {
plot [0:xmax][0:ymax] 'primes.dat' using 1:0 with lines title 'p(N)', \
                      'twins.dat' using 1:0 with lines title 'd(N)'
} else {
plot [0:xmax][0:ymax] 'primes.bin' binary format='%i' using 1:0 with lines title 'p(N)', \
                      'twins.bin' binary format='%i' using 1:0 with lines title 'd(N)'
}
set output
  • Pomáháme si proměnnými inicializovanými zkraje skriptu, abychom čelili výskytu magických konstant hlouběji ve skriptu.

  • Při popisu os příkazy set xtics, set ytics používáme formátové specifikace vyložené v nápovědě gnuplotu: help format specifiers. Pro horní indexy uvozené znakem ^ a podobná vylepšení je třeba podpora výstupního terminálu, zde: set terminal pngcairo enhanced.

  • Legendu lze přesouvat: set key left top, různě vylepšovat: help set key, i vypnout: unset key.

  • Větvíme podmíněným příkazem if podle hodnoty v proměnné binary (operátor ! neguje 0 na 1 a naopak) na kreslení z textových nebo binárních souborů.

  • Příkaz plot obsahuje explicitní nastavení dolní a horní meze v obou osách: [xmin:xmax][ymin:ymax], odkazuje na virtuální nultý sloupec čili číslo řádku: using 1:0, a mění default legendu: title string.

  • Příkaz plot s klauzulí binary čte data z binárních souborů. Jejich strukturu je vhodné popsat pomocí formátových specifikací podle návodu po help plot binary general a kliknutí na format. Stručně: píšeme format='%i' pro čtení 4bytového integer, %f (float) pro 4bytová reálná čísla a %lf (long float) pro 8bytová reálná čísla.

  • Skript s českými texty musíme pro gnuplot uložit v kódování UTF-8 (without BOM neboli byte order mark neboli signature). To je v dnešním Notepadu i VS Code default.

  • 32bitový gnuplot pro N rovno miliardě (tedy už pro 200MB-ový sloupec) selže. Použijte 64bitový gnuplot nebo N zmenšete.

_images/eratosthenes.png

Funkce p(N) má své jméno: prvočíselná funkce, své písmeno: \(\pi(n)\), i svou webovou stránku: Wikipedia. Tam je k porovnání i tabulka jejích hodnot.

Py a f90 zdrojové soubory s dnešními ukázkami: zip.

P. S. Goldbachova hypotéza

Goldbachova hypotéza praví, že pro sudé číslo n větší než 2 existuje rozklad na prvočísla p a r tak, že n = p + r. Význačný je pro dané n ten pár p a r, kde jedno z prvočísel je nejmenší možné a druhé největší možné. Alternativní domácí úkol (za čokoládu) zní takto: připravte program, který nalezne mezi sudými čísly n do 100 milionů taková, pro která je nejmenší z prvočísel ze všech jejich možných rozkladů větší než 1000 (neboli pro která neexistuje prvočíselný rozklad, který by obsahoval prvočíslo menší než 1000), a zjistěte, kolik takových n do 100 milionů je. (Jedním z nich je 60119912 = 1093 + 60118819, naopak nepatří mezi ně např. 503222 = 523 + 502699.) Odpověď byste měli získat do dalšího cvičení – ne každý algoritmus pro takto velká n za týden doběhne, při vhodně navrženém algoritmu však lze odpověď (a pak někdy čokoládu) získat během sekund.