forprime()

[Karim Belabas <Karim.Belabas@math.u-bordeaux1.fr>]

Hi pari-dev,

  adapting an implementation provided by Charles Greathouse (thanks!), I just
committed to 'master' a general iterator on arbitrary ranges of primes, which
is not limited by 'primelimit'. It uses a variety of strategies: the private
primetable up to primelimit, a sieve up to primelimit^2, and nextprime() from
then on.

  \\ f() compute the sum of primes in [a,b] in a silly way
  (18:22) gp > f(a, b = a+10^5) = my(s = 0); forprime(p = a, b, s += p); s;
  (18:22) gp > default(primelimit,2) \\ kill primetable
  (18:22) gp > f(10^20)
  time = 152 ms.
  %2 = 211500000000000105447503
  (18:22) gp > f(0, 10^9)
  time = 9,810 ms.
  %3 = 24739512092254535

It helps *a little* to have precomputed lots of primes, but not much:

  (18:23) gp > default(primelimit,10^9) \\ cache a huge primetable
  time = 2,364 ms.
  (18:23) gp > f(0,10^9)
  time = 7,464 ms.
  %4 = 24739512092254535

Total time (initialization + forprime iteration) is about the same, cached
primes help by saving about 25% in later iterations (7.4s vs. 9.8s), 
not a major saving [ in this case, the loop body is almost trivial; in non
trivial cases, the cost of prime generation will be masked by the time needed
to process the loop ]

Current libpari code still hardcodes in many places the old
NEXT_PRIME_VIA_DIFF() machinery (using exclusively the private prime table)
and is still limited by 'primelimit': if the prime table is too small, the
function raises an exception.

As I replace these by the new iterators -- forprime_next() in libpari --, I
expect that it will become mostly useless to precompute lots of primes; say
beyond 10^6, which already provides fast primes up to (10^6)^2 = 10^12.

N.B. There's no support for primes in arithmetic progressions yet
[ forprimestep() ? forstepprime() ? ], but we already have support for
infinite loops: omitting the upper bound as in
   forprime(p = 2, /* empty */, ...blah...)
runs through "all" primes in succession.

Cheers,

    K.B.
-- 
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite Bordeaux 1          Fax: (+33) (0)5 40 00 69 50
351, cours de la Liberation    http://www.math.u-bordeaux1.fr/~belabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux1.fr/  [PARI/GP]
`