Computing at compile time

49

nagfor revised times:
N=100,000 14 seconds
N=200,000 58 seconds
N=1,000,000 28 minutes 13 seconds

Primes are checked against FriCAS.

2 Likes
50

It is a quality-of-implementation issue and each vendor will take a view on how important it is for them to devote scarce resources on this or that and how difficult it is.

1 Like
51

I think the current restrictions make it impossible to initialise a series like Fibonacci numbers this way. The reason is that you need to have the results of a previous computation for the next one. I may be mistaken, but my experiments to achieve that have failed - the compilers I used rightfully determined that I was doing something illegal.

52

Here is a faster version of Arjenā€™s program for compile-time calculation of primes. This version is also faster (to compile and run) than the improved version that I posted two days ago, and is based on the following ideas (I describe the details for N = 1 million, sqrt(N) = 1000, int(sqrt(sqrt(N))) = 31:

  • Start with a seed set of primes, such as 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, all less than N^(1/4)

  • Using these seed primes as divisors, use Arjenā€™s PACK based algorithm to obtain primes between 31 and sqrt(N) = 1000.

  • Append these to the seed set to obtain a new and augmented set of primes < sqrt(N) = 1000.

  • Repeat the previous two steps to obtain a set of primes < N = 1 million.

Using NAG 7.0 on Windows (@themos, has the 7.1 version been released yet for Windows?), I observed a compilation time of 2.5 s, and a run time of 0.03 s (only printing is done at run time).

You may view the original algorithm as a special case of this new version with 2 as the sole member of the seed set, and a single, massive use of PACK to find all the remaining prime numbers < N.

Here is the revised program, but note that compilers may require special options or excessive CPU time and memory to compile and link this version.

program primesieve
    implicit none

    integer, parameter :: rtrtN = 31, rtN = 1000, N = 1000000
    integer            :: i, np
    integer, parameter :: candidates1(*) = [(i, i=rtrtN,rtN-1,2)]
    integer, parameter :: primes0(*)  = [2,3,5,7,11,13,17,19,23,29,31] ! primes <= rtrtN, used to seed calculation
    integer, parameter :: primes1(*) = &
       pack( candidates1, [(all(mod(candidates1(i),primes0) /= 0), &
          i = 1,size(candidates1))] )  ! newly identified primes < rtN
    integer, parameter :: candidates2(*) = [(i,i=rtN+1,N-1,2)]
    integer, parameter :: primes01(*) = [primes0, primes1]        ! primes <= rtN
    integer, parameter :: primes2(*) = pack(candidates2,[(all(mod(candidates2(i),primes01) /=0), &
                   i = 1,size(candidates2))])

    np = size(primes2)
    write(*,'(i5,A10,i4)') size(primes01),' primes < ',rtN
    write(*,'(i5,A10,i7)') size(primes01)+size(primes2),' primes < ',N
    print *
    print *,' Last 10 primes < ',N,' :'
    write(*,'(10I8)') primes2(np-9:np)

end program primesieve

If anyone wishes to modify this program, for instance to change N, they need to be careful not to introduce an even integer into the candidate arrays. Similarly, if a different seed set is used, make sure that no non-primes are present in that set!

I learned to pay attention to these points when attempting to change the current program to work for N = 10000 in order to test Gfortran and Ifort. The compilation times were 61 s for Gfortran 12, 1 s for Ifort and Ifx, for N = 10000.

6 Likes
53

I decided to check that hint and made a simple Python script reconstructing the original @mecej4 code, using lists. Indeed, the timing is pretty similar to what @themos found for nagfor. But then, I switched to Python sets and the results were amazingly fast.

   N       lists      sets
 100000      22.5     0.03
 200000    1:33.1     0.05
1000000   42:21.1     0.32
10000000    ???       4.4
Python lists code 
import sys
if len(sys.argv) > 1:
    N = int(sys.argv[1])
else:
    N=1000
rtN=int(N**0.5)
cand=[2]+list(range(3,N,2))
mult=[i*j for j in range(3,rtN+1,2) for i in range (j,N//j+1,2)]
primes=[c for c in cand if c not in mult]
np=len(primes)
print(np,primes[0],primes[-1])
Python sets code 
import sys
if len(sys.argv) > 1:
    N = int(sys.argv[1])
else:
    N=1000
rtN=int(N**0.5)
cand=set(range(3,N,2))
cand.add(2)
mult=set([i*j for j in range(3,rtN+1,2) for i in range (j,N//j+1,2)])
primes=sorted(cand-mult)
np=len(primes)
print(np,primes[0],primes[-1])

Shouldnā€™t we vote for intrinsic sets in Fortran? :slight_smile:

3 Likes
54

On my Linux machine, that takes 2 seconds whereas FriCASā€™s

[i for i in 2ā€¦1000000| prime? i];

takes 10 seconds.

2 Likes
55

If you want a fast pure Python implementation:

In [8]: from sympy import sieve

In [9]: %time b = list(sieve.primerange(1000000))
CPU times: user 178 ms, sys: 3.2 ms, total: 182 ms
Wall time: 198 ms

But if you want code similar to the FriCAS snippet, you can do:

In [18]: from sympy.ntheory.primetest import isprime

In [19]: %time c = [i for i in range(2,1000000) if isprime(i)]
CPU times: user 2.18 s, sys: 3.22 ms, total: 2.18 s
Wall time: 2.18 s
1 Like
56

I wanted to use not pure but raw Python, w/o importing any modules which may (and probably do) contains precompiled binary code.

57

SymPy does not seem to use precompiled binary code for this particular algorithm: sympy/generate.py at 306c0f2c02a8ec464b988ee898ae32f76fcaf3d1 Ā· sympy/sympy Ā· GitHub

58

Many years ago I asked for a mechanism to execute something that could produce initial values for variables, or named constants, as a step in compiling a program. The example I used at the time was computing the coefficients for Gauss quadrature, based upon the desired order. Otherwise, one would need to look them up and type them in, or write a program that computes them, which program outputs them as Fortran text that can be included.

2 Likes