0

In cases of low memory available I thought it was possible to implement a double-segment sieve to find prime numbers in a neighborhood of n. In practice, a segmented sieve has dimension dim_seg_1 for smaller primes than sqrt (n) and a sieve of size dim_seg_2 for larger numbers. This is the function to find prime numbers in the range [n_start, n_start + dim_seg_2]:

import numpy as np

def double_segmented_sieve(dim_seg_1,dim_seg_2,n_start):
    dim1=dim_seg_1//6
    dim2=dim_seg_2//6+1
    kmin=n_start//6-1
    sqrt_nmax=int(((n_start+dim_seg_2)**0.5+1)/6)+1
    if sqrt_nmax>(dim_seg_1**2)//6:
        return -1

    # vectors used to sieve prime numbers in the interval [n_start, n_start+dim_seg_2]
    sieve5m6 =np.ones((dim2+1), dtype=bool)
    sieve1m6 =np.ones((dim2+1), dtype=bool)
    
    # vectors used to sieve smaller primes of sqrt (n)
    sieve5m6_1 =np.ones((dim1+1), dtype=bool)
    sieve1m6_1 =np.ones((dim1+1), dtype=bool)
    for i in range(1,int((dim_seg_1**0.5+1)/6)+1):
        if sieve5m6_1[i]:
            sieve5m6_1[6*i*i::6*i-1]=[False]
            sieve1m6_1[6*i*i-2*i::6*i-1]=[False]
        if sieve1m6_1[i]:
            sieve5m6_1[6*i*i::6*i+1]=[False]
            sieve1m6_1[6*i*i+2*i::6*i+1]=[False]
    sieve5m6_2 =sieve5m6_1
    sieve1m6_2 =sieve1m6_1 
  
    for j in range(0,sqrt_nmax//dim1):
        for i in range(1,dim1+1):
            if sieve5m6_2[i]:
                sieve5m6[(-kmin+i+j*dim1)%(6*(i+j*dim1)-1)::6*(i+j*dim1)-1]=[False]
                sieve1m6[(-kmin-i-j*dim1)%(6*(i+j*dim1)-1)::6*(i+j*dim1)-1]=[False]
            if sieve1m6_2[i]:
                sieve5m6[(-kmin-i-j*dim1)%(6*(i+j*dim1)+1)::6*(i+j*dim1)+1]=[False]
                sieve1m6[(-kmin+i+j*dim1)%(6*(i+j*dim1)+1)::6*(i+j*dim1)+1]=[False]
        
        sieve5m6_2 =np.ones((dim1+1), dtype=bool)
        sieve1m6_2 =np.ones((dim1+1), dtype=bool)
        for i in range(1,int(((-1+(1+6*(j+1)*dim1))**0.5)/6)+1):
            if sieve5m6_1[i]:
                sieve5m6_2[(-(j+1)*dim1+i)%(6*i-1)::6*i-1]=[False]
                sieve1m6_2[(-(j+1)*dim1-i)%(6*i-1)::6*i-1]=[False]
            if sieve1m6_1[i]:
                sieve5m6_2[(-(j+1)*dim1-i)%(6*i+1)::6*i+1]=[False]
                sieve1m6_2[(-(j+1)*dim1+i)%(6*i+1)::6*i+1]=[False]

    sieve5m6[0:1:]=[False]  
    sieve1m6[0:0:]=[False]
    return np.r_[ (6 *(kmin+ np.nonzero(sieve5m6)[0]) - 1),(6 *(kmin+  np.nonzero(sieve1m6)[0]) + 1)]

P=np.sort(double_segmented_sieve(2**26,2**14,2**60))

Is it possible to speed up the function somehow?

user140242
  • 159
  • 6
  • 1
    What have you tried so far? Have you considered writing it in a compiled language, like Rust? Have you considered vectorizing your function using NumPy or JIT compiling it with a library like Numba? – dsillman2000 Jun 25 '22 at 18:04
  • @dsillman2000 I don't know Rust. I wrote the algorithm using numpy because it's easier to work with vectors. I'm interested in knowing if as an idea it can work. – user140242 Jun 25 '22 at 18:09
  • so is your question if a double-segmented prime sieve will work, or if it's possible to speed up the implementation you have? Because all of my above suggestions answer the latter. – dsillman2000 Jun 25 '22 at 18:47
  • 2
    It seems to me this question is more suited to be asked in the [Code Review Forum](https://codereview.stackexchange.com/). Code Review is a question and answer site for peer programmer code reviews. Please read the relevant guidance related to how to properly ask questions on this site before posting your question. – itprorh66 Jun 25 '22 at 19:20
  • Correctness first, efficiency later. `if sqrt_nmax>(dim_seg_1**2)//6:` doesn't look right: why are you comparing square root with a square? To sieve primes up to 2^60, you most certainly must use all primes up to 2^30. – Will Ness Jul 01 '22 at 15:01
  • Actually, 2^30 (approximately 10^9) is nothing major. My simple [odds-only non-segmented sieve in C](https://ideone.com/fapob) seems to be only taking about 4 seconds to get there - and that was in 2018! So today it's bound to be even faster. And it's cheap to also sieve an additional smaller segment at the same time. I've posted a working C code, and also some pseudocode that do this, [in a few of my answers](https://stackoverflow.com/search?q=user%3A849891+%22offset+sieve%22), for reference. – Will Ness Jul 01 '22 at 15:52
  • And of course this basic sieve to 2^30 can itself be made segmented, with the segment size of, say, 2^16, as fits in your cache. As you sieve each consecutive segment on your way to the 2^30 limit, you'd use the freshly found primes in it to also sieve the offset segment just under the 2^60 limit. You seem to be using base-6 (i.e. 2*3) wheel factorizarion. This can be extended as well, usually to 2*3*5*7. A Python code e.g. [is here](https://stackoverflow.com/a/30563958/849891). – Will Ness Jul 01 '22 at 16:06
  • @WillNess In the link at the end you find the C++ version which is much faster. In fact I have noticed that the size of the segments with this version does not vary much in speed. The sieve uses a version with two vectors considering only the numbers +-1 mod 6 for this the comparison is made with sqrt (n) / 6. https://codereview.stackexchange.com/questions/69791/sieve-of-eratosthenes-segmented-to-increase-speed-and-range/277754#277754 – user140242 Jul 01 '22 at 16:54
  • What exactly do you mean by "much faster", specifically? – Will Ness Jul 01 '22 at 18:07

0 Answers0