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, 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, and these days (quite unlike the first couple of versions) I very rarely suffer a sudden unexplained crash and almost always get a proper human-readable error message.

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.
Use the mpq_xxx() routines for rational numbers (aka fractions, ie n/m’ths, where n and m are mpz internally).

There are no hard-and-fast limits other than available memory; you should have no difficulty constructing a value that requires over 1GB of memory in the most compact possible form, as compared the the 4/8/10 bytes of native phix integers and atoms, though you might have trouble constructing a 2nd/3rd/4th that big at the same time. Printing numbers with more than about 500,000 digits is however neither very fast nor realistically speaking very sensible. See demo\rosetta\Fermat.exw for an example of how much further along you can get, provided you give up printing.

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

Note however that mpz/mpfr/mpq/randstate declare reference variables, as opposed to the copy-on-write semantics of standard phix types, eg after mpz a=mpz_init(3),b=a; mpz_add_ui(b,b,1) both a and b are 4, whereas after mpz a=mpz_init(3), b=mpz_init_set(a); mpz_add_ui(b,b,1) then a is 3 and b is 4. In short you need an mpz_init or similar for every distinct value that needs to be held simultaneously.

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_free(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 overlap, 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, and/or search demo/rosetta for examples.

Types mpfr, mpq, 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. They all allow NULL.

mpz (integer)

string res =
mpir_open_dll(string dll_name="", boolean 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.
sequence res =
mpir_get_versions(boolean 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).
If mpir_open_dll() was invoked with mpir_only=true, the first element of the result will be "mpfr_dll==NULL".
mpz res =
mpz_init(object v=0, integer bitcount=0) - Initialise. v can be integer/atom/string, but not mpz.
Note that a string value in v may not use exponentiation, such as "1e200" (as per the C api, sorry).
Unlike the nans of C’s mpfr_init(), the defaulted v=0 of mpz_init() matches the C api behaviour.
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 v may not be another existing mpz, for that use mpz_init_set() as described below.
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).
Invoking res = mpz_free(res) when no longer needed is recommended, but occurs automatically if forgotten.
sequence res =
mpz_inits(integer n, object v=0) - initialise n variables to v. (ditto mpz_free)
Note the method signature differs significantly from the C function of the same name.
Can also be coded as eg mpz {x,y,z} = mpz_inits(3).
v may be an integer, atom, string, or a sequence of length n of said.
mpz res =
mpz_init_set(mpz src) - Initialise to some prior mpz, eg (mpz x = mpz_init(...);) mpz y = mpz_init_set(x). (ditto mpz_free)
Passing an mpz directly to mpz_init() is a similar error to using pDword when peek4u(pDword) is required, and likewise there is no practical way to guard against that, you just get meaningless results. In fact mpz_init_set() can and does test it got an mpz, it is mpz_init() that simply cannot perform the inverse useful test that makes sure it has not been passed an mpz. (Even if it maintained a dirty great list, what’s to say an actual value will not coincide with an allocated memory address?)
object x =
mpz_free(object x) - Clear and deallocate any variables created using mpz_init[s|_set](), eg x = mpz_free(x) or {y,z} = mpz_free({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 (when the refcount drops to zero).
mpz_clear() is not publicly available, use mpz_free() instead, which releases both mpir.dll and mpfr.e allocated memory.
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=0) -- set x from s in the specified base. See also mpz_init.
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.
mpz_import(mpz rop, integer count, order, size, endian, nails, atom_string op) -- set rop from raw binary data.
The parameters specify the format of the data.
count: many words are read, each size bytes.
order can be 1 for most significant word first or -1 for least significant first.
Within each word endian can be 1 for most significant byte first, -1 for least, or 0 for the native endianness of the host CPU.
The most significant nails bits of each word are skipped, this can be 0 to use the full words.
There is no sign taken from the data, rop will simply be a positive integer.
An application can handle any sign itself, and apply it for instance with mpz_neg.
There are no data alignment restrictions on op, any address is allowed.
integer count =
mpz_export(atom pMem, integer order, size, endian, nails, mpz op) -- Fill pMem with word data from op.
The parameters specify the format of the data produced.
Each word will be size bytes and order can be 1 for most significant word first or -1 for least significant first.
Within each word endian can be 1 for most significant byte first, -1 for least, or 0 for the native endianness of the host CPU.
The most significant nails bits of each word are unused and set to zero, this can be 0 to produce full words.
The number of words produced is returned, pMem must have enough space for the data (see below).
(The second parameter of the C function, size_t *countp, is handled internally by mpfr.e)
If op is non-zero then the most significant word produced will be non-zero.
If op is zero then the count returned will be zero and nothing written to rop.
The sign of op is ignored, just the absolute value is exported. An application can use mpz_sign to get the sign and handle it as desired.
There are no data alignment restrictions on rop, any address is allowed.
The required space can be determined with a calculation like the following.
Since mpz_sizeinbase always returns at least 1, count here will be at least one, though if op is zero no space at all is actually needed (or written).
numb = 8*size - nail
count = floor((mpz_sizeinbase (op, 2) + numb-1) / numb)
pMem = allocate(count * size)
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_add_si(mpz rop, op1, integer op2) - (an extra) "" except op2 can be +/-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_ui_sub(mpz rop, integer op1, mpz op2) - "" except op1/op2 types are switched.
mpz_sub_si(mpz rop, op1, integer op2) - (an extra) as mpz_sub_ui() except op2 can be +/-1GB.
mpz_si_sub(mpz rop, integer op1, mpz op2) - (an extra) as mpz_ui_sub() except op1 can be +/-1GB.
mpz_abs(mpz rop, op) - rop := abs(op)
mpz_neg(mpz rop, op) - rop := -op
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_d(mpz rop, op1, atom op2) - (an extra) as mpz_mul_si() except op2 is a phix atom
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(mpz q, n, d) - q := floor(n/d)
mpz_fdiv_r(mpz r, n, d) - r := remainder(n,d)
mpz_fdiv_qr(mpz q, r, n, d) - {q,r} := {floor(n/d),remainder(n,d)}
integer res = mpz_fdiv_q_ui(mpz q, n, integer d) - {q,res} := {floor(n/d),remainder(n,d)}
mpz_fdiv_q_2exp(mpz q, n, integer b) - bitwise right shift, arithmetic if n -ve.
mpz_tdiv_q_2exp(mpz q, n, integer b) - q := trunc(n/2^b), rounds q towards zero
mpz_tdiv_r_2exp(mpz r, n, integer b) - r := remainder(n,2^b), r will have the same sign as n
mpz_cdiv_q(mpz q, n, d) - q := ceil(n/d)
For positive n mpz_tdiv_q_2exp is a simple bitwise right shift.
For negative n mpz_tdiv_q_2exp effectively treats n as sign and magnitude. [untested...]
fdiv rounds q down towards -inf, and r will have the same sign as d. The f stands for "floor".
tdiv rounds q towards zero, and r will have the same sign as n. The t stands for "truncate".
cdiv rounds q up towards +inf, and r will have the opposite sign to d. The c stands for "ceil".
In all cases q and r will satisfy n = q(d) + r, and r will satisfy 0 <= |r| < |d|, where d is 2^b.
boolean res =
mpz_divisible_p(mpz n, d) - returns non-zero if n is exactly divisible by d.
n is divisible by d if there exists an integer q satisfying n = qd.
Unlike the other division functions, d = 0 is accepted and following the rule it can be seen that only 0 is considered divisible by 0.
boolean res = mpz_divisible_ui_p(mpz n, integer d) - "" except d is a phix integer, 0..1GB
boolean res = mpz_divisible_2exp_p(mpz n, integer b) - "" except tests if n is divisible by 2 b
integer res =
mpz_remove(mpz rop, op, f) - removes all occurrences of the factor f from op and stores the result in rop.
Typically only invoked after mpz_divisible_p() returns non-zero. The return value is how many occurrences were removed.
integer res =
mpz_fdiv_ui(mpz n, integer d) - returns mod(n,d) - n and d remain unaltered, d is a phix integer, 0..1GB.
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
(In practice mpz_mod_ui() utilises the mpz_fdiv_r_ui() C entry point.)
mpz_pow_ui(mpz rop, base, integer exponent) - rop := base^exponent, where exponent is a phix integer, 0..+1GB.
mpz_ui_pow_ui(mpz rop, integer base, exponent) - "", except base is also a phix integer. The case 0^0 yields 1.
mpz_powm(mpz rop, base, exponent, modulus) - rop := mod(base^exponent,modulus)
mpz_powm_ui(mpz rop, base, integer exponent, mpz modulus) - "", except exponent is a phix integer, 0..+1GB
bool bExact =
mpz_root(mpz rop, op, integer n) - rop:=trunc(power(op,1/n)), returns (rop===op^n).
Set rop to the truncated integer part of the nth root of op.
Return true if the computation was exact, i.e. rop===op^n.
mpz_sqrt(mpz rop, op) - rop := floor(sqrt(op))
mpz_sqrtrem(mpz rop1, rop2, op) - {rop1,rop2} := {floor(sqrt(op)),op-rop1^2}
rop2 will be zero if op is a perfect square. If rop1 and rop2 are the same variable, the results are undefined.
mpz_fib_ui(mpz fn, integer n) - fn := fibonacci(n), where n is a phix integer, 0..+1GB
mpz_fib2_ui(mpz fn, fnsub1, integer n) - as "", plus fnsub1 := fibonacci(n-1) (and must have the prev value on input).
See demo\rosetta\fibonacci.exw: note that mpz_fib_ui() is intended for isolated values only; while an isolated mpz_fib_ui(4784969) takes just 0.1s, invoking mpz_fib_ui(1..4784969) would take about 3½ days(!!) whereas the version in that file gets the same job done in about 2mins 40s, averaging ˜30,000/s - not bad for million-digit numbers!
mpz_gcd(mpz rop, op1, op2) - Set rop to the greatest common divisor of op1 and op2.
The result is always positive even if one or both input operands are negative.
atom res =
mpz_gcd_ui(mpz rop, op1, integer op2) - Compute the greatest common divisor of op1 and op2.
If rop is not NULL, store the result there.
If the result is small enough to fit in an mpir_ui, it is returned.
If the result does not fit, 0 is returned, and the result is equal to the argument op1.
Note that the result will always fit if op2 is non-zero.
mpz_lcm(mpz rop, op1, op2) - Set rop to the least common multiple of op1 and op2.
mpz_lcm_ui(mpz rop, op1, integer op2) - "" except op2 is a phix integer, 0..1GB.
rop is always positive, irrespective of the signs of op1 and op2. rop will be zero if either op1 or op2 is zero.
mpz_fac_ui(mpz rop, integer n) - Set rop to the factorial of n.
mpz res =
mpz_binom(integer n, k) - calculate n choose k, ie n!/((n-k)!*k!)
Equivalent, for small n and k, to choose(n,k)
integer res =
mpz_cmp(mpz op1, op2) - Compare op1 and op2. Return +1 if op1>op2, 0 if op1=op2, or -1 if op1<op2.
Note the mpfr.e wrapper explicitly converts the C +ve/0/-ve result to +1/0/-1 using sign().
integer res = mpz_cmp_si(mpz op1, integer op2) - "", except op2 is a phix integer, -1GB..+1GB
integer res = mpz_sign(mpz op1) - -1: op1 -ve, 0: op1=0, +1: op1 +ve. Compatibility shim for the C macro mpz_sgn
boolean res = mpz_odd(mpz op1) - returns true is op1 is odd. Compatibility shim for the C macro mpz_odd_p
boolean res = mpz_even(mpz op1) - returns true is op1 is even. Compatibility shim for the C macro mpz_even_p
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.
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, boolean comma_fill=false) - Return x as a string in the specified base (2..62).
Note that mpz_free_string() is taken care of automatically for you by mpfr.e.
boolean 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.

extras

The following routines are mpfr.e-specific and not part of the C api, hence there is probably not much point in googling them.
Some were largely written for simple convenience, rather than absolute performance, although I would say that mpz_prime_factors is (I think) exceptionally efficient at what it does best, unless/until you push it too far, as detailed below. Also, some routines (the integer/atom ones) are deliberately renamed because they have slightly different upper limits.

mpz_add_si(mpz rop, op1, integer op2) - a simple wrapper that invokes either mpz_add_ui, or mpz_sub_ui with -op2, and will crash if passed -(maxint+1).
mpz_sub_si(mpz rop, op1, integer op2) - a simple wrapper that invokes either mpz_sub_ui, or mpz_add_ui with -op2, and will crash if passed -(maxint+1).
mpz_si_sub(mpz rop, integer op1, mpz op2) - a simple wrapper that invokes mpz_sub_si(rop,op2,op1) then mpz_mul_si(rop,rop,-1).
mpz_mul_d(mpz rop, op1, atom op2) - as mpz_mul_si() except op2 is a phix atom
mpz/integer res =
mpz_min(sequence s, boolean return_index=false) - return the smallest element of s, or it’s index.
s must be a flat and non-empty sequence of mpz.
This routine is similar to minsq() [as opposed to min()], but applies mpz_cmp() internally.
The result is a shared reference, modifying it directly (ie w/o an mpz_[init_]set() call) will modify the corresponding element of s.
The corresponding mpz_max() routine is identical, but obviously returns the largest element/index of s.
boolean res = mpz_fits_integer(mpz op) - Return non-zero if the value of op fits in a (signed) phix integer, otherwise, return zero.
Note that the C routines mpz_fits_slong_p() and mpz_fits_ulong_p() are (deliberately) not available in mpfr.e, and
this routine tests for slightly smaller limits, in fact (on 32 bit) +/-#3FFFFFFF with -#40000000 being deemed out-of-range.
boolean res = mpz_fits_atom(mpz op, boolean tztrim=false) - Return non-zero if the value of op fits in a (signed) phix atom, otherwise, return zero.
Note that (on 32 bit) 9,007,199,254,740,992 is deemed out-of-range, since that is the first value that "accidentally" fits, by ending in a binary 0, that is when tztrim is false.
When tztrim is true it returns true for that and 9,007,199,254,740,994 but obviously false for 9,007,199,254,740,993.
Obviously that 994 is technically out of range, but would not be altered by a round trip into and back out of an atom, which is what that tztrim (trim trailing binary zeroes) parameter effectively implements.
While on 32-bit there is a huge difference between 31-bit integers and 53-bit atoms, on 64-bit there is only 1 bit of difference between mpz_fits_integer (63 bits) and mpz_fits_atom (64 bits), at least when tztrim is false.
integer res =
mpz_get_integer(mpz op) - Return op as a phix integer. Use mpz_fits_integer() to find out if the value will fit/be meaningful.
Note there are (deliberately) no mpz_get_si() or mpz_get_ui() functions in mpfr.e, use this instead and be aware that a phix integer is one bit smaller than a C signed long.
atom res = mpz_get_atom(mpz op) - Return op as a phix atom. Use mpz_fits_atom() to find out if the value will fit/be meaningful.
type mpz_or_string - I trust this is self explanatory, and obviously integer/atom are not allowed.
sequence res =
mpz_prime_factors(mpz_or_string s, integer maxprime=100) - attempt to decompose the "integer" s into powers of small primes.
returns eg 108 ==> {{2,2},{3,3}} (ie 2^2*3^3==4*27==108), or 10080 ==> {{2,5},{3,2},{5,1},{7,1}}, or 1196836 ==> {{2,2},{"299209"}}
Each element of the result is a {prime,power} pair, except for the last which may be a lone string:
The default 100th prime is 541, so at that setting this is exact/complete for all inputs <= 541^2 == "292681", and you can easily raise (or lower) that limit, within reason.
However, factors of even a 72-digit number is sufficiently hard that almost all internet security is based on it being a really hard problem. (72 digits is now rather ancient 256-bit security, more recently we have seen 2048-bit security which is ~= factors of 575 decimal digits)
Hence this is designed to "give up" early/in a sensible fashion, eg: mpz_prime_factors(sprintf("%d",power(2*get_prime(101),2)),100) yields {{2,2},{"299209"}}.
Increasing maxprime to 101 obviously makes a big difference, ie: mpz_prime_factors(sprintf("%d",power(2*get_prime(101),2)),101) yields {{2,2},{547,2}}.
Increasing maxprime to say 100,000,000 is likely to make this excruciatingly slow, yet still only fully handle 20-digit numbers.
A length(res[$]) of 1 means it either failed mpz_fits_atom() or is greater than get_prime(maxprime)^2, returned as a (lone) final string.
Also, while all other elements of res are almost certainly phix integer [pairs], res[$][1] may be atom (1GB..4GB), independently of res[$] being a lone string.
Finally, two special cases exist: s=0 yields {}, and s=1 yields {{2,0}}, aka 2^0.
string res =
mpz_factorstring(sequence s) - converts eg {{2,2},{3,3}} to "2^2*3^3", or {{2,5},{3,2},{5,1},{7,1}} to "2^5*3^2*5*7".
s is typically from mpz_prime_factors(), but does not have to be, and s[$] may be {string} (ie unfactored/able).
mpz_re_compose(mpz rop, sequence s) - takes eg s of {{2,2},{3,3}} and sets rop to 108, ie 2^2*3^3 === 4*27 === 108.
The same origins and contents of parameter s as noted with mpz_factorstring() above apply equally here.


mpfr (floating point)

integer precision =
mpfr_get_default_prec(boolean decimal=false) - yields the default floating point precision in bits, initially 53 to match IEEE-754 (approx 15dp).
mpfr_set_default_prec(integer precision) - set the default precision in binary bits or (if -ve) decimal places.
Internally the decimal precision always adds a couple extra to cater for various intermediate values, sometimes you may need a few more.
If decimal is true, mpfr_get_default_prec() returns the minimum number of decimal digits that would be required to distinguish all possible values in the actual binary representation, rather than the number of decimal places that are actually held. For example if we were storing in 8 bits, we could hold -128..127 and need 3 d.p, but could not hold 128..999 [at least, that is, not exactly - remember that it is using extended IEEE-754 internally, and would never go as low as 8 bits]. As above, the +2 on setting means that things should always be sufficient to obtain the desired results/accuracy.
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. v can be integer/atom/string, but not mpfr.
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 below.
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).
Invoking res = mpfr_free(res) when no longer needed is recommended, but occurs automatically if forgotten.
sequence res =
mpfr_inits(integer n, object v=0, integer precision=default_precision, rounding=default_rounding) - initialise n variables to v. (ditto mpfr_free)
Note the method signature differs significantly from the C function of the same name. Can also be coded as eg mpfr {x,y,z} = mpfr_inits(3).
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). (ditto mpfr_free)
Passing an mpfr directly to mpfr_init() is a similar error to using pDword when peek4u(pDword) is required, and likewise there is no practical way to guard against that, you just get meaningless results. In fact mpfr_init_set() can and does test it got an mpfr, it is mpfr_init() that simply cannot perform the inverse useful test that makes sure it has not been passed an mpfr. (Even if it maintained a dirty great list, what’s to say an actual value will not coincide with an allocated memory address?)
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_clear[s]() - use mpfr_free() instead, which releases both mpfr.dll and mpfr.e allocated memory.
integer precision =
mpfr_get_prec(mpfr x, boolean decimal=false) - retrieve the current precision of a specific variable.
See mpfr_get_default_precision() above for more details about the decimal parameter.
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_set_q(mpfr x, mpq q, integer rounding=default_rounding) - set x from an mpq (rational, see below).
mpfr_set_z(mpfr rop, mpz op, integer rounding=default_rounding) - set the mpfr rop from an mpz (integer, see above).
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_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 (2..62) 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("%.1000Rf",x) 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 (just a trivial puts(fn,mpfr_sprintf(fmt,x)) wrapper).
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_floor(mpfr rop, op) - rop := floor(op)
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_si_sub(mpfr rop, integer op1, mpfr op2, integer rounding=default_rounding) - "" except op1 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_si_div(mpfr rop, integer op1, mpfr op2, integer rounding=default_rounding) - "" except op1 is a phix integer, -1GB..+1GB
mpfr_sqr(mpfr rop, op, integer rounding=default_rounding) - rop := op squared with specified rounding
mpfr_sqrt(mpfr rop, op, integer rounding=default_rounding) - rop := sqrt(op) with specified rounding
mpfr_sqrt_ui(mpfr rop, integer op, rounding=default_rounding) - "", except rop is a phix integer 0..1GB
mpfr_pow(mpfr rop, op1, op2, integer rounding=default_rounding) - rop := op1^op2 with specified rounding
mpfr_pow_si(mpfr rop, op1, integer op2, integer rounding=default_rounding) - "" except op2 is a phix integer, -1GB..+1GB
mpfr_ui_pow(mpfr rop, integer op1, mpfr op2, integer rounding=default_rounding) - "" except op1 is a phix integer, 0..+1GB
mpfr_ui_pow_ui(mpfr rop, integer op1, op2, rounding=default_rounding) - "" except op1 and op2 are phix integer, 0..+1GB
Set rop to op1 raised to op2, rounded in the specified direction.
Note that mpfr_pow_si must be used instead of mpfr_pow_ui (no benefit to be had), whereas there are no C mpfr_si_pow[_si] functions.
0^0 is 1 (as per the standard phix power() function), see mpfr.e for more specific details of +/-0/Inf/NaN handling.
mpfr_neg(mpfr rop, op, integer rounding=default_rounding) - rop := -op with specified rounding
Just changes the sign if rop and op are the same variable, otherwise rounding can occur when the precision of rop is less than that of op.
mpfr_sin(mpfr rop, op, integer rounding=default_rounding) - rop := sin(op) with specified rounding
mpfr_log(mpfr rop, op, integer rounding=default_rounding) - rop := log(op) with specified rounding
Set rop to the natural logarithm of op, rounded in the direction rnd.
Set rop to +0 if op is 1 (in all rounding modes), for consistency with the ISO C99 and IEEE 754-2008 standards.
Set rop to -Inf if op is +/-0 (i.e., the sign of the zero has no influence on the result).
mpfr_exp(mpfr rop, op, integer rounding=default_rounding) - rop := exponential(op) with specified rounding
mpfr_gamma(mpfr rop, op, integer rounding=default_rounding) - rop := Gamma(op) with specified rounding
atom res = mpfr_get_si(mpfr op, integer rounding=default_rounding) -- res := op as a machine-word-sized signed integer.
integer res =
mpfr_cmp(mpfr op1, op2) -- Compare op1 and op2. Return +1 if op1>op2, 0 if op1=op2, or -1 if op1<op2.
mpfr_cmp_si(mpfr op1, integer op2) - "" except op2 is a phix integer, -1GB..+1GB
Note that mpfr_cmp_si() should always be used instead of mpfr_cmp_ui(), since the latter would not increase the permitted range.


mpq (rational, aka fractions (mpz/mpz))

mpq res = mpq_init() - Initialise a single mpq, set to 0/1. Invoking res = mpq_free(res) when no longer needed is recommended, but occurs automatically if forgotten.
sequence res = mpq_inits(integer count) - Returns a sequence of count mpq ({} if count=0), set to 0/1, eg mpq {y,z} = mpq_inits(2). (ditto mpq_free)
mpq res = mpq_init_set(mpq op) - Initialise to some prior mpq.
Likewise mpq_init_set_z(), mpq_init_set_si(), and mpq_init_set_str() are as follows, but returning res rather than being passed rop.
mpq_set(mpq rop, op) - Set from another mpq.
mpq_set_z(mpq rop, mpz op) - Set from an mpz (with an implied denominator of 1).
mpq_set_si(mpq rop, integer op1, op2=1) - Set to op1/op2. op1 is -1GB..+1GB, op2 cannot be negative (0..+1GB).
Note that mpq_set_si must be used instead of mpq_set_ui, and the upper limit is therefore 1GB, not 4GB.
mpq_set_str(mpq rop, string str, integer base=0) - Set rop from a string, such as "41" or "41/152".
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.
Note that is done separately for the numerator and denominator, so for instance "0xEF/100" is 239/100, whereas "0xEF/0x100" is 239/256.
object x =
mpq_free(object x) - Clear and deallocate any variables created from mpq_init(), eg x = mpq_free(x), or {y,z} = mpq_free({y,z}).
NB: lhs==rhs intentionally. Should you forget to invoke this routine, mpfr.e will do so automatically.
mpq_clear[s] is not publicly available, since incorrect use can lead to sudden and mysterious program termination (admittedly that can still happen if you break the lhs==rhs rule) - use mpfr_free() instead, which releases both mpir.dll and mpfr.e allocated memory.
mpq_get_num(mpz numerator, mpq rational) - get the numerator of a rational.
mpq_get_den(mpz denominator, mpq rational) - get the denominator of a rational.
mpq_add(mpq rsum, addend1, addend2) - set rsum to addend1 + addend2.
mpq_sub(mpq rdifference, minuend, subtrahend) - set rdifference to minuend - subtrahend.
mpq_mul(mpq rproduct, multiplier, multiplicand) - set rproduct to multiplier * multiplicand.
mpq_div(mpq rquotient, dividend, divisor) - set rquotient to dividend / divisor.
mpq_mul_2exp(mpq rop, op, integer bits) - set rop to op * 2^bits.
mpq_div_2exp(mpq rop, op, integer bits) - set rop to op / 2^bits.
mpq_neg(mpq negated_operand, operand) - set negated_operand to -operand.
mpq_abs(mpq rop, op) - set rop to the absolute value of op.
mpq_inv(mpq inverted_number, number) - set inverted_number to 1/number. If the new denominator is zero, this routine will divide by zero.
mpq_canonicalize(mpq op) - Remove any factors that are common to the numerator and denominator of op, and make the denominator positive.
integer res =
mpq_cmp_si(mpq op1, integer num2, den2) - compare op1 and num2/den2. Return +1|0|-1 for >|=|< respectively.
num2 and den2 are allowed to have common factors. mpq_cmp_si must be used instead of mpq_cmp_ui, den2 may not be negative.
Note the mpfr.e wrapper explicitly converts the C +ve/0/-ve result to +1/0/-1 using sign().

More routines (many more) will be added as needed - there are several hundred (maybe even thousands) of them and, at least in my opinion, it is just simply not worthwhile without something that actually uses and therefore tests them. Each one is usually quite straightforward, but of course this document and Edita/Edix syntax and lookup files also need to be updated.

Note that mpz/mpfr/mpq_clear() are not publically available, use xxx_free() instead, which free both mpf/ir.dll and mpfr.e allocated memory. Since the xxx_clear() have a different calling convention to xxx_free, they are syntax coloured illegal (and linked here). In the early stages, several sudden mysterious crashes were traced to misuse of xxx_clear(), that now occur far less often after migrating to xxx_free()[lhs==rhs].
Also, in contrast, mpr_free_str() and mpz_free_string() are invoked internally for you, so when translating code you can just delete those statements.
Lastly note that all mpf_xxx functions have deliberately been expunged in favour of mpfr_ ones.