mpfr / gmp

Arbitrary Precision Arithmetic: mpir is a fork of gmp (gnu multiple precision), and mpfr is a floating-point extension of that.

It has been wrapped/is available for many programs languages, including Ada, C++, Fortran, Haskell, Java, Julia, Lisp, .NET, OCaml, Perl, PHP, Pike, Prolog, Python, R, Racket, Ruby, Rust, Scheme, and now Phix - which means there should be no shortage of examples.

I have seen some performance improvements in my time, but this absolutely has to take the (chocolate covered triple-choc) biscuit.
For demo\rosetta\Ackermann, translating the Go code to calculate ack(3,1e6), which has over 300,000 digits:
Using bigatom.e (as originally submitted by cargoan) took so long that I just had to kill it....
Using bigint.e (wot I wrote meself): 5 minutes to calculate and then 32 minutes to print (gulp).
Using gmp: 0.1s for the lot! (Do I feel humiliated and then some or what?!)

Included in the distribution as builtins\mpfr.e (not an autoinclude), but you have to download the dlls from PCAN or perhaps http://www.atelierweb.com/mpir-and-mpfr
You should get a readable and straightforward error message when the required dlls cannot be opened.
It is entirely up to you whether to install (simple copy) them to system32/syswow64, builtins, or the application directory.
A future release of phix may offer to download and install them (to builtins) as part of installation, however you would still be responsible for shipping them along with a finished application.
Let me know if any newer pre-built dlls are uploaded to atelierweb, or elsewhere.

Version: mpfr: 3.1.5, mpir:2.7.2, two dlls, 733K (32 bit, 308K zipped), 936K (64 bit, 373K zipped)
Version 3.1.0-p3 was already installed my Ubuntu box, though I also needed to install mpfr-dev, in order to compile the recommended version check, which was just the one click via the ubuntu software centre, at least after typing "mpfr" into the search box and picking the one that I needed.

Be warned: mpfr/mpir/gmp are built for speed and are generally very unforgiving, hanging or terminating silently without any clues as to what went wrong - good old-fashioned C for ya!
Making (for example) mpfr_clear() private/internal, and wrapping the C mpfr_init() functions and hence effectively making them non-optional is the least/best I could do to help.

Use the mpfr_xxx() routines for your ridiculously accurate floats (with ridiculously big exponents).
Use the mpz_xxx() routines for your ridiculously big (and perfectly exact) integers.

Once you understand how to create and display mpfr/mpz variables, the rest is all pretty straightforward.

Example:

include mpfr.e

mpz n = mpz_init()
for e = 1 to 3 do
    mpz_ui_pow_ui(n, 10, e*30)
    mpz_sub_ui(n, n, 1)
    printf(1, "10^%d-1 = %s\n",{e*30,mpz_get_str(n,comma_fill:=true)})
end for
n = mpz_clear(n)

-- output:
--  10^30-1 = 999,999,999,999,999,999,999,999,999,999
--  10^60-1 = 999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999
--  10^90-1 = 999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999
For comparision, here is the (accurate to ˜15/17dp) behaviour of the standard builtin atom: (Manually re-aligned. Aside: 10^60 has seven more digits of accuracy than expected on 32 bit, no idea why, or even how.)

for e = 1 to 3 do
    printf(1, "10^%d-1 = %,d\n",{e*30,power(10,e*30)-1})
end for

-- output (32-bit):
--  10^30-1 = 1,000,000,000,000,000,424,684,240,284,426
--  10^60-1 = 1,000,000,000,000,000,000,000,082,442,268,486,026,840,024,008,840,604,006,400,824
--  10^90-1 =   999,999,999,999,999,686,426,868,868,000,042,200,008,468,264,024,420,280,668,688,280,420,208,406,626,688,068,446,402,024
-- output (64-bit):
--  10^30-1 =   999,999,999,999,999,998,264,846,264,264
--  10^60-1 =   999,999,999,999,999,996,186,866,828,000,402,842,426,462,280,408,402,286,040,882
--  10^90-1 = 1,000,000,000,000,000,004,244,424,008,082,284,008,804,242,048,002,866,280,462,022,042,286,202,060,886,602,048,286,628,466
Many operations are performed in-situ via procedure calls, eg mpfr_mul(x,y,z) sets x to y*z.
Source and targets may overlay, eg mpfr_mul(x,x,x) etc is perfectly fine and sets x to x*x.

Since phix does not have unsigned integers, if an _ui routine has an _si variant it should be used instead. It also makes sense to utilise the standard integer typechecking within the mpfr.e wrapper, hence some of the _ui routines are syntax coloured illegal/red and linked to their _si counterparts, and where no such exists in C the _ui name is kept, but with fatal errors for args<0.
In the rare case some code actually uses an unsigned value above the normal 30 bit integer limit of 1GB, it must be changed to use xxx_d or xxx_str.
Also take heed of the -1GB..+1GB vs 0..1GB below (obviously not technically accurate for 64-bit, but still a necessity for any 32 and 64 bit compatible code).
As ever, this document is just a quick cribsheet: for specific details refer to the mpfr.e source code and/or the official mpfr/mpir/gmp documentation.

Types mpfr and mpz and randstate are optional, you can just use atom if preferred, though when debug_types is true (it is) these also perform some extra anti-corruption checks.

mpfr (floating point)

string res =
mpir_open_dll(string dll_name="", bool mpir_only=false) - optional: returns the name of a dll/so file it cannot open, or "" on success.
If only using mpz routines, you can get by with only opening/distributing the mpir dll(s), and should set the optional flag to true, and in that case mpir_open_dll() becomes non-optional. Conversely, attempts to load the mpfr dll/so, without loading the mpir one first, fail at the OS level.
If dll/so are missing, and you have not used this, any and all of the following routines below will terminate with a fatal error (including mpz_’s).
sequence res =
mpfr_get_versions(bool bAsNumSeq=false) - yields a 3-element sequence of three versions (mpfr, mpir, and gmp), eg {"3.1.5","2.7.2","5.1.3"} or {{3,1,5},{2,7,2},{5,1,3}}.
Note that on windows 64 bit, the mpfr.e wrapper assumes mpir is version 2.6.0 or later, with 64 bit integer arguments for _ui and _si routines, but does not explicitly check that, which may cause subtle problems should you obtain and use an older dll build (not from me).
integer precision =
mpfr_get_default_prec() - yields the default floating point precision in bits, initially 53 to match IEEE-754.
mpfr_set_default_prec(integer precision) - set the default precision in bits.
You can determine the required precision using log2(), or perhaps, probably in some prior test run:
    string nines = repeat('9',<number_of_decimal_digits_required>)
    mpz ndp = mpz_init(nines)
    integer precision = mpz_sizeinbase(ndp,2)
integer rounding =
mpfr_get_default_rounding_mode() - yields MPFR_RNDx where x is N/Z/U/D/A (0..4, initially N==0) for nearest/toward zero/+inf/-inf/away from zero.
mpfr_set_default_rounding_mode(integer rounding) - sets the rounding mode, as above.
mpfr res =
mpfr_init(object v=0, integer precision=default_precision, rounding=default_rounding) - Initialise to integer/atom/string.
Covers the C functions mpfr_init, mpfr_init2, mpfr_init_set_si, mpfr_init_set_d, and mpfr_init_set_str (note the latter three can take a precision in mpfr.e, and that comes before the rounding mode), however it does not cover mpfr_inits or mpfr_inits2, and it must not be used for mpfr_init_set(mpfr), ie init from another existing mpfr, see next.
Note that (as per v=0) the default is to initialise to 0, unlike the C api which creates variables set to nan (not a number).
mpfr res = mpfr_init_set(mpfr src, integer rounding=default_rounding) - Initialise to some prior mpfr, eg (mpfr x = mpfr_init(...);) mpfr y = mpfr_init_set(x).
Passing an mpfr directly to mpfr_init() is a similar error to using pPtr when peek4u(pPtr) is required, and likewise there is no practical way to guard against that.
integer precision =
mpfr_get_prec(mpfr x) - retrieve the current precision of a specific variable.
mpfr_set_prec(mpfr x, integer precision) - precision is the number of bits required for the mantissa, see mpfr_set_default_prec above.
mpfr_set(mpfr tgt, src, integer rounding=default_rounding) - set x from an existing mpfr variable.
mpfr_set_si(mpfr x, integer v, rounding=default_rounding) - set x from a (machine-word-sized) integer.
In mpfr.e, v is (deliberately) limited to +/-1 billion, rather than trying to extend the range a bit by using an atom parameter.
mpfr_set_d(mpfr x, atom v, integer rounding=default_rounding) - set x from a normal phix atom.
mpfr_set_str(mpfr x, string s, integer base=0, rounding=default_rounding) - set x from a string.
If base is zero, bases 2 and 16 can be auto-detected (after a leading +/-) by '0b', '0B', '0x' or '0X', otherwise base 10 is assumed.
The exponent prefix can be 'e' or 'E' for bases up to 10, or '@' in any base, or 'p' or 'P' for bases 2 and 16.
The value of an exponent is always written in base 10. A '.' is always accepted as a decimal point (as well as the current locale).
mpfr_const_pi(mpfr x, integer rounding=default_rounding) - x := PI
As precise as needed/specified in some prior mpfr_set_prec(x) or the x = mpfr_init() call.
mpfr_clear[s]() - use mpfr_free() instead (which clears all memory alocated by mpfr/mpir.dll and mpfr.e).
object x =
mpfr_free(object x) - Clear and deallocate any variables created from mpfr_init(), eg x = mpfr_free(x) or {y,z} = mpfr_free({y,z})
NB: lhs==rhs intentionally. Should you forget to invoke this routine, mpfr.e will do so automatically.
mpfr_free_str() - invoked automatically within mpfr_get_str() and hence not publically available (there is no need for it to be).
sequence res =
mpfr_get_str(mpfr x, integer base=10, n=0, rounding=default_rounding) -- print x in the specified base to precision n(0=full).
res is {string digits, integer exponent}. Does not print a decimal point; mpfr_sprintf() [see next] may be a better choice.
Also note the parameters are quite different, C: mpfr_get_str(?, ?, b, n, x, r) ==> phix: mpfr_get_str(x, b, n, r).
For the obvious reason of a FILE* argument, use this or mpfr_sprintf instead of mpfr_out_str.
string res =
mpfr_sprintf(string fmt, mpfr_or_mpz x) - formatted print using 'R', eg mpfr_sprintf(x,"%.1000Rf") prints x to 1000 d.p.
Code using mpfr_Xprintf() where X!="s" must be changed to use mpfr_sprintf, and any C-style length modifiers (eg h, ll) removed.
Since x is typechecked to be of type mpfr_or_mpz, only a single 'R' or 'Z' conversion specification is supported.
An 'R' can be followed by an optional rounding mode (one of NZUDA), and then it must be followed by one of abefgAEFG.
A 'Z' should be followed by one of "dxX", or possibly (untested) one of "aceEfigGo".
Currently the only other routine of similar note that has been wrapped is mpz_get_str(), which can be used as an alternative for 'Z' conversion specifications.
mpfr_printf(integer fn, string fmt, mpfr_or_mpz x) - formatted print to file.
Note this differs from the C function of the same name, in having a file number as the first parameter (like C’s mpfr_fprintf), plus the phix version only accepts/permits a single mpfr or mpz variable at a time.
mpfr_add(mpfr rop, op1, op2, integer rounding=default_rounding) - rop := op1+op2 with specified rounding
mpfr_add_si(mpfr rop, op1, integer op2, integer rounding=default_rounding) - "" except op2 is a phix integer, -1GB..+1GB
mpfr_sub(mpfr rop, op1, op2, integer rounding=default_rounding) - rop := op1-op2 with specified rounding
mpfr_sub_si(mpfr rop, op1, integer op2, integer rounding=default_rounding) - "" except op2 is a phix integer, -1GB..+1GB
mpfr_mul(mpfr rop, op1, op2, integer rounding=default_rounding) - rop := op1*op2 with specified rounding
mpfr_mul_si(mpfr rop1, op1, integer op2, integer rounding=default_rounding) - "" except op2 is a phix integer, -1GB..+1GB
mpfr_div(mpfr rop, op1, op2, integer rounding=default_rounding) - rop := op1/op2 with specified rounding
mpfr_div_si(mpfr rop, op1, integer op2, integer rounding=default_rounding) - "" except op2 is a phix integer, -1GB..+1GB
mpfr_sqr(mpfr rop, op, integer rounding=default_rounding) - rop := op squared with specified rounding
mpfr_sin(mpfr rop, op, integer rounding=default_rounding) - rop := sin(op) with specified rounding
atom res = mpfr_get_si(mpfr op, integer rounding=default_rounding) -- res := op as a machine-word-sized signed integer.

mpz (integer)

mpz res =
mpz_init(object v=0, integer bitcount=0) - Initialise. v can be integer/atom/string.
Covers the C functions mpz_init, mpz_init2, mpz_init_set_d, mpz_init_set_si, and mpz_init_set_str (note the latter three can take a bitcount in mpfr.e), however it must not be used for mpz_init_set(mpz), ie init from another existing mpz, see next.
(Aside: unlike mpfr_init(), the defaulted v=0 matches the C api.)
The bitcount is only the initial space, ie/eg mpz_init(1) happily stores it in one dword, a later store of 2^16000 reallocates it to 16000 bits (about/whatever 500 dwords). A non-zero bitcount makes it possible to avoid such reallocations if a maximum size is known in advance.
You could use mpz_sizeinbase(x,2), in some prior test run, to determine a suitable (initial) bitcount (and see whether or not it actually makes the thing any faster).
mpz res = mpz_init_set(mpz src) - usage eg (mpz x = mpz_init(...);) mpz y = mpz_init_set(x).
Passing an mpz directly to mpz_init() is a similar error to using pPtr when peek4u(pPtr) is required, and likewise there is no practical way to guard against that.
object x =
mpz_clear(object x) - Clear and deallocate any variables created using mpz_init[_set](), eg x = mpz_clear(x) or {y,z} = mpz_clear({y,z}).
NB: lhs==rhs intentionally. Also covers the C mpz_clears() routine. Should you forget to invoke this routine, mpfr.e will do so automatically.
mpz_set(mpz rop, op) -- rop := op. See also mpz_init_set.
mpz_set_si(mpz rop, integer op) -- "" except op is a phix integer, -1GB..+1GB.
mpz_set_d(mpz rop, atom op) -- "" except op is a phix atom.
mpz_set_str(mpz rop, string s, integer base=10) -- set x from s in the specified base. See also mpz_init.
mpz_add(mpz rop, op1, op2) - rop := op1 + op2
mpz_add_ui(mpz rop, op1, integer op2) - "" except op2 is a phix integer, 0..1GB
mpz_sub(mpz rop, op1, op2) - rop := op1 - op2
mpz_sub_ui(mpz rop, op1, integer op2) - "" except op2 is a phix integer, 0..1GB
mpz_mul(mpz rop, op1, op2) - rop := op1 * op2
mpz_mul_si(mpz rop, op1, integer op2) - "" except op2 is a phix integer, -1GB..+1GB
mpz_mul_2exp(mpz rop, op1, integer op2) - rop := op1 * 2^op2. This operation can also be defined as a left shift by op2 bits.
mpz_fdiv_q_2exp(mpz q, n, integer bit_count) - bitwise right shift, arithmetic if n -ve.
mpz_mod(mpz r, n, d) - r := mod(n,d)
mpz_mod_ui(mpz r, n, integer d) - "" except op2 is a phix integer, 0..1GB
integer res =
mpz_cmp(mpz op1, op2) - Compare op1 and op2. Return a positive value if op1 > op2, zero if op1 = op2, or a negative value if op1 < op2.
integer res = mpz_cmp_si(mpz op1, integer op2) - "", except op2 is a phix integer, -1GB..+1GB
integer res = mpz_tstbit(mpz op, integer bit_index) - Test bit bit_index in op and return 0 or 1 accordingly.
integer res = mpz_scan0(mpz op, integer starting_bit) - Find first 0 in op >= starting_bit.
integer res = mpz_scan1(mpz op, integer starting_bit) - Find first 1 in op >= starting_bit.
mpz_powm(mpz rop, base, exponent, modulus) - rop := mod(base^exponent,modulus)
mpz_powm_ui(mpz rop, base, integer exponent, mpz modulus) - "", except exponent and modulus are a phix integers, 0..+1GB
mpz_ui_pow_ui(mpz rop, integer base, exponent) - rop := base^exponent. The case 0^0 yields 1.
bool res = mpz_fits_slong_p(mpz op) - Return non-zero iff the value of op fits in a (signed) long, otherwise, return zero.
bool res = mpz_fits_ulong_p(mpz op) - Return non-zero iff the value of op fits in an unsigned long, otherwise, return zero.
atom res = mpz_get_si(mpz op) - Return op as an machine-word-sized signed integer. Use mpz_fits_slong_p() to find out if the value will fit/be meaningful.
atom res = mpz_get_ui(mpz op) - Return op as an machine-word-sized unsigned integer. Use mpz_fits_ulong_p() to find out if the value will fit/be meaningful.
integer res = mpz_sizeinbase(mpz op, integer base) - Return the size of op measured in number of digits in the given base.
string res =
mpz_get_str(mpz x, integer base=10, bool comma_fill=false) - Return x as a string in the specified base.
Note that mpz_free_string() is taken care of automatically for you by mpfr.e.
bool res =
mpz_probable_prime_p(mpz n, randstate state, integer prob=5, div=0)
Determine whether n is a probable prime with the chance of error being at most 1 in 2^prob.
The return value is 1 if n is probably prime, or 0 if n is definitely composite.
This function does some trial divisions to speed up the average case, then some probabilistic primality tests to achieve the desired level of error.
div can be used to inform the function that trial division up to div has already been performed on n and so n has NO divisors <= div.
(The default of 0 informs the function that no trial division has been done.)
The variable state must have been initialized by calling one of the gmp_randinit functions.
This function interface is preliminary and may change in the future.

random integers

randstate state = gmp_randinit_mt() - Initialize state for a Mersenne Twister algorithm, and invokes gmp_randseed on it (with a NULL mpz_seed).
gmp_randseed(randstate state, atom mpz_seed=NULL) - Set an initial seed value into state.
If mpz_seed is NULL then a random 200-digit seed is generated internally (and properly disposed of after use), otherwise mpz_seed must be a valid/init mpz variable (and absolutely not a phix integer/atom/string).
Just like set_rand(), using a fixed value gives repeatable results.
randstate state = gmp_randclear(randstate state) - Free all memory occupied by state. Should you forget to invoke this routine, mpfr.e will do so automatically.
mpz_urandomm(mpz rop, randstate state, mpz n) - Generate a uniform random integer in the range 0 to n-1 inclusive, and store it in rop.

More routines (many more) will be added as needed - there are several hundred of them and, at least in my opinion, it is just simply not worthwhile without something that actually uses and therefore tests them.

Note that mpfr_clear() is not publically available, use mpfr_free() instead (which frees both mpfr.dll and mpfr.e allocated memory). Since mpfr_clear() has a different calling convention (to mpfr_free), it is syntax coloured illegal (and linked here). In contrast, mpr_free_str() is invoked internally for you so when translating code you can just delete that statement.
Also note that all mpf_xxx functions have deliberately been expunged in favour of mpfr_ ones.