Trial-division interface

[Niels Möller]

Torbjorn Granlund <tg at gmplib.org> writes:

> (See mpn/generic/trialdiv.c.)

I looked briefly at that, and I had some difficulty before I understood
why it works...

Inputs:

  r: Number to be factored, 0 <= r < B
  binv: mod B inverse of a prime number p
  lim: floor((B-1)/p)

Then you do 

  q = r * binv (mod B)
  if (q <= lim)
    return p;  // p is a factor of r

It's obvious that if r *is* divisible by p, then q is the correct
quotient, and it's smaller or equal to lim. But if r is not divisible by
p, why do one always get q > lim? Before thinking about it, I would
really have expected to get something more or less random.

I think I can make this argument: The number of r in the range 0 <= r <
B which are multiples of p are lim+1, and multiplicatiion by binv (mod
B) maps these to the corresponding quotient, also in the range 0, ..., lim.

And then since multiplication by binv (mod B) is a one-to-one operation,
the other values must be mapped to elements > lim. (Is there any other
argument for this than counting elements?)

This is curious fact (to me) which I haven't seen or thought of before:
p^{-1} (mod B) has the property that multiplication mod B maps the p
multiples in the range 0 <= r < B onto the interval 0, ...,
floor((B-1)/p), in order, and the non-multiples onto the interval
floor((B-1)/p) + 1, ..., B-1, but presumably(?) in scrambled order.

And as long as p is invertible mod B, p need not be a prime and B need
not be a power of two for this property to hold.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.