New function mpn_mulmod_bnm1
bodrato at mail.dm.unipi.it
bodrato at mail.dm.unipi.it
Sat Nov 28 17:14:12 CET 2009
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/
More information about the gmp-devel
mailing list