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

[Torbjorn Granlund]

Jim Meyering <address@hidden> writes:

  The existing code can factor arbitrarily large numbers quickly, as long
  as they have no large prime factors.  We should retain that capability.
  
OK, so we'll put the GMP demos program into this one.

This opens another technical concern:

We have moved towards proving primality, since for 128 bit numbers it
can be done quickly.  But if we allow arbitrary large numbers, it is
expensive.

We might want an option for this choosing probabilistic testing, perhaps
--prp (common abbreviation for PRobabilistic Prime).  By default, we
should prove primality, I think.

My current devel version if GMP's demos/factorize has Lucas code.

  > * We think a --verbose option would be nice, although we don't have
  >   one in the present version.  It would output information on
  >   algorithm switches and bounds reached.  Opinions?
  
  I think it would be worthwhile, especially to give an idea of what progress
  is being made when factoring very large numbers, but hardly something
  that need be done now.
  
  E.g., currently this doesn't print much:
  
      $ M8=$(echo 2^31-1|bc) M9=$(echo 2^61-1|bc) M10=$(echo 2^89-1|bc)
      $ factor --verbose $(echo "$M8 * $M9 * $M10" | bc)
      [using arbitrary-precision arithmetic][trial division (32761)] [is number 
prime?] [pollard-rho (1)]
  
  Ideally it'd print something every second or two.
  
I'll let Niels worry about this, since he was the one to ask for it.

  > Portability caveats:
  >
  > * We rely on POSIX.1 getchar_unlocked for a performance advantage.
  >
  > * 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
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".)

-- 
Torbjörn