mpz_prevprime
Seth Troisi
braintwo at gmail.com
Thu Oct 15 09:20:18 UTC 2020
On Thu, Oct 15, 2020 at 12:13 AM Niels Möller <nisse at lysator.liu.se> wrote:
> Seth Troisi <braintwo at gmail.com> writes:
>
> > I modified the patch a tiny bit. Still hoping to get this in for an
> > upcoming prime-gap search project.
>
> Looks pretty good to me.
>
I updated and attached a new patch with the tests discussed on Oct 3rd
>
> > +static int
> > +findnext (mpz_ptr p,
> > + unsigned long(*nextmod_func)(const mpz_t, unsigned long),
> > + void(*nextseq_func)(mpz_t, const mpz_t, unsigned long))
>
> I'd name the function pointer arguments to be more similar to functions
> they refer to, and without the "next" and "_func" parts. Maybe mod_ui
> and add_ui/incr_ui/update_ui?
>
> And is there a good reason you need different mod functions? I haven't
> been following along closely, so I don't know why the current code uses
> mpz_cdiv_ui rather than mpz_fdiv_ui. It would make things a bit simpler
> if we could use mpz_fdiv_ui (the standard mathematical mod operation)
> always.
>
I use cdiv in next_prime to find the offset to the next multiple of a prime
<=> (-n) % p, for previous prime this is n % p (e.g. fdiv)
for nextprime sieving the interval 1,000,000 to 1,000,500
p is the start of the interval (1,000,000) in this case
let prime = 401
(-p) % prime = mpz_cdiv_ui(p, prime) = 94
the 94th number in the sieve (p + 94 = 100094) is the first number
divisible by 401 in the interval
for previous prime we can reuse almost all of the logic but with fdiv_ui
sieving the interval 1,000,000 to 999,500
set p to the "start" of the interval (1,000,000) in this case
p % prime = mpz_fdiv_ui(p, prime) = 307
the 307th number in the sieve (p - 307 = 999693) is the first number
divisible by 401 in the interval.
Technically I can restructure the code to make both use with mpz_fdiv_ui
My quick and dirty attempt
diff -r 1e638b839cbe mpz/nextprime.c
--- a/mpz/nextprime.c Thu Oct 15 01:13:43 2020 -0700
+++ b/mpz/nextprime.c Thu Oct 15 01:30:03 2020 -0700
@@ -156,7 +156,7 @@
static int
findnext (mpz_ptr p,
- unsigned long(*nextmod_func)(const mpz_t, unsigned long),
+ char direction,
void(*nextseq_func)(mpz_t, const mpz_t, unsigned long))
{
char *composite;
@@ -234,11 +234,17 @@
unsigned long incr, prime;
int primetest;
+ /* set p to the top of the interval */
+ if (direction == 1) {
+ mpz_add_ui(p, p, 2 * (odds_in_composite_sieve-1));
+ }
+
memset (composite, 0, odds_in_composite_sieve);
prime = 3;
for (i = 0; i < prime_limit; i++)
{
- m = nextmod_func(p, prime);
+ /* m = p % prime which is the distance from top of sieve
interval to previous multiple of p */
+ m = mpz_fdiv_ui(p, prime);
/* Only care about odd multiplies of prime. */
if (m & 1)
m += prime;
@@ -250,10 +256,15 @@
prime += primegap[i];
}
+ /* reset p to old value */
+ if (direction == 1) {
+ mpz_sub_ui(p, p, 2 * (odds_in_composite_sieve-1));
+ }
+
difference = 0;
for (incr = 0; incr < odds_in_composite_sieve; difference += 2, incr
+= 1)
{
- if (composite[incr])
+ if (direction < 0 ? composite[incr] :
composite[odds_in_composite_sieve - 1 - incr])
continue;
nextseq_func(p, p, difference);
@@ -288,7 +299,7 @@
mpz_add_ui (p, n, 1);
mpz_setbit (p, 0);
- findnext(p, mpz_cdiv_ui, mpz_add_ui);
+ findnext(p, 1, mpz_add_ui);
}
int
@@ -309,7 +320,7 @@
mpz_sub_ui (p, n, 2);
mpz_setbit (p, 0);
- return findnext(p, mpz_fdiv_ui, mpz_sub_ui);
+ return findnext(p, -1, mpz_sub_ui);
}
> > +int
> > +mpz_prevprime (mpz_ptr p, mpz_srcptr n)
>
> Interface looks good to me. And if we later add a function to find the
> first prime in an arithmetic progression, I think that should fit well.
>
> Regards,
> /Niels
>
> --
> Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
> Internet email is subject to wholesale government surveillance.
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: prevprime_20201015.patch
Type: text/x-patch
Size: 16680 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-devel/attachments/20201015/0203c2c5/attachment-0001.bin>
More information about the gmp-devel
mailing list