I've attached an updated diff, thanks for the comments Marco.
There doesn't seem to be a way to set s->reps without setting a high
precision so
I reused SPEED_BLOCK_SIZE from speed_gcd and ran the inner loop that many
times
Would folks rather I time
1. s->reps next_primes chained (next_prime(p, p); next_prime(p, p);)
2. s->reps of next_prime on random s->size bit numbers?
1. will be slightly slower (as it will have to search the full gap each
time) and more likely
for intermediate results to overflow s->size bits (which is seems fine)
than 2.
Currently it takes about 20 seconds to test -s 16-128, and 15 minutes to
test -s 16-1028
with -t 1.04 the resulting graph is much smoother than the old graph.
On Sat, Aug 24, 2019 at 8:33 AM Marco Bodrato <bodr...@mail.dm.unipi.it>
wrote:
> 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/
>
>
diff --git a/tune/common.c b/tune/common.c
index 9fefb0d..545abfa 100644
--- a/tune/common.c
+++ b/tune/common.c
@@ -1722,6 +1722,11 @@ speed_mpn_gcd_1N (struct speed_params *s)
SPEED_ROUTINE_MPN_GCD_1N (mpn_gcd_1);
}
+double
+speed_mpz_nextprime (struct speed_params *s)
+{
+ SPEED_ROUTINE_MPZ_NEXTPRIME (mpz_nextprime);
+}
double
speed_mpz_jacobi (struct speed_params *s)
diff --git a/tune/speed.c b/tune/speed.c
index 82e9c1a..588b7eb 100644
--- a/tune/speed.c
+++ b/tune/speed.c
@@ -307,6 +307,9 @@ const struct routine_t {
#if 0
{ "mpn_gcdext_lehmer", speed_mpn_gcdext_lehmer },
#endif
+
+ { "mpz_nextprime", speed_mpz_nextprime },
+
{ "mpz_jacobi", speed_mpz_jacobi },
{ "mpn_jacobi_base", speed_mpn_jacobi_base },
{ "mpn_jacobi_base_1", speed_mpn_jacobi_base_1 },
diff --git a/tune/speed.h b/tune/speed.h
index 6ea000b..0272f4e 100644
--- a/tune/speed.h
+++ b/tune/speed.h
@@ -396,6 +396,7 @@ double speed_mpz_fib_ui (struct speed_params *);
double speed_mpz_fib2_ui (struct speed_params *);
double speed_mpz_init_clear (struct speed_params *);
double speed_mpz_init_realloc_clear (struct speed_params *);
+double speed_mpz_nextprime (struct speed_params *);
double speed_mpz_jacobi (struct speed_params *);
double speed_mpz_lucnum_ui (struct speed_params *);
double speed_mpz_lucnum2_ui (struct speed_params *);
@@ -3099,6 +3100,46 @@ int speed_routine_count_zeros_setup (struct speed_params *, mp_ptr, int, int);
return t; \
}
+/* Calculate nextprime(n) for random n of s->size bits (not limbs). */
+#define SPEED_ROUTINE_MPZ_NEXTPRIME(function) \
+ { \
+ unsigned i, j; \
+ mpz_t wp, n; \
+ double t; \
+ \
+ SPEED_RESTRICT_COND (s->size >= 10); \
+ \
+ mpz_init (wp); \
+ mpz_init_set_n (n, s->xp, s->size); \
+ /* limit to s->size bits, as this function is very slow */ \
+ mpz_tdiv_r_2exp (n, n, s->size); \
+ /* set high bits so operand and result are genaral s->size bits */ \
+ mpz_setbit (n, s->size - 1); \
+ mpz_clrbit (n, s->size - 2); \
+ \
+ speed_starttime (); \
+ i = s->reps; \
+ do \
+ { \
+ /* nextprime timing is variable, so average over many calls */ \
+ j = SPEED_BLOCK_SIZE - 1; \
+ /* starts on random, after measures prime to next prime */ \
+ function (wp, n); \
+ do \
+ { \
+ function (wp, wp); \
+ } \
+ while (--j != 0); \
+ } \
+ while (--i != 0); \
+ t = speed_endtime (); \
+ \
+ mpz_clear (wp); \
+ mpz_clear (n); \
+ s->time_divisor = SPEED_BLOCK_SIZE; \
+ return t; \
+ }
+
#define SPEED_ROUTINE_MPZ_JACOBI(function) \
{ \
mpz_t a, b; \
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
