Ach so, tut mir leid. ich hatte den Code bisher noch nicht gepostet, da das ursprüngliche Problem ja nichts mit einem konkreten Anwendungsfall zu tun hatte.
Die komische Festellung, von der ich schrieb, hat in der Tat mehr mit Fortran zu tun als mit Python. Es gehört hier vielleicht gar nicht rein.
Das Skript soll ein Array von ungeraden zahlen auf die Teilbarkeit durch alle Primzahlen bis zu deren abgerundeter Wurzel testen, um kontinuierlich eine Primzahlliste zu bauen.
Das hatte ich bisher mit NumPy und dem kartesischen Produkt gelöst. Eine Interpretationsfunktion setzt mittels np.where() anschließend eine Null für die Teilbarkeit mit Rest und eine 1 für die restlose Teilbarkeit.
Durch spaltenweises Summieren mittels np.matmul() erkennt man nun Primzahlen anhand der 0 als Summe. Das ist algorithmisch gesehen ineffizient, da die meisten Zahlen bereits durch kleine Primzahlen geteilt werden, sodass sehr viele unnötige Divisonen durchgeführt werden. Die Alternative wären Python-Loops gewesen, doch über deren Performance braucht man nicht viel zu sagen.
Ein weiteres Problem war, dass ich nicht immer vollständige Primquadratintervalle benutzen konnte, da irgendwann die Primzahllücken zu groß werden und es in Verbindung mit dem kartesischen Produkt bei etwa 2+e08 zum Memory-Error kommt. Zudem bremst das ständige Lesen und Schreiben in den RAM großer Datenmengen unnötig. Der Fortran-Loop löst dieses Probelm, da er kaum RAM braucht und schneller ist.
Code: Alles auswählen
import numpy as np
import os as os
import sieve
x_max = int(input('Insert range to compute prime list up to:'))
x_root = int(2 * np.floor((1 + np.sqrt(x_max)) / 2)) - 1
prm_def = 0
x_stop = -1
x_oddsq0 = pow(np.linspace(1, x_root+2, int((x_root+3)/2)), 2).astype(int)
z_oddsq0 = np.zeros(int((x_root+3)/2))
xrs = pow(x_root+2, 2)
prime_list = np.zeros(int(np.divide(xrs, np.log(xrs))*1.2)).astype(int)
prime_list[:] = 2 * xrs
prime_list[0] = 2
prm_sum = 1
def f_i(a, b):
#return (np.floor(a/b) - np.ceil(a/b) + 1)
return (np.where(np.floor(a/b) == a/b, 1, 0))
#return (np.where((a % b) == 0, 1, 0))
#def f_z0(a):
#return (np.where(pdiv_sum == prm_def, 1, 0))
#return ((np.absolute(a+1-prm_def) + np.absolute(a-1-prm_def) - 2*np.absolute(a-prm_def)) / 2)
for xr_tmp in range(1, x_root+1, 2):
x_start = x_stop + 4
x_stop = x_start + 4 * xr_tmp
primelist_tmp = prime_list[np.where(prime_list <= xr_tmp)[0]]
prm_sum_tmp = np.sum(primelist_tmp * 0 + 1)
# x_arr = (np.linspace(x_start, x_stop, 2*xr_tmp+1)
# .reshape(2*xr_tmp+1, 1))
# prmic = f_i(x_arr, primelist_tmp)
# mtp = np.ones(prm_sum_tmp).reshape(prm_sum_tmp, 1)
# pdiv_sum = np.matmul(prmic, mtp)[:, 0]
x_arr = np.linspace(x_start, x_stop, 2*xr_tmp+1)
pdiv_sum = sieve.sieve(2*xr_tmp+1, prm_sum_tmp, x_arr, primelist_tmp)
primelist_new = (np.where(pdiv_sum == prm_def))[0] * 2 + x_start
prm_sum_new = len(primelist_new)
prime_list[prm_sum : prm_sum + prm_sum_new] = primelist_new
print(xr_tmp)
#print('Prime list is completed from 0 to', x_stop, '.')
prm_sum += prm_sum_new
z_oddsq0[int((xr_tmp+1)/2)] = prm_sum
z_oddsq0 = z_oddsq0.astype(int)
prime_list = prime_list.copy()[np.where(prime_list != 2*xrs)[0]]
prm_sum = len(prime_list)
prime_list = np.array(prime_list.copy(), int)
prime_list[0] = 2
prm_num = len(prime_list)
func_zeta = np.array(np.linspace(0, 0, x_stop+3), int)
func_zeta[prime_list] = 1
print("There're", prm_num, "primes.")
print('Restoring prime list to disk...')
np.save(f'./prime_list_{pow(x_root+2, 2)}.npy', prime_list)
print('Restoring prime list to disk has been completed.')
#print('Restoring uncumulated zeta function to disk...')
#np.savetxt('./../arrays.d/func_zeta.txt', func_zeta)
#print('Restoring uncumulated zeta function to disk is completed.')
#exec(open('./zf_cum_exact.py').read())
pl = prime_list[np.where(prime_list <= x_max)[0]]
print("There're", sum(prime_list*0+1), "primes from", 1, "to",
pow(x_root+2, 2), ".")
Code: Alles auswählen
subroutine sieve(n_num, prm_sum_tmp, x_arr_int, primelist_tmp, pdiv_sum)
implicit none
integer(kind=8), intent(in) :: n_num
integer(kind=8), intent(in) :: prm_sum_tmp
integer(kind=8), intent(in) :: x_arr_int(:)
integer(kind=8), intent(in) :: primelist_tmp(:)
integer(kind=8), intent(out) :: pdiv_sum(n_num)
integer(kind=8) :: i, j
real(kind=8) :: x_arr_float(n_num)
! write(*,*) "n_num =", n_num
x_arr_float = real(x_arr_int)
pdiv_sum(:) = 1
do i = 1, n_num, 1
do j = 2, prm_sum_tmp, 1
if ((x_arr_int(i) / primelist_tmp(j)) == (x_arr_float(i) / primelist_tmp(j))) then
goto 10
end if
end do
pdiv_sum(i) = 0
10 continue
end do
end subroutine sieve