mpz_prevprime

[Marco Bodrato]

Ciao,

Il Gio, 22 Agosto 2019 8:48 am, Marco Bodrato ha scritto:
> Il Gio, 22 Agosto 2019 7:24 am, Niels Möller ha scritto:
>> Maybe it would make sense to generalize a bit further, and write a
>> function to find the first prime (if any) in an arithmetic sequence.
>> I.e., find the smallest k >= such that A + kB is prime (if any). A

... here too we should decide what to do with negative numbers.
which prime should we return if:
 1) A = -123; B= 42; (I'd suggest A+3*B = +3)
 2) A = -123; B= 30; (I'd suggest A+4*B = -3)
 3) A = -123; B= 26; (I'd suggest A+1*B = -97)

> I know that my proposal is an incompatible change for the interface... but
> I'd suggest to extend the range for mpz_nextprime. Can we consider both
> positive and negative primes?

By the way... our documentation does not mention "positive" primes :-) and
I'm quite sure that there aren't programs using the fact that the current
implementation of nextprime returns 2 for negative numbers...

Il Gio, 22 Agosto 2019 10:32 am, Niels Möller ha scritto:
> I think that's a bit unusual. We can compare with the corresponding
> functions in pari/gp. These differ from gmp by nextprime(p) == p if p is
> prime. And they don't consider negative numbers to be primes.

Il Mer, 21 Agosto 2019 10:02 am, Seth Troisi ha scritto:
> Finally for mpz_prevprime(rop, op <= 2) it's not clear what rop should be
> set to, so I use the smallest prime (2).

Well, in pari/gp there is a function
precprime(x): largest pseudoprime <= x, 0 if x<=1.

I'm not sure I like that interface...

My proposal to extend the range of nextprime to negative numbers:

diff -r 643c931da9bd mpz/nextprime.c
--- a/mpz/nextprime.c   Thu Aug 22 15:14:09 2019 +0200
+++ b/mpz/nextprime.c   Sat Aug 24 16:52:21 2019 +0200
@@ -56,21 +56,24 @@
   mp_size_t pn;
   mp_bitcnt_t nbits;
   unsigned incr;
+  mpz_t tp;
   TMP_SDECL;

   /* First handle tiny numbers */
-  if (mpz_cmp_ui (n, 2) < 0)
+  if (mpz_cmp_ui (n, 2) < 0 && mpz_cmpabs_ui (n, 4) < 0)
     {
       mpz_set_ui (p, 2);
+      if (mpz_cmpabs_ui (n, 3) == 0)
+       mpz_neg (p, p);
       return;
     }
   mpz_add_ui (p, n, 1);
   mpz_setbit (p, 0);

-  if (mpz_cmp_ui (p, 7) <= 0)
+  if (mpz_cmpabs_ui (p, 7) <= 0)
     return;

-  pn = SIZ(p);
+  pn = ABSIZ(p);
   MPN_SIZEINBASE_2EXP(nbits, PTR(p), pn, 1);
   if (nbits / 2 >= NUMBER_OF_PRIMES)
     prime_limit = NUMBER_OF_PRIMES - 1;
@@ -88,7 +91,7 @@
       prime = 3;
       for (i = 0; i < prime_limit; i++)
        {
-         moduli[i] = mpz_tdiv_ui (p, prime);
+         moduli[i] = mpz_fdiv_ui (p, prime);
          prime += primegap[i];
        }

@@ -115,7 +118,7 @@
          difference = 0;

          /* Miller-Rabin test */
-         if (mpz_millerrabin (p, 25))
+         if (mpz_millerrabin (mpz_roinit_n (tp, PTR (p), ABSIZ (p)), 25))
            goto done;
        next:;
          incr += 2;


After those changes, the _precprime function can be implemented with:

void
mpz_precprime (mpz_ptr p, mpz_srcptr n)
{
  mpz_t tn;

  mpz_nextprime (p, mpz_roinit_n (tn, PTR (n), -SIZ (n)));
  mpz_neg (p, p);
}


Ĝis,
m

-- 
http://bodrato.it/papers/