mpfr / gmp
Arbitrary Precision Arithmetic: mpir is a fork of gmp (gnu multiple precision), and mpfr is a floatingpoint 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 triplechoc) 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/mpirandmpfr
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 prebuilt 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.0p3 was already installed my Ubuntu box, though I also needed to install mpfrdev, 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 oldfashioned C for ya!
Making (for example) mpfr_clear() private/internal, and wrapping the C mpfr_init() functions and hence effectively making them nonoptional 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 humanreadable 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 hardandfast 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 copyonwrite semantics of standard phix types, eg after
For comparision, here is the (accurate to ˜15/17dp) behaviour of the standard builtin atom:
(
Manually realigned. Aside: 10^60 has seven more digits of accuracy than expected on 32 bit, no idea why, or even how.)
Many operations are performed insitu 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 64bit, 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 anticorruption checks. They all allow NULL.
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.
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.
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 triplechoc) 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/mpirandmpfr
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 prebuilt 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.0p3 was already installed my Ubuntu box, though I also needed to install mpfrdev, 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 oldfashioned C for ya!
Making (for example) mpfr_clear() private/internal, and wrapping the C mpfr_init() functions and hence effectively making them nonoptional 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 humanreadable 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 hardandfast 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 copyonwrite 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^%d1 = %s\n",{e*30,mpz_get_str(n,comma_fill:=true)}) end for n = mpz_free(n)  output:  10^301 = 999,999,999,999,999,999,999,999,999,999  10^601 = 999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999,999  10^901 = 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 e = 1 to 3 do printf(1, "10^%d1 = %,d\n",{e*30,power(10,e*30)1}) end for  output (32bit):  10^301 = 1,000,000,000,000,000,424,684,240,284,426  10^601 = 1,000,000,000,000,000,000,000,082,442,268,486,026,840,024,008,840,604,006,400,824  10^901 = 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 (64bit):  10^301 = 999,999,999,999,999,998,264,846,264,264  10^601 = 999,999,999,999,999,996,186,866,828,000,402,842,426,462,280,408,402,286,040,882  10^901 = 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
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 64bit, 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 anticorruption checks. They all allow NULL.
mpz (integer)
string res = 



mpz res = 



mpz res = 

object x = 





integer count = 





mpz_abs(mpz rop, op)  rop := abs(op)  
mpz_neg(mpz rop, op)  rop := op  
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)} 


boolean res = 

integer res = 

integer res = 



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. 



bool bExact = 

mpz_sqrt(mpz rop, op)  rop := floor(sqrt(op))  






atom res = 





mpz res = 

integer res = 

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 

mpz_tstbit(mpz op, integer bit_index)  Test bit bit_index in op and return 0 or 1 accordingly. 

mpz_scan0(mpz op, integer starting_bit)  Find first 0 in op >= starting_bit. 

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 = 

boolean res = 

random integers

gmp_randinit_mt()  Initialize state for a Mersenne Twister algorithm, and invokes gmp_randseed on it (with a NULL mpz_seed). 


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 n1 inclusive, and store it in rop. 
extras
The following routines are mpfr.especific 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 



boolean res = 
mpz_fits_integer(mpz op)  Return nonzero 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 outofrange. 
boolean res = 
mpz_fits_atom(mpz op, boolean tztrim=false)  Return nonzero 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 outofrange, 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 32bit there is a huge difference between 31bit integers and 53bit atoms, on 64bit 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 = 

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. 


string res = 


mpfr (floating point)


integer rounding = 

mpfr res = 



mpfr res = 

object x = 

integer precision = 





mpfr_free_str()  invoked automatically within mpfr_get_str() and hence not publically available (there is no need for it to be).  
sequence res = 

string res = 













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_neg(mpfr rop, op, integer rounding=default_rounding)  rop := op with specified rounding 

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 +0 if op is 1 (in all rounding modes), for consistency with the ISO C99 and IEEE 7542008 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 machinewordsized signed integer. 
integer res = 

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

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. 

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_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. 


object x = 

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 = 

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.