mpfr / gmp
Arbitrary Precision Arithmetic: mpfr is a floating-point extension of gmp (gnu multiple precision).
It has been wrapped/is available for many programming languages, including Ada, C, C++, Fortran, Haskell, Java, Julia, Lisp, .NET, 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), with the required windows dlls in builtins\(32|64), but you may have to install gmp/mpfr manually on Linux.
There is also a hand-crafted drop-in replacement in Phix/pwa/builtins/mpfr.js which is handled automatically for you by pwa/p2js, and obviously that does not need any dlls.
Note that mpfr are effectively implemented as (limited precision) mpq and transcendental functions such as mpfr_sin(), mpfr_log(), mpfr_exp(), and mpfr_gamma() are not supported (though, perhaps suprisingly, mpfr_sqrt() and mpfr_const_pi() are already implemented). Also note that mpfr.js (effectively) assumes precision means the number of decimal places after the decimal point. The mechanisms for precision and rounding are quite different in the gmp dlls and mpfr.js, in other words I have probably recreated/interpreted them imperfectly or even quite wrongly, so you should expect a few little glitches, especially in the first few releases. All I know is that right now it gets several hundred digits right on several dozen examples, and in fact I have no more examples left to test. A couple of debug sessions have also suggested that the "discard precision" loop may be ridiculously expensive, but I have not figured out any way to separate the time of that from the time spent elsewhere, so it may just be more cheap iterations, and 50% or even 70% is perfectly fine whereas 95% or more is probably not, and 99% definately not!
Version: mpfr: 4.1.0, gmp:5.1.2, two dlls, 3.5MB, 1.5MB zipped (same/twice that for 32 and 64-bit)
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. Note that eg mpfr_gamma_inc() requires 4.0.0+, I think.
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 that where z is an mpz/mpfr/mpq variable, on desktop/Phix atom(z) will return true (and integer(z) sometimes true) whereas on pwa/p2js atom(z) will always return false. Otherwise the types are "strong" since they are properly tagged on both systems, in other words the usual unsafe warning does not apply to explicit invocations of mpz(), mpfr(), or mpq().
Also mpz/mpfr/mpq declare reference variables, as opposed to the copy-on-write semantics of standard phix types, eg after
Output:
Output:
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, and mpz 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.
Some were largely written for simple convenience, as opposed to any improvements in 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. Note that mpz_set_v() as documented above more strictly belongs in this section.
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. Likewise there may be a few still missing from mpfr.js but so far any such have proved remarkably easy to add, with the trivial updates to js.syn and p2js_keywords.e being as fiddly if not moreso than the new code itself.
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 (the latter can I believe start small and grow whereas the former are always allocated their maximum size) and likewise the low-level mpn_xxx functions (which require the exact number of "limbs" to be specified on every call and as started in the gmp docs have a deliberately "non-coherent calling interface") are not wrapped either.
It has been wrapped/is available for many programming languages, including Ada, C, C++, Fortran, Haskell, Java, Julia, Lisp, .NET, 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), with the required windows dlls in builtins\(32|64), but you may have to install gmp/mpfr manually on Linux.
There is also a hand-crafted drop-in replacement in Phix/pwa/builtins/mpfr.js which is handled automatically for you by pwa/p2js, and obviously that does not need any dlls.
Note that mpfr are effectively implemented as (limited precision) mpq and transcendental functions such as mpfr_sin(), mpfr_log(), mpfr_exp(), and mpfr_gamma() are not supported (though, perhaps suprisingly, mpfr_sqrt() and mpfr_const_pi() are already implemented). Also note that mpfr.js (effectively) assumes precision means the number of decimal places after the decimal point. The mechanisms for precision and rounding are quite different in the gmp dlls and mpfr.js, in other words I have probably recreated/interpreted them imperfectly or even quite wrongly, so you should expect a few little glitches, especially in the first few releases. All I know is that right now it gets several hundred digits right on several dozen examples, and in fact I have no more examples left to test. A couple of debug sessions have also suggested that the "discard precision" loop may be ridiculously expensive, but I have not figured out any way to separate the time of that from the time spent elsewhere, so it may just be more cheap iterations, and 50% or even 70% is perfectly fine whereas 95% or more is probably not, and 99% definately not!
Version: mpfr: 4.1.0, gmp:5.1.2, two dlls, 3.5MB, 1.5MB zipped (same/twice that for 32 and 64-bit)
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. Note that eg mpfr_gamma_inc() requires 4.0.0+, I think.
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 that where z is an mpz/mpfr/mpq variable, on desktop/Phix atom(z) will return true (and integer(z) sometimes true) whereas on pwa/p2js atom(z) will always return false. Otherwise the types are "strong" since they are properly tagged on both systems, in other words the usual unsafe warning does not apply to explicit invocations of mpz(), mpfr(), or mpq().
Also mpz/mpfr/mpq 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. Note that you should never
pass an mpfr variable to deep_copy(): doing so would almost guarantee being left with a
reference to reclaimed memory, hence the latter now [1.0.5+] has a guard in place to prevent such attempts.
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)
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,999For 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
-- 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 -- 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,466Many 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, and mpz 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.
general
string res = |
|
|
|
mpz (integer)
mpz res = |
|
|
|
mpz res = |
|
object x = |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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 = |
|
boolean res = |
|
boolean res = |
|
integer res = |
|
integer res = |
|
|
|
|
|
|
|
bool bExact = | |
|
|
mpz_sqrt(mpz rop, op) - rop := floor(sqrt(op)) | |
|
|
|
|
integer res = |
|
|
|
|
|
|
|
bool res = | |
|
|
integer res = |
|
|
|
integer res = |
|
integer res = |
mpz_cmp_si(mpz op1, integer op2) - "", except op2 is a phix integer, -1GB..+1GB - use this instead of mpz_cmp_ui(). |
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 if op1 is odd. Compatibility shim for the C macro mpz_odd_p |
boolean res = | mpz_even(mpz op1) - returns true if 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. |
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. Note there is no mpz_size() function in mpfr.e ("" in limbs), since it is not very useful and would not be p2js compatible. |
string res = |
|
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, as opposed to any improvements in 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. Note that mpz_set_v() as documented above more strictly belongs in this section.
string res = |
|
mpz_add_si(mpz rop, op1, integer op2) - invokes mpz_add_ui() or mpz_sub_ui() with -op2, crashes if op2 is -(maxint+1). mpz_sub_si(mpz rop, op1, integer op2) - invokes mpz_sub_ui() or mpz_add_ui() with -op2, crashes if op2 is -(maxint+1). mpz_si_sub(mpz rop, integer op1, mpz op2) - 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_addmul_si(mpz rop, op1, integer op2) - as mpz_addmul_ui() except op2 is +/-1GB mpz_submul_si(mpz rop, op1, integer op2) - as mpz_submul_ui() except op2 is +/-1GB mpfr_addmul_si(mpfr rop, op1, integer op2, rounding=default_rounding) - to match mpz_addmul_si() (note there is no mpfr_addmul_ui) |
|
|
|
boolean res = |
|
boolean res = |
|
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. |
|
|
sequence res = |
|
sequence res = |
|
string res = |
|
boolean res = |
|
|
|
|
|
|
mpfr (floating point)
|
|
|
|
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). | |
string res = |
|
string 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_abs(mpfr rop, op, integer rounding=default_rounding) - rop := abs(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 754-2008 standards. Set rop to -Inf if op is +/-0 (i.e., the sign of the zero has no influence on the result). Aside: unlike log/ln and friends, there is no mpfr_ln(), since that would not be a valid gmp/mpfr name. Not supported under pwa/p2js. |
|
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 mpfr_gamma_inc(mpfr rop, op, integer rounding=default_rounding) - incomplete gama function on op and op2 mpfr_zeta(mpfr rop, op, integer rounding=default_rounding) - rop := Riemann Zeta function on op, with specified rounding mpfr_zeta_ui(mpfr rop, integer op, rounding=default_rounding) - same, except op is a Phix integer, 0..1GB |
|
|
mpfr_get_d_2exp(mpfr op, integer rounding=default_rounding)
(where op is rounded to double precision using the given rounding mode). If op is zero, then a zero of the same sign is returned, and e is set to 0. If op is NaN or an infinity, then the corresponding d is returned, and e is undefined. Not supported under pwa/p2js. |
atom res = |
mpfr_get_si(mpfr op, integer rounding=default_rounding) -- res := op as a machine-word-sized signed integer. mpfr_get_d(mpfr op, integer rounding=default_rounding) -- res := op as a double. As usual, use mpfr_get_si() instead of mpfr_get_ui(), or possibly mpfr_get_d(). |
integer res = |
|
mpq (rational, aka fractions (mpz numerator/mpz denominator))
|
|
|
|
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 tgt. |
object x = |
|
|
|
string res = |
|
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. |
|
atom res = |
mpq_get_d(mpq op) - get the value as a phix atom, rounded to ~15|19 significant digits, Inf/NaN/0 if too big/small. |
|
|
mpq_add(mpq rsum, addend1, addend2) - set rsum to addend1 + addend2. mpq_add_si(mpq rsum, addend1, integer n, d=1) - set rsum to addend1 + n/d. 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 rop, op) - set rop to -op. mpq_abs(mpq rop, op) - set rop to the absolute value of op. mpq_inv(mpq rop, op) - set rop to 1/op. If the new denominator is zero, this routine will divide by zero. |
|
integer res = |
|
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. Likewise there may be a few still missing from mpfr.js but so far any such have proved remarkably easy to add, with the trivial updates to js.syn and p2js_keywords.e being as fiddly if not moreso than the new code itself.
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 (the latter can I believe start small and grow whereas the former are always allocated their maximum size) and likewise the low-level mpn_xxx functions (which require the exact number of "limbs" to be specified on every call and as started in the gmp docs have a deliberately "non-coherent calling interface") are not wrapped either.