Trial division for double-word input in ``fmpz_factor``

[fredrik-johansson]

For double-word integers, i.e. in the 65-128 bit range, one can get 5-30% speedup in fmpz_factor by doing efficient trial division by all the 16-bit primes, using the table of primes and inverses available in ulong_extras.

For double-word by single-word divisibility testing and subsequent factor removal by divexact, one can use something like this code:

static
ulong n_mulhi(ulong a, ulong b)
{
    ulong h, l;
    umul_ppmm(h, l, a, b);
    return h;
}

static
ulong n_subc(ulong * c, ulong a, ulong b)
{
    ulong d = a - b;
    *c = d > a;
    return d;
}

/* 2^64-adic division adapted from GMP's mpn_modexact_1_odd specialized for
two limbs; dinv is the inverse mod 2^64 */
static
int n_ll_divisible_odd(ulong ahi, ulong alo, ulong d, ulong dinv)
{
    ulong c, h, l, t;

    l = alo * dinv;
    c = n_mulhi(l, d);

    l = n_subc(&c, ahi, c);
    l = l * dinv;
    c += n_mulhi(l, d);
    return (c == 0);
}

static
void n_ll_divexact_odd(ulong * rhi, ulong * rlo, ulong xhi, ulong xlo, ulong d, ulong dinv)
{
    ulong s, l, c;

    s = xlo;
    l = s * dinv;
    *rlo = l;
    c = n_mulhi(l, d);
    s = xhi;
    l = s - c;
    l = l * dinv;
    *rhi = l;
}

An alternative which I have not tested is to compute a remainder by a product of four primes at a time and checking the single-word remainder for divisibility.

Some related inefficiencies in fmpz_factor:

• Checking cofactors for primality after the trial division calls fmpz_is_prime which does redundant trial division; we should have a fmpz_is_prime_no_trial.
• Similarly, when we reach a single-limb cofactor, we pass this to n_factor which may do more redundant trial division.
• Like in n_factor, it might make sense to throw in a primality test after a few rounds of trial division rather than after all trial division.