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
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)
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
-- 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, 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 =
|
-
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.
The pwa/p2js version simply always returns "" (and ignores any parameters).
|
|
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".
The pwa/p2js version (effectively) returns {"JavaScript","mpfr.js","n/a"}, ignoring bAsNumSeq,
and is therefore acceptable for display but not version checking (obviously an "if platform()!=JS then" is advised).
|
mpz (integer)
|
mpz res =
|
-
mpz_init(object v=0, integer bitcount=0) - Initialise. v can be integer/atom/string, but not mpz, for that use mpz_init_set().
-
As of 1.0.0 you can use exponentiation, such as "1e200" or even "2.5e1" (but not "2.5" or "2.55e1" aka non-integer,
or bases>10). 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 (equivalent of the) latter three can take a bitcount in mpfr.e (albeit completely ignored in mpfr.js).
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 500 dwords). A non-zero bitcount makes it possible to avoid reallocations if a suitable size is known in advance.
You c/should use (eg) mpz_sizeinbase(x,2) in some prior test run to determine an initial bitcount (and check it actually makes the
thing 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, (but not mpz), 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.
(The problem being that an actual value could 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 (not actually enforced but you should pretend it is). 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 - use this instead of mpz_set_ui().
mpz_set_d(mpz rop, atom op) -- "" except op is a phix atom.
mpz_set_q(mpz rop, mpq op) -- "" except op is an mpq (truncated to make it an integer).
mpz_set_str(mpz rop, string s, integer base=0) -- set x from s in the specified base. See also mpz_init.
A base of 0 allows bases 2 and 16 to be auto-detected (after a leading +/-) by '0b', '0B', '0x' or '0X', otherwise base 10 is assumed.
A Phix-specific enhancement is to allow eg "1e200" (only base 10 is officially supported, though some others may work).
General support for bases other than 2/8/10/16 (2..62 matching desktop/Phix) added to pwa/p2js in 1.0.2.
Note that desktop/Phix quietly ignores whitespace in s, whereas pwa/p2js (probably) won’t, largely due to js/BigInt details
beyond my control, and/or not really being worthwhile implementing fiddly cleanup/workarounds, at least so far not yet.
mpz_set_v(mpz rop, object v) -- a Phix-specific shorthand for the above three (with said decimal assumption)
|
|
|
-
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.
Not supported under pwa/p2js.
|
|
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)
Not supported under pwa/p2js.
|
|
|
-
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_addmul(mpz rop, op1, op2) - rop := rop + op1 * op2
-
mpz_addmul_ui(mpz rop, op1, integer op2) - "" except op2 is a phix integer, 0..1GB.
mpz_addmul_si(mpz rop, op1, integer op2) - (an extra) "", +/-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_submul(mpz rop, op1, op2) - rop := rop - op1 * op2
-
mpz_submul_ui(mpz rop, op1, integer op2) - "" except op2 is a phix integer, 0..1GB.
mpz_submul_si(mpz rop, op1, integer op2) - (an extra) "", +/-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 - use this instead of mpz_mul_ui()
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.
|
|
|
-
mpz_divexact(mpz r, n, d) - r := n/d
-
mpz_divexact_ui(mpz r, n, integer d) - "" except d is a phix integer, 1..1GB
q = n/d, but only when it is known in advance that d divides n exactly.
The mpz_divexact[_ui] routines are much faster than the other division functions, and are the best choice
when exact division is known to occur, for example reducing a rational to lowest terms.
|
|
boolean res =
|
-
mpz_divisible_p(mpz n, d) - returns non-zero if n is exactly divisible by d.
|
|
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 2b
-
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 that rule (only) 0 is considered
divisible by 0.
|
|
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 d 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_nthroot(mpz rop, op, integer n) - rop:=trunc(power(op,1/n)).
-
Set rop to the truncated integer part of the nth root of op.
|
|
|
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_and(mpz rop, op1, op2) - rop := and_bits(op1,op2)
|
|
|
-
mpz_xor(mpz rop, op1, op2) - rop := xor_bits(op1,op2)
|
|
integer res =
|
-
mpz_popcount(mpz op) - res := count_bits(op)
-
If op >= 0, return the population count of op, which is the number of 1 bits in the binary
representation. If op < 0, the number of 1s is infinite, and the return value is ULONG MAX,
the largest possible mp_bitcnt_t.
|
|
|
-
mpz_fac_ui(mpz rop, integer n) - Set rop to the factorial of n.
-
Equivalent, for small n, and bar the mpz result, to factorial(n).
|
|
|
-
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_bin_uiui(mpz rop, integer n, k) - calculate n choose k, ie n!/((n-k)!*k!)
-
Equivalent, for small n and k, and bar the mpz parameter, to choose(n,k)
(replaces a non-standard mpz_binom() function in the previous release)
|
|
bool res =
|
-
mpz_invert(mpz rop, op1, op2) - Set rop to the inverse of op1 modulo op2.
-
Compute the inverse of op1 modulo op2 and put the result in rop.
If the inverse exists, the return value is true and rop will satisfy 0 <= rop < op2.
If an inverse doesn’t exist the return value is false and rop is undefined.
|
|
|
-
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.
|
|
integer 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, though when translating code you may hit problems
with op2 being limited (in this Phix wrapper) to 0..+1GB rather than the C api of mpir_ui being 0..+4GB.
It would of course not be difficult to add a phix-specific mpz_gcd_d() with an atom parameter to ease that.
Also note this wrapper implodes if op2<0, whereas the C api would quietly cast it to an unsigned integer, and
likewise I forsee no real difficultly performing a proper 2’s-compliment abs(), bar perhaps minint.
|
|
|
-
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.
|
|
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 - 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
|
|
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.
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 =
|
-
mpz_get_str(mpz op, integer base=10, boolean comma_fill=false)
- Return op as a string in the specified base (2..62).
-
Note that mpz_free_string() is taken care of automatically for you by mpfr.e.
See also mpz_get_short_str().
|
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_get_short_str(mpz op, integer ml=20, base=10, boolean comma_fill=false, string what="digits")
-
Equivalent to shorten(mpz_get_str(op,base,comma_fill),ml:=ml))
but much faster, since it does not construct potentially hundreds of thousands of middle digits and just throw them away.
Sadly no real impact under pwa/p2js, oh well.
Note that ml is expressed in characters, so if you specify comma_fill and a length of 4 then you’ll
get (say) "12,3...,678 (nn digits)" back, the commas where they w/should rightly be on a full-length result,
and obviously the trailing part will probably look nicer when ml is 4k-1, ie 3,7,11,15 etc.
|
|
|
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)
|
|
mpz/int res =
|
-
mpz_min(sequence s, boolean return_index=false) - return the smallest element of s, or its 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 just-out-of-range, since that is the first value that "accidentally" fits, by ending in
a binary 0, that is when tztrim is false, otherwise it returns true for that and
9,007,199,254,740,994 but false for 9,007,199,254,740,993, in other words tztrim=true
makes it start returning true for every 2nd, 4th, 8th, etc integer value.
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.
While mpz_fits_integer() is somewhat of a compatibility fudge, in that under pwa/p2js it artificially recreates a boundary that does not
technically exist, mpz_fits_atom() is in contrast an absolutely perfect fit: you could not design a more appropriate routine for JavaScript
even if you tried.
|
|
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 effectively one bit smaller than a C signed long,
and two bits smaller than a C unsigned 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.
Note that while prime_factors() accepts a -1 in maxprime there is no such facility here.
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.
|
|
sequence res =
|
-
mpz_pollard_rho(mpz_or_string s, bool bAsStrings=false)
- find factors of s, usually significantly faster than mpz_prime_factors().
-
Uses a pollard rho method for s greater than 100,000,000 else mpz_prime_factors() ie primes<=10,000.
The result is a slightly flattened version of that from mpz_prime_factors(), eg {{2,2},"837568909977835437"}, or all strings if asked for,
in other words mpz_pollard_rho("8",true) yeilds {"2","2","2"} whereas mpz_pollard_rho("8") yields {{2,3}}.
Any factors which fail mpz_fits_atom() will certainly be returned as strings, however some smaller (but > 100,000) may also slip through,
In other words you might get "456773" instead of 456773, but I’m sure you’ll cope, and unlike mpz_prime_factors() there may be
multiple strings in the result, even when bAsStrings is false.
Since this routine invokes mpz_prime() in the main loop anyway, there is not likely to be much benefit from
calling it first.
Obviously mpz_prime_factors() remains a better choice should you not be interested in factors over 10,000 and maybe 100,000.
|
|
sequence res =
|
-
mpz_factors(object s, object include1=0, bool bAsAtmStr=true,
bSort=true)
-
Constructs the equivalent of factors(n,include1) from the prime powers.
s can be an integer atom (within the native precision limits), a string, a dword-sequence result from mpz_pollard_rho, or a "tagged" mpz
which looks like this: {"mpz",z}. Take special care to wrap an mpz like that; there is a high risk of treating the raw memory address of
an unwrapped mpz as an integer (and pretend to work, as opposed to terminating with a more helpful error).
if bAsAtmStr is false, the result is a sequence of mpz, otherwise (the default) result is a mixture of atoms
(when they fit) and strings (when they exceed the 53/64 bit precision limit of normal Phix atoms).
A bAsAtmStr parameter of 2 means "return all as strings".
Just as mpz_pollard_rho() can return strings for prime factors that do not fit in a phix atom, any strings in a dword_sequence s are assumed
to be (such) primes, and the result may be a mixture of integers and strings.
Sorting probably adds very little overhead, but you never know...
?mpz_factors(12,bSort:=false) ==> {1,2,4,3,6,12}
?mpz_factors(12,bSort:=true) ==> {1,2,3,4,6,12}
?mpz_factors("12345678999999977",bAsAtmStr:=true) ==> {1,35604581,346744117,"12345678999999977"}
?mpz_factors("12345678999999977",bAsAtmStr:=2) ==> {"1","35604581","346744117","12345678999999977"}
Note this is not expected to handle duplicate input strings correctly: should you have (somehow) managed to factor a 32-digit+ number, good
on you, but you have long since passed the bounds of what is realistically feasible, unless, that is, all the prime factors are reasonably
small integers.
|
|
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.
|
|
boolean res =
|
-
mpz_prime(mpz p, integer prob=5) - Determine whether p is a probable prime with the chance of error being at most 1 in 2^prob.
-
The return value is true (1) if p is probably prime, or false (0) if p is definitely composite.
Replaces mpz_probable_prime_p() which has been deprecated in version 1.0.0, as have the
randstate functions gmp_randinit_mt(), gmp_randseed(),
and gmp_randclear().
Additionally there was a js-compatible mpz_prime_mr() for a while, with desktop/Phix still using the
official (mpfr.dll) version by default, but the performance of the latter fell off a cliff (quite seriously) with a newer mpfr.dll
(mid-2022) and hence desktop/Phix now uses the same code as pwa/p2js, which is
deterministic to 3,317,044,064,679,887,385,961,981 whereas the gmp docs do(/did) not appear to make any such guarantee.
By deterministic I mean that all numbers up to that limit yield perfectly correct results, and there is no point or gain/loss
specifying a prob other than 1 for any p below or equal to that limit.
|
|
|
-
mpz_nextprime(mpz rop, op) - set rop to the next prime greater than op.
-
I needed this for a rosettacode task but found it horribly slow and don’t recommend it at all.
Note that it no longer tries to use either of the (deprecated) mpz_nextprime or (its presumed successor)
mpz_next_prime_candidate C functions and instead just throws successive odd numbers at mpz_prime(), so it may be
even slower, but in my tests at least, not by much.
|
|
|
-
mpz_rand(mpz n, range) - set n to a random number 0..range-1
|
|
|
-
mpz_rand_ui(mpz n, integer range) - set n to a random number 0..range-1
-
mpz_rand(n,r) is equivalent to mpz_set_si(n,rand(mpz_get_integer(r))-1), and
mpz_rand_ui(n,i) is actually implemented as mpz_set_si(n,i); mpz_rand(n,n).
Obviously you should just use rand() for small numbers
(<= machine-word) you don’t actually need to be stored in an mpz.
Note that mpz_rand[_ui](n,n) yields 0..n-1 whereas rand(n) yields 1..n.
Also, rand() is maybe 64 times faster than mpz_rand(), and neither should
be considered cryptographically secure. In the case of rand() that would
be down to a trivial algorithm, in the case of mpz_rand() that would be
down to non-linear timings, with both possibly containing other faults.
Lastly, mpz_urandomm() has been deprecated, due to
pwa/p2js difficulties, and
mpz_rand[_ui]() should now be used instead (admittedly from a gmp perspective
that means reverting to a single shared randstate, and conversely p2js.js relies
on Math.random, the lack of a randstate being the difficultly mentioned above).
|
mpfr (floating point)
|
integer precision =
|
-
mpfr_get_default_precision(boolean decimal=false)
- yields the default floating point precision in bits or decimal places.
mpfr_set_default_precision(integer precision) - set the default precision in binary bits or (if -ve) decimal places.
-
Initially the precision is 53 to match IEEE-754 (approx 15dp).
It always adds four to a decimal precision to cover intermediate values, sometimes you may need a few more.
If decimal is true, mpfr_get_default_precision() 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 +4 on setting helps ensure it is sufficient to obtain the desired results/accuracy, but no such fixed safety margin
could ever offer any guarantee for all algorithms/intermediates needed.
mpfr_get_default_prec() has been renamed mpfr_get_default_precision() in 1.0.0
mpfr_set_default_prec() has been renamed mpfr_set_default_precision() in 1.0.0
Note that pwa/p2js emulates floating points using BigInt rationals, and as such there is
no real limited precision handling at all, as yet (intuitively I imagine/guess there is a point that I haven’t figured out where
you can just lop digits off (/2 or /10) both numerator and divisor without losing too much accuracy, or perhaps better yet hold/calculate
rationals with a logical 10(or 2)^k where k can go quite -ve or +ve and [adjust k and] just lop the numerator when it gets too big).
|
|
integer rounding =
|
-
mpfr_get_default_rounding_mode() - yields MPFR_RNDx (0..4, initially N==0), as below.
mpfr_set_default_rounding_mode(integer rounding) - sets the rounding mode.
-
MPFR_RND(N/Z/U/D/A) stand for Nearest, toward Zero/+inf(Up)/-inf(Down), and Away from zero.
No effort has yet been made to implement rounding modes in pwa/p2js and they are
currently (almost?) completely ignored.
|
|
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_ui, 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,
and of course that recommendation/handling is exactly the same for the following two functions as well.
|
|
sequence res =
|
-
mpfr_inits(integer n, object v=0, integer precision=default_precision, rounding=default_rounding) - initialise n variables to v.
-
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).
v may be an integer, atom, string, (but not mpfr), or a sequence of length n of said.
|
|
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 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 (realistically) perform
the inverse useful test that ensures it has not accidentally been passed an mpfr.
(Since an actual value could coincide with an allocated memory address)
The functions mpfr_init_set_q() and mpfr_init_set_z() behave as mpfr_set_q and mpfr_set_z, but
returning res rather than being passed tgt.
|
|
|
-
mpfr_set(mpfr tgt, src, integer rounding=default_rounding) - set tgt from an existing mpfr variable.
-
mpfr_set_si(mpfr tgt, integer i, rounding=default_rounding) - set tgt from a (machine-word-sized) integer.
In mpfr.e, i is (deliberately) limited to +/-1GB, use this instead of mpfr_set_ui, or
mpfr_set_d to cope with 1..4GB.
mpfr_set_d(mpfr tgt, atom a, integer rounding=default_rounding) - set tgt from a normal phix atom.
mpfr_set_str(mpfr tgt, string s, integer base=0, rounding=default_rounding) - set tgt from a string.
A base of 0 allows bases 2 and 16 to 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 (see also p/P and base 0 caveat below).
Note that desktop/Phix uses the gmp function, which has accumulated hacks and tweaks from the past several decades, such as using 'p'
and 'P' for the exponent prefix in bases 2 and 16, or "-@nAn@(This_Is_Not_17)" being a "valid" representation of NaN in base 17,
whereas pwa/p2js uses a hand-crafted replacement, so some inconsistencies are inevitable,
though by and large any normal/sane use should be and is of course perfectly consistent.
You should (I think) replace p/P with @ and always use a '.' for the decimal point - gmp also accepts the current locale but pwa/p2js does not (knowingly).
Non-base 10 exponents are not offically supported/adequately tested (under pwa/p2js), and likewise the above-mentioned auto-detection prevents
proper handling of '@' (ie it looks for 'e'/'E' when base<=10, and '@' only when base>10). Also, some browsers may not honour bases 37..62.
A Phix-specific enhancement is to allow eg "1/10", which arose because pwa/p2js implements mpfr in an mpq-like fashion, hence "1/10" is spot-on
whereas "0.1" is or w/could be inaccurate, and that has now (as of 1.0.2) been supported on desktop/Phix as well
(despite the fact that said inaccuracy persists within the IEEE-754-based gmp library on the latter platform).
Fairly obviously such differences also "bleed through" to the likes of mpfr_init(), etc. I will of course happily fix any further discrepancies
that may arise between desktop/Phix and pwa/p2js, as long as they seem to be actually worthwhile fixing and/or are indeed even possible.
mpfr_set_q(mpfr tgt, mpq q, integer rounding=default_rounding) - set tgt from an mpq (rational, see below).
mpfr_set_z(mpfr tgt, mpz z, integer rounding=default_rounding) - set tgt from an mpz (integer, see above).
|
|
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 (not actually enforced but you should pretend it is).
Should you forget to invoke this routine, mpfr.e will do so automatically.
Neither mpfr_clear() nor mpfr_clears() are publicly available,
use mpfr_free() instead, which releases both mpfr.dll and mpfr.e allocated memory (and if you donb’t
understand the importance or significance of that, don’t ask. No, really, just don’t).
|
|
integer precision =
|
-
mpfr_get_precision(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_precision(mpfr x, integer precision) - precision is the number of bits required for the mantissa, see mpfr_set_default_precision() above.
mpfr_get_prec() has been renamed mpfr_get_precision() in 1.0.0
mpfr_set_prec() has been renamed mpfr_set_precision() in 1.0.0
Note that on desktop/Phix, mpfr_set_precision() sets x to NaN, whereas pwa/p2js sets it to 0, not likely to be an issue since it
would normally be set immediately afterwards anyway. As noted above, there isn’t much of anything in the way of precision
handling within pwa/p2js as yet.
|
|
|
mpfr_free_str() - invoked automatically within mpfr_get_str() and hence not publically available (there is no need for it to be).
|
|
string res =
|
-
mpfr_get_str(mpfr x, integer base=10, n=0, rounding=default_rounding, bool comma_fill=false)
- print x in base (2..62) to precision n(0=full).
-
Prior to version 1.0.0 this routine returned {string digits, integer exponent} and digits did not contain a decimal point.
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.
Note this routine is not supported under pwa/p2js, use mpfr_get_fixed() [see next] instead.
(The main reason being that I’m not really getting the n parameter, and certainly not capabable of accurately reproducing that behaviour.)
|
|
string res =
|
-
mpfr_get_fixed(mpfr x, integer dp=6, base=10, bool comma_fill=false, integer maxlen=0)
- get a string representation of x
-
Returns a string representation of x in base to dp digits after the "decimal" point.
A dp of 0 prints up to 16 "decimal" places, obviously if you want "zero" dp then you’ll have to mpfr_floor(x,x) or similar first.
A negative dp is allowed, which causes it to print those as trailing zeroes, so eg 11111.111 with -3 yeilds "11000".
A non-zero maxlen (introduced in 1.0.2) makes it resort to scientific notation when necessary. Once you make that flip, there are never any commas
in the result, however due to "bouncy length syndrome" there may be some cases where maxlen and comma_fill at the same time run interference with
each other, eg maxlen 8, "1,234,567" -> "1.2345e6", but internally that works via "1234567?", padding a comma-less result with '?' characters
before patching in the dot and trailing "e6", and in some cases some of those '?' may leak through. Quite probably a future release will simply
replace the '?' padding with '0' and assume that will always be ok, unless of course some counter-examples are found/reported first.
As of 1.0.4 the last digit is rounded in the same fashion as bankers_rounding().
|
|
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.
Note this routine is not supported under pwa/p2js, use
mpfr_get_fixed() [see above] instead.
|
|
|
-
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.
Obviously this routine is not supported under pwa/p2js, since you cannot
perform file/io from within a browser (although you might be able to ask a server to do that for you).
|
|
|
-
mpfr_floor(mpfr rop, op) - rop := floor(op)
|
|
|
-
mpfr_ceil(mpfr rop, op) - rop := ceil(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, rounding=default_rounding) - "" except op2 is a phix integer, -1GB..+1GB
mpfr_add_d(mpfr rop, op1, atom op2, rounding=default_rounding) - "" except op2 is a phix atom.
mpfr_addmul_si(mpfr rop, op1, integer op2, rounding=default_rounding) - (an extra) rop += op1*op2
|
|
|
-
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_sub_d(mpfr rop, op1, atom op2, integer rounding=default_rounding) - "" except op2 is a phix atom.
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_mul_d(mpfr rop, op1, atom op2, integer rounding=default_rounding) - "" except op2 is a phix atom.
mpfr_mul_z(mpfr rop1, op1, mpz op2, integer rounding=default_rounding) - "" except op2 is an mpz
|
|
|
-
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 ""
mpfr_div_d(mpfr rop, op1, atom op2, integer rounding=default_rounding) - "" except op2 is a phix atom.
mpfr_div_z(mpfr rop, op1, mpz op2, integer rounding=default_rounding) - "" except rop := mpfr/mpz
As usual, use mpfr_div_si() instead of mpfr_div_ui(), or possibly mpfr_div_d().
|
|
|
mpfr_fmod(mpfr r, x, y, integer rounding=default_rounding) - r := mod(x,y)
Set r to the value of x - ny, rounded as specified, where n is the integer quotient of x divided by y, rounded toward zero.
|
|
|
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
mpfr_si_pow_si(mpfr rop, integer op1, op2, rounding=default_rounding) - "" except op1 and op2 can be -1GB..+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, that last one is a Phix-specific extra.
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
mpfr_abs(mpfr rop, op, integer rounding=default_rounding) - rop := abs(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_const_pi(mpfr x, integer rounding=default_rounding) - x := PI
-
As precise as needed/specified in some prior mpfr_set_precision(x) or the x = mpfr_init() call.
Unfortunately trying to print more than about 42,000 digits of pi crashes (in the dll) on desktop/Phix.
Note that in pwa/p2js, mpfr_const_pi() uses a very simple
BPP based method, which yields 12,000 digits in 4.5s, without
any attempt at optimisation, however it should probably not be used for more than say 1000 decimal digits
(and in fact I have limited use on rosettacode to 1000 digits), plus it makes no attempt to cache results.
In contrast, the gmp routine invoked by desktop/Phix uses a state-of-the-art
Chudnovsky method - should anyone want to contribute a (normal
hll Phix, or even JavaScript) version of that, feel free. Note however that generating a 100,000 digit string
in a browser is going to cause it significantly more grief than a 100,000,000 digit string does desktop/Phix,
so there are other already fairly close limits to how useful any such effort could be, somewhere between 4 and
8 times as many digits at best I would guess.
|
|
|
-
mpfr_const_euler(mpfr x, integer rounding=default_rounding) - x := gamma (approx 0.57721566490153286)
-
Set x to Eulers constant aka C aka gamma, not to be confused with Eulers number aka EULER, 2.7182818284590452.
As precise as needed/specified in some prior mpfr_set_precision(x) or the x = mpfr_init() call.
Similar limits to mpfr_const_pi of printing excessive digits likewise apply here.
Not supported under pwa/p2js.
|
|
|
mpfr_sin(mpfr rop, op, integer rounding=default_rounding) - rop := sin(op) with specified rounding
Not supported under pwa/p2js.
|
|
|
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 specified direction.
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
Not supported under pwa/p2js.
|
|
|
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
Not supported under pwa/p2js.
|
|
{atom d, integer e} =
|
mpfr_get_d_2exp(mpfr op, integer rounding=default_rounding)
Return d and e such that 0.5<=abs(d)<1 and d*2e equals op
(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 =
|
-
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
As usual, 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 numerator/mpz denominator))
|
mpq res =
|
-
mpq_init(object v=0) - Initialise a single mpq, by default set to 0/1.
-
v may be an integer or string, but not a non-integer atom or mpq.
Note that mpq_init_set_si() can set a non-1 denominator in a way this routine (when passed an integer) simply cannot.
Invoking res = mpq_free(res) when no longer needed is recommended, but occurs automatically if forgotten,
and obviously that recommendation/handling applies aidentically to half a dozen or so of the following routines.
|
|
sequence res =
|
-
mpq_inits(integer n, object v) - Returns a sequence of n mpq ({} if n=0), set to 0/1, eg mpq {y,z} = mpq_inits(2). (ditto mpq_free)
-
v may be an integer or string, or a sequence of length n of said (ditto no non-integer atom(s) or mpq elements).
|
|
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 =
|
-
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 (not actually enforced but you should pretend it is).
Should you forget to invoke this routine, mpfr.e will do so automatically. Neither mpq_clear() nor
mpq_clears() are 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 mpq_free() instead, which releases both mpir.dll and mpfr.e allocated memory.
|
|
|
-
mpq_set(mpq tgt, src) - Set tgt from another mpq.
mpq_set_z(mpq tgt, mpz n, d=NULL) - Set tgt from an mpz (d=NULL results in an implied denominator of 1).
mpq_set_si(mpq tgt, integer n, d=1) - Set to n/d. n is -1GB..+1GB, d cannot be zero or negative (1..+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 tgt, string s, integer base=0) - Set tgt 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.
That is done separately for the numerator and denominator, so for instance "0xEF/100" is 239/100, whereas "0xEF/0x100" is 239/256.
|
|
string res =
|
-
mpq_get_str(mpq op, integer base=10, boolean comma_fill=false)
- Return op as a string in the specified base (2..62).
-
The result will be of the form "num/den", or if the denominator is 1 then just "num".
|
|
|
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_canonicalize(mpq op) - Remove any factors common to numerator and denominator, and make the latter positive.
-
This is performed automatically by mpq_set_si/str and mpq_set_z with a non-null d, and by most of the following routines.
|
|
|
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 =
|
-
mpq_cmp(mpq op1, op2) - compare op1 and op2. Return +1|0|-1 for >|=|< respectively.
-
Note the mpfr.e wrapper explicitly converts the C +ve/0/-ve result to +1/0/-1 using sign().
|
|
integer res =
|
-
mpq_cmp_si(mpq op1, integer n, d=1) - compare op1 and n/d. Return +1|0|-1 for >|=|< respectively.
-
n and d are allowed to have common factors. mpq_cmp_si must be used instead of mpq_cmp_ui,
d may not be zero or 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. 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.
