New function mpn_mulmod_bnm1

[bodrato at mail.dm.unipi.it]

Hi,

> I wrapped and integrated Niels' code for multiplication mod B^n-1.  The

> The interface is mpn_mulmod_bnm1 (rp, ap, bp, n', tp) where ap and bp
> are pointers to the input operands, and rp should point to the result.
> The common size is in n'.  Scratch space is passed in tp, its size
> should be mpn_mulmod_bnm1_itch(n').
>
> Note that this is an internal interface, and it will as such change.  We
> actually plan to change it to take at least one more size argument,
> accepting both n and n', making padding of input arguments unnecessary.

Some time ago I changed the interface, now it is
mpn_mulmod_bnm1 (rp, rn, ap, an, bp, bn, tp) and it computes
{rp,rn} = {ap,an}*{bp,bn} mod(B^rn-1); with any 0<an<=rn, 0<bn<=rn.
Scratch space is always in tp, its size: mpn_mulmod_bnm1_itch(rn).

I also added a test function, that helped me in discovering a subtle bug,
and forced me to change the implementation again. Now, the two basecase
functions (_bnm1 and _bnp1) have the old simpler interface.

Two notes:
 - the function does never return the ZERO value, unless one of the two
operand was zero. Otherwise the class [0] modulo B^rn-1 is represented by
the value B^rn-1.
 - the code can be simplified a little bit if we impose some restrictions
to an, bn. For example one could impose an >= bn > rn/2. The current test
uses (an >= rn/2 && bn >= rn/2) || (an >= rn/4 && an >= bn), and it
should be updated accordingly.

> Next we should make the many places that now uses FFT wraparound to use
> this function instead.  This means Newton iterations everywhere, the gcd
> code, and large-operands division code.

Remember that B^rn-1 < B^rn, and one can need one more limb to obtain full
precision...

Regards,
Marco

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