More accurate calculation of sieve array size

Marco Bodrato bodrato at mail.dm.unipi.it
Thu Aug 11 16:32:36 CEST 2022


Ciao Adrien,

Il 2022-08-10 15:37 Adrien Prost-Boucle ha scritto:
> Looking into oddfac and primorial computation code,
> I think there is a suboptimality in the computation of the prime sieve 
> size.
> Could anyone confirm?

Maybe there are some too-conservative assumptions, but we must be 
careful with
the code.

> --- mpz/oddfac_1.old	2022-08-10 15:14:04.511298108 +0200
> +++ mpz/oddfac_1.c	2022-08-10 15:31:34.721354647 +0200
> @@ -381,8 +381,8 @@
>  	  TMP_MARK;
> 
>  	  flag--;
> -	  size = n / GMP_NUMB_BITS + 4;
> -	  ASSERT (primesieve_size (n - 1) <= size - (size / 2 + 1));
> +	  size = (n / 3 / GMP_NUMB_BITS + 1) * 2;
> +	  ASSERT (primesieve_size (n - 1) <= size / 2);
>  	  /* 2-multiswing(n) < 2^(n-1)*sqrt(n/pi) < 2^(n+GMP_NUMB_BITS);
>  	     one more can be overwritten by mul, another for the sieve */
>  	  MPZ_TMP_INIT (mswing, size);

Here, mswing is not just used for the sieve, it will contain the result
of the computation of the "2-multiswing.factorial of n". A portion of
the same allocated space is used for the sieve (the ASSERT checks, is
the size needed for the sieve less than half the size we allocate? Yes,
ok we can go on...).

The comment says that /* 2-multiswing(n) < 2^(n+GMP_NUMB_BITS) */, then
it makes sense to allocate n / GMP_NUMB_BITS + 4 limbs.

Surely this is an over-estimate, and maybe we can estimate it better.
But to reduce the memory usage, we have to estimate the maximal size
of the value that will be stored in mswing, not the sieve only.

> --- mpz/primorial_ui.c.old	2022-08-10 14:44:04.437867901 +0200
> +++ mpz/primorial_ui.c	2022-08-10 14:55:41.721238758 +0200
> @@ -82,8 +82,7 @@
>        mp_limb_t prod;
>        TMP_DECL;
> 
> -      size = n / GMP_NUMB_BITS;
> -      size = size + (size >> 1) + 1;
> +      size = n / 3 / GMP_NUMB_BITS + 1;
>        ASSERT (size >= primesieve_size (n));
>        sieve = MPZ_NEWALLOC (x, size);
>        size = (gmp_primesieve (sieve, n) + 1) / log_n_max (n) + 1;

The same here. sieve is stored at the same address as x, but we need
to estimate the size of the value that x will store, not the sieve.

x will store the result, primorial (n). The question is: is the
allocated size reasonable to store the result x=primorial(n)?
The other question, "Can it also contain the sieve?", is secondary.

Ĝis,
m

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


More information about the gmp-devel mailing list