bug#12350: Composites identified as primes in factor.c (when HAVE_GMP)

[Jim Meyering]

Torbjorn Granlund wrote:
...
>   > * We have some hardwired W_TYPE_SIZE settings for the code interfacing
>   >   to longlong.h.  It is now 64 bits.  It will break on systems where
>   >   uintmax_t is not a 64-bit type.  Please see the beginning of
>   >   factor.c.
>
>   I wonder how many types of systems would be affected.
>
> It is not used currently anywhere in coreutils?  Perhaps coreutils could

uintmax_t is used throughout coreutils, but nowhere (that comes to mind)
does it fail when UINTMAX_MAX happens to be different than 2^64-1.
What I was wondering is how many systems have a uintmax_t that is
only 32 bits wide.  Now that I reread, I suppose this code would be
ok (albeit slower) with uintmax_t wider than 64.

> use autoconf for checking this?  (If we're really crazy, we could speed
> the factor program by an additional 20% by using blocked input with
> e.g. fread.)
>
> Please take a look at the generated code for factor_using_division,
> towards the end where 8 imulq should be found (on amd64).  The code uses
> mov, imul, cmp, jbe for testing the divisibility of a prime; the branch
> is taken when the prime divides the number being factored, thus highly
> non-taken.  (I suppose we could do a better job at describing the maths,
> with some references.  This particular trick is from "Division by
> invariant integers using multiplication".)

Any place you can add a reference would be most welcome.

Here's one where I'd appreciate a reference in a comment:

  #define MAGIC64 ((uint64_t) 0x0202021202030213ULL)
  #define MAGIC63 ((uint64_t) 0x0402483012450293ULL)
  #define MAGIC65 ((uint64_t) 0x218a019866014613ULL)
  #define MAGIC11 0x23b

  /* Returns the square root if the input is a square, otherwise 0. */
  static uintmax_t
  is_square (uintmax_t x)
  {
    /* Uses the tests suggested by Cohen. Excludes 99% of squares before
       computing the square root. */
    if (((MAGIC64 >> (x & 63)) & 1)
        && ((MAGIC63 >> (x % 63)) & 1)
        /* Both 0 and 64 are squares mod (65) */
        && ((MAGIC65 >> ((x % 65) & 63)) & 1)
        && ((MAGIC11 >> (x % 11) & 1)))
      {
        uintmax_t r = isqrt (x);
        if (r*r == x)
          return r;
      }
    return 0;
  }