# GMP Development Projects

```Copyright 2000-2006, 2008-2011
Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or (at

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.
```

If you want to work on any of the projects below, please let gmp-devel@gmplib.org know. If you want to help with a project that already somebody else is working on, you will get in touch through gmp-devel@gmplib.org. (There are no email addresses of volunteers below, due to spamming problems.)

• Faster multiplication
1. Work on the algorithm selection code for unbalanced multiplication.
2. Implement an FFT variant computing the coefficients mod m different limb size primes of the form l*2^k+1. i.e., compute m separate FFTs. The wanted coefficients will at the end be found by lifting with CRT (Chinese Remainder Theorem). If we let m = 3, i.e., use 3 primes, we can split the operands into coefficients at limb boundaries, and if our machine uses b-bit limbs, we can multiply numbers with close to 2^b limbs without coefficient overflow. For smaller multiplication, we might perhaps let m = 1, and instead of splitting our operands at limb boundaries, split them in much smaller pieces. We might also use 4 or more primes, and split operands into bigger than b-bit chunks. By using more primes, the gain in shorter transform length, but lose in having to do more FFTs, but that is a slight total save. We then lose in more expensive CRT.

[We now have two implementations of this algorithm, one by Tommy Färnqvist and one by Niels Möller.]

3. Work on short products. Our mullo and mulmid are probably K, but we lack mulhi.

Another possibility would be an optimized cube. In the basecase that should definitely be able to save cross-products in a similar fashion to squaring, but some investigation might be needed for how best to adapt the higher-order algorithms. Not sure whether cubing or further small powers have any particularly important uses though.

• Assembly routines

Write new and improve existing assembly routines. The tests/devel programs and the tune/speed.c and tune/many.pl programs are useful for testing and timing the routines you write. See the README files in those directories for more information.

Please make sure your new routines are fast for these three situations:

1. Small operands of less than, say, 10 limbs.
2. Medium size operands, that fit into the cache.
3. Huge operands that does not fit into the cache.

The most important routines are mpn_addmul_1, mpn_mul_basecase and mpn_sqr_basecase. The latter two don't exist for all machines, while mpn_addmul_1 exists for almost all machines.

Standard techniques for these routines are unrolling, software pipelining, and specialization for common operand values. For machines with poor integer multiplication, it is sometimes possible to remedy the situation using floating-point operations or SIMD operations such as MMX (x86) (x86), SSE (x86), VMX (PowerPC), VIS (Sparc).

Using floating-point operations is interesting but somewhat tricky. Since IEEE double has 53 bit of mantissa, one has to split the operands in small pieces, so that no intermediates are greater than 2^53. For 32-bit computers, splitting one operand into 16-bit pieces works. For 64-bit machines, one operand can be split into 21-bit pieces and the other into 32-bit pieces. (A 64-bit operand can be split into just three 21-bit pieces if one allows the split operands to be negative!)

• Faster sqrt

The current code uses divisions, which are reasonably fast, but it'd be possible to use only multiplications by computing 1/sqrt(A) using this iteration:

```				    2
x   = x  (3 − A x )/2
i+1	  i	    i  ```
The square root can then be computed like this:
```		     sqrt(A) = A x
n  ```

That final multiply might be the full size of the input (though it might only need the high half of that), so there may or may not be any speedup overall.

We should probably allow a special exponent-like parameter, to speed computations of a precise square root of a small number in mpf and mpfr.

• Nth root

Improve mpn_rootrem. The current code is not too bad, but its time complexity is a function of the input, while it is possible to make the average complexity a function of the output.

• Fat binaries

Add more functions to the set of fat functions.

The speed of multipliciaton is today highly dependent on combination functions like `addlsh1_n`. A fat binary will never use any such functions, since they are classified as optional. Ideally, we should use them, but making the current compile-time selections of optional functions become run-time selections for fat binaries.

If we make fat binaries work really well, we should move away frm tehe current configure scheme (at least by default) and instead include all code always.

• Exceptions

Some sort of scheme for exceptions handling would be desirable. Presently the only thing documented is that divide by zero in GMP functions provokes a deliberate machine divide by zero (on those systems where such a thing exists at least). The global `gmp_errno` is not actually documented, except for the old `gmp_randinit` function. Being currently just a plain global means it's not thread-safe.

The basic choices for exceptions are returning an error code or having a handler function to be called. The disadvantage of error returns is they have to be checked, leading to tedious and rarely executed code, and strictly speaking such a scheme wouldn't be source or binary compatible. The disadvantage of a handler function is that a `longjmp` or similar recovery from it may be difficult. A combination would be possible, for instance by allowing the handler to return an error code.

Divide-by-zero, sqrt-of-negative, and similar operand range errors can normally be detected at the start of functions, so exception handling would have a clean state. What's worth considering though is that the GMP function detecting the exception may have been called via some third party library or self contained application module, and hence have various bits of state to be cleaned up above it. It'd be highly desirable for an exceptions scheme to allow for such cleanups.

The C++ destructor mechanism could help with cleanups both internally and externally, but being a plain C library we don't want to depend on that.

A C++ `throw` might be a good optional extra exceptions mechanism, perhaps under a build option. For GCC `-fexceptions` will add the necessary frame information to plain C code, or GMP could be compiled as C++.

Out-of-memory exceptions are expected to be handled by the `mp_set_memory_functions` routines, rather than being a prospective part of divide-by-zero etc. Some similar considerations apply but what differs is that out-of-memory can arise deep within GMP internals. Even fundamental routines like `mpn_add_n` and `mpn_addmul_1` can use temporary memory (for instance on Cray vector systems). Allowing for an error code return would require an awful lot of checking internally. Perhaps it'd still be worthwhile, but it'd be a lot of changes and the extra code would probably be rather rarely executed in normal usages.

A `longjmp` recovery for out-of-memory will currently, in general, lead to memory leaks and may leave GMP variables operated on in inconsistent states. Maybe it'd be possible to record recovery information for use by the relevant allocate or reallocate function, but that too would be a lot of changes.

One scheme for out-of-memory would be to note that all GMP allocations go through the `mp_set_memory_functions` routines. So if the application has an intended `setjmp` recovery point it can record memory activity by GMP and abandon space allocated and variables initialized after that point. This might be as simple as directing the allocation functions to a separate pool, but in general would have the disadvantage of needing application-level bookkeeping on top of the normal system `malloc`. An advantage however is that it needs nothing from GMP itself and on that basis doesn't burden applications not needing recovery. Note that there's probably some details to be worked out here about reallocs of existing variables, and perhaps about copying or swapping between "permanent" and "temporary" variables.

Applications desiring a fine-grained error control, for instance a language interpreter, would very possibly not be well served by a scheme requiring `longjmp`. Wrapping every GMP function call with a `setjmp` would be very inconvenient.

Another option would be to let `mpz_t` etc hold a sort of NaN, a special value indicating an out-of-memory or other failure. This would be similar to NaNs in mpfr. Unfortunately such a scheme could only be used by programs prepared to handle such special values, since for instance a program waiting for some condition to be satisfied could become an infinite loop if it wasn't also watching for NaNs. The work to implement this would be significant too, lots of checking of inputs and intermediate results. And if `mpn` routines were to participate in this (which they would have to internally) a lot of new return values would need to be added, since of course there's no `mpz_t` etc structure for them to indicate failure in.

Stack overflow is another possible exception, but perhaps not one that can be easily detected in general. On i386 GNU/Linux for instance GCC normally doesn't generate stack probes for an `alloca`, but merely adjusts `%esp`. A big enough `alloca` can miss the stack redzone and hit arbitrary data. GMP stack usage is normally a function of operand size, which might be enough for some applications to know they'll be safe. Otherwise a fixed maximum usage can probably be obtained by building with `--enable-alloca=malloc-reentrant` (or `notreentrant`). Arranging the default to be `alloca` only on blocks up to a certain size and `malloc` thereafter might be a better approach and would have the advantage of not having calculations limited by available stack.

Actually recovering from stack overflow is of course another problem. It might be possible to catch a `SIGSEGV` in the stack redzone and do something in a `sigaltstack`, on systems which have that, but recovery might otherwise not be possible. This is worth bearing in mind because there's no point worrying about tight and careful out-of-memory recovery if an out-of-stack is fatal.

Operand overflow is another exception to be addressed. It's easy for instance to ask `mpz_pow_ui` for a result bigger than an `mpz_t` can possibly represent. Currently overflows in limb or byte count calculations will go undetected. Often they'll still end up asking the memory functions for blocks bigger than available memory, but that's by no means certain and results are unpredictable in general. It'd be desirable to tighten up such size calculations. Probably only selected routines would need checks, if it's assumed say that no input will be more than half of all memory and hence size additions like say `mpz_mul` won't overflow.

• Performance Tool

It'd be nice to have some sort of tool for getting an overview of performance. Clearly a great many things could be done, but some primary uses would be,

1. Checking speed variations between compilers.
2. Checking relative performance between systems or CPUs.

A combination of measuring some fundamental routines and some representative application routines might satisfy these.

The tune/time.c routines would be the easiest way to get good accurate measurements on lots of different systems. The high level `speed_measure` may or may not suit, but the basic `speed_starttime` and `speed_endtime` would cover lots of portability and accuracy questions.

• Using `restrict`

There might be some value in judicious use of C99 style `restrict` on various pointers, but this would need some careful thought about what it implies for the various operand overlaps permitted in GMP.

Rumour has it some pre-C99 compilers had `restrict`, but expressing tighter (or perhaps looser) requirements. Might be worth investigating that before using `restrict` unconditionally.

Loops are presumably where the greatest benefit would be had, by allowing the compiler to advance reads ahead of writes, perhaps as part of loop unrolling. However critical loops are generally coded in assembler, so there might not be very much to gain. And on Cray systems the explicit use of `_Pragma` gives an equivalent effect.

One thing to note is that Microsoft C headers (on ia64 at least) contain `__declspec(restrict)`, so a `#define` of `restrict` should be avoided. It might be wisest to setup a `gmp_restrict`.

• Factorial

Rewrite for simplicty and speed. Work is in progress.

• Binomial Coefficients

Rewrite for simplicty and speed. Work is in progress.

• Prime Testing

GMP is not really a number theory library and probably shouldn't have large amounts of code dedicated to sophisticated prime testing algorithms, but basic things well-implemented would suit. Tests offering certainty are probably all too big or too slow (or both!) to justify inclusion in the main library. Demo programs showing some possibilities would be good though.

The present "repetitions" argument to `mpz_probab_prime_p` is rather specific to the Miller-Rabin tests of the current implementation. Better would be some sort of parameter asking perhaps for a maximum chance 1/2^x of a probable prime in fact being composite. If applications follow the advice that the present reps gives 1/4^reps chance then perhaps such a change is unnecessary, but an explicitly described 1/2^x would allow for changes in the implementation or even for new proofs about the theory.

`mpz_probab_prime_p` always initializes a new `gmp_randstate_t` for randomized tests, which unfortunately means it's not really very random and in particular always runs the same tests for a given input. Perhaps a new interface could accept an rstate to use, so successive tests could increase confidence in the result.

`mpn_mod_34lsub1` is an obvious and easy improvement to the trial divisions. And since the various prime factors are constants, the remainder can be tested with something like

```#define MP_LIMB_DIVISIBLE_7_P(n) \
((n) * MODLIMB_INVERSE_7 <= MP_LIMB_T_MAX/7)
```
Which would help compilers that don't know how to optimize divisions by constants, and is even an improvement on current gcc 3.2 code. This technique works for any modulus, see Granlund and Montgomery "Division by Invariant Integers" section 9.

The trial divisions are done with primes generated and grouped at runtime. This could instead be a table of data, with pre-calculated inverses too. Storing deltas, ie. amounts to add, rather than actual primes would save space. `udiv_qrnnd_preinv` style inverses can be made to exist by adding dummy factors of 2 if necessary. Some thought needs to be given as to how big such a table should be, based on how much dividing would be profitable for what sort of size inputs. The data could be shared by the perfect power testing.

Jason Moxham points out that if a sqrt(-1) mod N exists then any factor of N must be == 1 mod 4, saving half the work in trial dividing. (If x^2==-1 mod N then for a prime factor p we have x^2==-1 mod p and so the jacobi symbol (-1/p)=1. But also (-1/p)=(-1)^((p-1)/2), hence must have p==1 mod 4.) But knowing whether sqrt(-1) mod N exists is not too easy. A strong pseudoprime test can reveal one, so perhaps such a test could be inserted part way though the dividing.

Jon Grantham "Frobenius Pseudoprimes" (www.pseudoprime.com) describes a quadratic pseudoprime test taking about 3x longer than a plain test, but with only a 1/7710 chance of error (whereas 3 plain Miller-Rabin tests would offer only (1/4)^3 == 1/64). Such a test needs completely random parameters to satisfy the theory, though single-limb values would run faster. It's probably best to do at least one plain Miller-Rabin before any quadratic tests, since that can identify composites in less total time.

Some thought needs to be given to the structure of which tests (trial division, Miller-Rabin, quadratic) and how many are done, based on what sort of inputs we expect, with a view to minimizing average time.

It might be a good idea to break out subroutines for the various tests, so that an application can combine them in ways it prefers, if sensible defaults in `mpz_probab_prime_p` don't suit. In particular this would let applications skip tests it knew would be unprofitable, like trial dividing when an input is already known to have no small factors.

For small inputs, combinations of theory and explicit search make it relatively easy to offer certainty. For instance numbers up to 2^32 could be handled with a strong pseudoprime test and table lookup. But it's rather doubtful whether a smallnum prime test belongs in a bignum library. Perhaps if it had other internal uses.

An `mpz_nthprime` might be cute, but is almost certainly impractical for anything but small n.

• Intra-Library Calls

On various systems, calls within libgmp still go through the PLT, TOC or other mechanism, which makes the code bigger and slower than it needs to be.

The theory would be to have all GMP intra-library calls resolved directly to the routines in the library. An application wouldn't be able to replace a routine, the way it can normally, but there seems no good reason to do that, in normal circumstances.

The `visibility` attribute in recent gcc is good for this, because it lets gcc omit unnecessary GOT pointer setups or whatever if it finds all calls are local and there's no global data references. Documented entrypoints would be `protected`, and purely internal things not wanted by test programs or anything can be `internal`.

Unfortunately, on i386 it seems `protected` ends up causing text segment relocations within libgmp.so, meaning the library code can't be shared between processes, defeating the purpose of a shared library. Perhaps this is just a gremlin in binutils (debian packaged 2.13.90.0.16-1).

The linker can be told directly (with a link script, or options) to do the same sort of thing. This doesn't change the code emitted by gcc of course, but it does mean calls are resolved directly to their targets, avoiding a PLT entry.

Keeping symbols private to libgmp.so is probably a good thing in general too, to stop anyone even attempting to access them. But some undocumented things will need or want to be kept visible, for use by mpfr, or the test and tune programs. Libtool has a standard option for selecting public symbols (used now for libmp).

• Math functions for the mpf layer

Implement the functions of math.h for the GMP mpf layer! Check the book "Pi and the AGM" by Borwein and Borwein for ideas how to do this. These functions are desirable: acos, acosh, asin, asinh, atan, atanh, atan2, cos, cosh, exp, log, log10, pow, sin, sinh, tan, tanh.

Note that the mpfr functions already provide these functions, and that we usually recommend new programs to use mpfr instead of mpf.