Two 10-lines functions

[Niels Möller]

"gianrico.fini at terra.es" <gianrico.fini at terra.es> writes:

>>There is no function mpn_mulmod_2expp1 in any version of GMP.  (But an
>>early version of mpn_mulmod_bnm1 was called something like that.)
>
> I know that the function is not there yet. What is mpn_mulmod_bnm1 supposed to 
> do?

It takes two inputs of sizes 0 < bn <= an <= rn limbs (with some
additional restrictions), and computes the product mod
2^{GMP_NUMB_BITS*rn} - 1.

It's used by various functions where wraparound arithmetic is acceptable
because the high or low part of the product is known apriori, e.g., the
Newton iterations used for inversion. For these applications, there's no
problem that the wraparound has to be on a limb boundary, since it can
be rounded up to a suitable size (and usually, its not just a limb
boundary, but rn is also rounded upp so that it is divisible by some
(small) power of two.

The reason it's useful is that it does not (except for small sizes or
odd rn) compute the full product and then add the high and low halves,
instead the product is computed mod 2^{GMP_NUMB_BITS*rn/2} - 1
(recursively) and mod 2^{GMP_NUMB_BITS*rn/2} + 1 (full mul + wrapping
for small and moderate sizes, FFT multiplication with "built in"
wraparound for large sizes), and these two values are then CRT:ed
together. Here's a benchmark comparing it to a full multiplication:

$ ./speed -s1-100000 -f 2 -r mpn_mul_n mpn_mulmod_bnm1     
overhead 0.000000010 secs, precision 10000 units of 1.67e-09 secs, CPU freq 600.00 MHz
            mpn_mul_n mpn_mulmod_bnm1
1        #0.000000066        2.9454
2        #0.000000108        2.1145
4        #0.000000223        1.6569
8        #0.000000657        1.2435
16        0.000002380       #0.9734
32        0.000007972       #0.7814
64        0.000025025       #0.6732
128       0.000074633       #0.6113
256       0.000205307       #0.6064
512       0.000542622       #0.6284
1024      0.001401250       #0.6535
2048      0.003667448       #0.6116
4096      0.009647013       #0.5399
8192      0.025509497       #0.4498
16384     0.024001000       #0.5000
32768     0.052003000       #0.5384
65536     0.108006000       #0.4815

(And to be clear, mpn_mulmod_bnm1 is an *internal* function, it can
change or even disappear between any two releases).

> If you need my functions, I'm happy to give them to you. If you don't, no 
> problems. I was challenged and I wrote them.

Too me, they seem too applicaton specific to fit in GMP, but maybe
there are some relevant use cases?

Regards,
/Niels

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