Category talk:Wren-gmp: Difference between revisions
m
→Source code (Wren): Now uses Wren S/H lexer.
(Added preamble and source code for new Wren-gmp module.) |
m (→Source code (Wren): Now uses Wren S/H lexer.) |
||
| (24 intermediate revisions by the same user not shown) | |||
Line 1:
==Arbitrary precision arithmetic==
Although this is already supported in an easy to use way by the [
I have therefore written a Wren wrapper for the C library, GMP, to deal with such cases. This supports arithmetic not just on integers and rational numbers (as Wren-big does) but also on floating point numbers of arbitrary size. Most of GMP's functions (around 190) have been wrapped though some - notably those concerned with input/output and random numbers - have been excluded as I didn't think they would be very useful from Wren's perspective. I have also included some additional 'convenience' methods which can be coded in a few lines using the GMP functions as well as some prime factorization routines.
Line 23:
However, it's also far more complicated to wrap and doesn't integrate quite so well as MPF with the rest of GMP - as well as an additional rounding mode it uses a default value of NaN rather than zero for new objects which is not ideal from Wren's perspective.
I have therefore decided to stick with MPF for basic arithmetic, for which it is perfectly adequate, but use MPFR for the transcendental functions. There are
===How fast?===
Line 32:
==Source code (Wren)==
<
import "./trait" for Comparable
Line 52:
static four { from(4) }
static five { from(5) }
static six { from(6) }
static seven { from(7) }
static eight { from(8) }
static nine { from(9) }
static ten { from(10) }
Line 102 ⟶ 106:
// Returns the smaller of two Mpz objects.
static min(op1, op2) {
if (!(op1 is Mpz)) op1 = from(op1)
if (!(op2 is Mpz)) op2 = from(op2)
if (op1 < op2) return op1
return op2
Line 108 ⟶ 114:
// Returns the greater of two Mpz objects.
static max(op1, op2) {
if (!(op1 is Mpz)) op1 = from(op1)
if (!(op2 is Mpz)) op2 = from(op2)
if (op1 > op2) return op1
return op2
Line 114 ⟶ 122:
// Returns the positive difference of two Mpz objects.
static dim(op1, op2) {
if (!(op1 is Mpz)) op1 = from(op1)
if (!(op2 is Mpz)) op2 = from(op2)
if (op1 >= op2) return op1 - op2
return Mpz.zero
Line 173 ⟶ 183:
// Converts the current instance to:
foreign toNum
foreign toString(b)
toString { toString(10) }
toMpq { Mpq.
toMpf { Mpf.
/* Methods which assign their result to the current instance ('this').
Line 208 ⟶ 218:
foreign tdiv(n, d) // divides one Mpz object by another (rounds towards zero)
foreign tdivUi(n, d) // divides an Mpz object by a uint (rounds towards zero)
addSi(op1, op2) { (op2 >= 0) ? addUi(op1, op2) : subUi(op1, -op2) } // op2 is a sint
subSi(op1, op2) { (op2 >= 0) ? subUi(op1, op2) : addUi(op1, -op2) } // op2 is a sint
div(n, d) { tdiv(n, d) } // alias for tdiv
divUi(n, d) { tdivUi(n, d) } // alias for tdivUi
divSi(n, d) { (d >= 0) ? tdivUi(n, d) : tdivUi(n, -d).neg } // d is a sint
foreign crem(n, d) // sets the remainder after 'cdiv' by another Mpz object
Line 222 ⟶ 236:
rem(n, d) { trem(n,d) } // alias for trem
remUi(n, d) { tremUi(n, d) }
remSi(n, d) { (d >= 0) ? tremUi(n, d) : tremUi(n, -d) } // d is a sint
foreign mod(n, d) // sets to n (Mpz) mod d (Mpz) ignoring the sign of d
Line 262 ⟶ 277:
foreign modInv(op1, op2) // modular inverse of op1 mod op2
foreign remove(op, f) // removes all factors f from op and returns number removed
foreign nextPrime(op) // next prime > op (Mpz)
prevPrime(op) { // previous prime < op (Mpz) or aborts if no such prime
if (op < Mpz.three) {
Mpz.abort_()
}
if (op == Mpz.three) {
this.setMpz(Mpz.two)
return
}
var n = Mpz.new().setMpz(op.isEven ? op - Mpz.one : op - Mpz.two)
while (true) {
if (n.probPrime(15) > 0) {
this.setMpz(n)
return
}
n.sub(Mpz.two)
}
}
/* Convenience versions of the above methods where the first (or only) argument is 'this'.
Unless otherwise noted, any other argument must either be another Mpz object or a uint.
Other Nums
add(op) { (op is Mpz) ? add(this, op) : ((op is Num) ?
sub(op) { (op is Mpz) ? sub(this, op) : ((op is Num) ?
mul(op) { // sint
if (op is Mpz) return mul(this, op)
if (op is Num) {
Line 283 ⟶ 316:
fdiv(d) { (d is Mpz) ? fdiv(this, d) : ((d is Num) ? fdivUi(this, d) : Mpz.abort_()) }
tdiv(d) { (d is Mpz) ? tdiv(this, d) : ((d is Num) ? tdivUi(this, d) : Mpz.abort_()) }
div(d) { (d is Mpz) ? tdiv(this, d)
crem(d) { (d is Mpz) ? crem(this, d) : ((d is Num) ? cremUi(this, d) : Mpz.abort_()) }
frem(d) { (d is Mpz) ? frem(this, d) : ((d is Num) ? fremUi(this, d) : Mpz.abort_()) }
trem(d) { (d is Mpz) ? trem(this, d) : ((d is Num) ? tremUi(this, d) : Mpz.abort_()) }
rem(d) { (d is Mpz) ? trem(this, d) : ((d is Num) ? remSi (this, d) : Mpz.abort_()) } //
mod(d) { (d is Mpz) ? mod (this, d) : ((d is Num) ? modUi (this, d) : Mpz.abort_()) }
Line 325 ⟶ 358:
modInv(op) { modInv(this, op) } // op is an Mpz object
remove(f) { remove(this, f) } // f is an Mpz object
nextPrime { nextPrime(this) }
prevPrime { prevPrime(this) }
/* As above methods where the first (or only) argument is 'this'
Line 334 ⟶ 368:
~ { copy().com }
+(op) { copy().add(op)
-(op) { copy().sub(op)
*(op) { copy().mul(op)
/(op) { copy().
%(op) { copy().
<<(b) { copy().lsh(b)
>>(b) { copy().rsh(b)
&(op) { copy().and(op)
|(op) { copy().ior(op)
^(op) { copy().xor(op)
/* Methods which may mutate the current instance or simply return it. */
Line 351 ⟶ 385:
max(op) { (op > this) ? set(op) : this } // maximum of this and op
clamp(min, max) { // clamps this to the interval [min, max]
if (
if (this < min) set(min) else if (this > max)
return this.copy()
}
Line 401 ⟶ 435:
foreign isSquare // true if 'this' = a ^ 2, for some integer a
foreign probPrime(reps) //
// or definitely prime respectively after 'reps' repetitions
foreign sign(op) // the sign of an Mpz object
Line 414 ⟶ 449:
foreign fibonacci(n) // the n'th (uint) Fibonacci number
foreign lucas(n) // the n'th (uint) Lucas number
multinomial(n, f) { // the multinomial coefficient of n over a list f where sum(f) == n
if (!(n is Num && n >= 0)) {
Fiber.abort("First argument must be a non-negative integer.")
}
if (!(f is List)) Fiber.abort("Second argument must be a list.")
var sum = f.reduce { |acc, i| acc + i }
if (n != sum) {
Fiber.abort("The elements of the list must sum to 'n'.")
}
var prod = Mpz.one
var fact = Mpz.new()
for (e in f) {
if (e < 0) Fiber.abort("The elements of the list must be non-negative integers.")
if (e > 1) prod.mul(fact.factorial(e))
}
return factorial(n).div(prod)
}
foreign popCount // number of 1 bits in this
Line 423 ⟶ 476:
foreign tstBit(index) // test bit at 'index' in this and return its value
foreign sizeInBase(base) // number of digits in the given base 2..62 ignoring any sign
// if the base is not a power of 2 may be one too big
foreign digitsInBase(base) // as sizeInBase, always exact but slower
sizeInBits { sizeInBase(2) } // number of binary digits in this ignoring any sign
/* Prime factorization methods. */
// Private worker method for Pollard's Rho algorithm.
static pollardRho_(m, seed, c) {
var g = Fn.new { |x|
var x = Mpz.from(seed)
var y = Mpz.from(seed)
Line 456 ⟶ 511:
// Returns Mpz.zero in the event of failure.
static pollardRho(m, seed, c) {
if (m < 2) return
if (m is Num) m = Mpz.from(m)
if (m.probPrime(15) > 0) return m.copy()
if (m.isSquare) return m.copy().sqrt
for (p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]) {
Line 475 ⟶ 530:
var factors = []
var k = Mpz.from(37)
var i = 0
while (k * k <= n) {
if (n.isDivisible(k)) {
Line 489 ⟶ 544:
}
// Private worker method (recursive) to obtain the prime factors of a number (Mpz).
static primeFactors_(m, trialDivs) {
if (m.probPrime(15) > 0) return [m.copy()]
Line 508 ⟶ 563:
while (n > 1) {
if (checkPrime && n.probPrime(15) > 0) {
factors.add(n.copy())
break
}
if (n >= threshold) {
var d = pollardRho_(n, seed, c)
if (d != 0) {
Line 528 ⟶ 583:
checkPrime = false
}
} else if (c < 101) {
c = c + 1
Line 536 ⟶ 590:
factors.addAll(primeFactorsWheel_(n))
break
}
} else {
factors.addAll(primeFactorsWheel_(n))
Line 548 ⟶ 602:
// Returns a list of the primes factors (Mpz) of 'm' (an Mpz object or an integral Num)
// using the wheel based factorization and/or Pollard's Rho algorithm as appropriate.
static primeFactors(m) {
if (m < 2) return []
Line 554 ⟶ 607:
return primeFactors_(m, true)
}
// Returns all the divisors of 'n' including 1 and 'n' itself.
static divisors(n) {
if (n < 1) return []
if (n is Num) n = from(n)
var divs = []
var divs2 = []
var i = one
var k = n.isEven ? one : two
var sqrt = n.copy().sqrt
while (i <= sqrt) {
if (n.isDivisible(i)) {
divs.add(i.copy())
var j = n / i
if (j != i) divs2.add(j)
}
i.add(k)
}
if (!divs2.isEmpty) divs = divs + divs2[-1..0]
return divs
}
// Returns all the divisors of 'n' excluding 'n'.
static properDivisors(n) {
var d = divisors(n)
var c = d.count
return (c <= 1) ? [] : d[0..-2]
}
// As 'divisors' method but uses a different algorithm.
// Better for large numbers with a small number of prime factors.
static divisors2(n) {
if (n is Num) n = from(n)
var pf = primeFactors(n)
if (pf.count == 0) return (n == 1) ? [one] : pf
var arr = []
if (pf.count == 1) {
arr.add([pf[0].copy(), 1])
} else {
var prevPrime = pf[0]
var count = 1
for (i in 1...pf.count) {
if (pf[i] == prevPrime) {
count = count + 1
} else {
arr.add([prevPrime.copy(), count])
prevPrime = pf[i]
count = 1
}
}
arr.add([prevPrime.copy(), count])
}
var divisors = []
var generateDivs
generateDivs = Fn.new { |currIndex, currDivisor|
if (currIndex == arr.count) {
divisors.add(currDivisor.copy())
return
}
for (i in 0..arr[currIndex][1]) {
generateDivs.call(currIndex+1, currDivisor)
currDivisor = currDivisor * arr[currIndex][0]
}
}
generateDivs.call(0, one)
return divisors.sort()
}
// As 'properDivisors' but uses 'divisors2' method.
static properDivisors2(n) {
var d = divisors2(n)
var c = d.count
return (c <= 1) ? [] : d[0..-2]
}
// Returns the sum of all the divisors of 'n' including 1 and 'n' itself.
static divisorSum(n) {
if (n < 1) return zero
n = from(n)
var total = one
var power = two
while (n.isDivisible(two)) {
total.add(power)
power.mul(2)
n.div(2)
}
var i = three
while (i * i <= n) {
var sum = one
power.set(i)
while (n.isDivisible(i)) {
sum.add(power)
power.mul(i)
n.div(i)
}
total.mul(sum)
i.add(2)
}
if (n > 1) total.mul(n + 1)
return total
}
// Returns the number of divisors of 'n' including 1 and 'n' itself.
static divisorCount(n) {
if (n < 1) return 0
n = from(n)
var count = 0
var prod = 1
while (n.isDivisible(two)) {
count = count + 1
n.div(2)
}
prod = prod * (1 + count)
var i = three
while (i * i <= n) {
count = 0
while (n.isDivisible(i)) {
count = count + 1
n.div(i)
}
prod = prod * (1 + count)
i.add(2)
}
if (n > 2) prod = prod * 2
return prod
}
/* Methods which apply to a sequence of Mpz objects. */
Line 575 ⟶ 755:
foreign class Mpq is Comparable {
// Creates small commonly used
static minusOne { from(-1, 1) }
static zero { new() }
Line 583 ⟶ 763:
static four { from(4, 1) }
static five { from(5, 1) }
static six { from(6, 1) }
static seven { from(7, 1) }
static eight { from(8, 1) }
static nine { from(9, 1) }
static ten { from(10, 1) }
static half { from(1, 2) }
Line 588 ⟶ 772:
static quarter { from(1, 4) }
static fifth { from(1, 5) }
static sixth { from(1, 6) }
static seventh { from(1, 7) }
static eighth { from(1, 8) }
static ninth { from(1, 9) }
static tenth { from(1, 10) }
Line 595 ⟶ 783:
// Returns the smaller of two Mpq objects.
static min(op1, op2) {
if (!(op1 is Mpq)) op1 = from(op1)
if (!(op2 is Mpq)) op2 = from(op2)
if (op1 < op2) return op1
return op2
Line 601 ⟶ 791:
// Returns the greater of two Mpq objects.
static max(op1, op2) {
if (!(op1 is Mpq)) op1 = from(op1)
if (!(op2 is Mpq)) op2 = from(op2)
if (op1 > op2) return op1
return op2
Line 607 ⟶ 799:
// Returns the positive difference of two Mpq objects.
static dim(op1, op2) {
if (!(op1 is Mpq)) op1 = from(op1)
if (!(op2 is Mpq)) op2 = from(op2)
if (op1 >= op2) return op1 - op2
return Mpq.zero
Line 692 ⟶ 886:
foreign abs(op) // sets to the absolute value of op (Mpq)
foreign inv(op) // sets to 1/op (Mpq)
// rounding methods, similar to those in Mpf class
ceil(op) { // higher integer
if (op.den == Mpz.one) return setMpq(op)
var div = op.num / op.den
if (op.sign >= 0) div.inc
return setMpz(div)
}
floor(op) { // lower integer
if (op.den == Mpz.one) return setMpq(op)
var div = op.num / op.den
if (op.sign < 0) div.dec
return setMpz(div)
}
trunc(op) { op.sign < 0 ? ceil(op) : floor(op) } // lower integer, towards zero
round(op) { // nearer integer, 1/2 away from zero
if (op.den == Mpz.one) return setMpq(op)
var op2 = op.copy()
if (op2.sign < 0) op2.sub(Mpq.half) else op2.add(Mpq.half)
return trunc(op2)
}
roundUp (op) { op.sign >= 0 ? ceil(op) : floor(op) } // higher integer, away from zero
/* As above methods where the first (or only) argument is 'this'. */
Line 704 ⟶ 924:
abs { abs(this) }
inv { inv(this) }
ceil { ceil(this) }
floor { floor(this) }
trunc { trunc(this) }
round { round(this) }
roundUp { roundUp(this) }
/* As above methods where the first (or only) argument is 'this'
Line 723 ⟶ 949:
max(op) { (op > this) ? set(op) : this } // maximum of this and op
clamp(min, max) { // clamps this to the interval [min, max]
if (
if (this < min) set(min) else if (this > max)
return this.copy()
}
Line 754 ⟶ 980:
sign { sign(this) } // the sign of the current instance
isZero { cmpUi(0) == 0 } // returns 'true' if the current instance is zero
fraction { // the fractional part of this with the same sign, as a new Mpq object
var t = copy().trunc
return t.sub(this, t)
}
/* Methods which apply to a sequence of Mpq objects. */
Line 778 ⟶ 1,009:
// Creates small commonly used Mpf objects with default precision.
static minusOne { from(-1) }
static zero { new() }
static one { from(1) }
static two { from(2) }
static three { from(3) }
static four { from(4) }
static five { from(5) }
static
static
static eight { from(8) }
static nine { from(9) }
static ten { from(10) }
static half { from(0.5) }
static third { new().inv(three) }
static quarter { from(0.25) }
static fifth { from(0.2) }
static
static seventh { new().inv(seven) }
static eighth { from(0.125) }
static ninth { new().inv(nine) }
static tenth { from(0.1) }
// Creates transcendental constants with a specified precision.
Line 801 ⟶ 1,039:
static phi(prec) { from(5, prec).sqrt.add(1).div(2) } // Phi (golden ratio)
static sqrt2(prec) { from(2, prec).sqrt } // Square root of 2
static euler(prec) { from(1, prec).digamma.neg } // Euler's constant
// Convenience versions of the above constants which use the default precision.
static e { from(1).exp }
static ln2 { from(2).log }
static ln10 { from(10).log }
static pi { new().acos(Mpf.zero).mul(2) }
static tau { new().acos(Mpf.zero).mul(4) }
static phi { from(5).sqrt.add(1).div(2) }
static sqrt2 { from(2).sqrt }
static euler { from(1).digamma.neg }
// Swaps the values of two Mpf objects.
Line 815 ⟶ 1,064:
// Returns the smaller of two Mpf objects.
static min(op1, op2) {
if (!(op1 is Mpf)) op1 = from(op1)
if (!(op2 is Mpf)) op2 = from(op2)
if (op1 < op2) return op1
return op2
Line 821 ⟶ 1,072:
// Returns the greater of two Mpf objects.
static max(op1, op2) {
if (!(op1 is Mpf)) op1 = from(op1)
if (!(op2 is Mpf)) op2 = from(op2)
if (op1 > op2) return op1
return op2
Line 827 ⟶ 1,080:
// Returns the positive difference of two Mpf objects.
static dim(op1, op2) {
if (!(op1 is Mpf)) op1 = from(op1)
if (!(op2 is Mpf)) op2 = from(op2)
if (op1 >= op2) return op1 - op2
return Mpf.zero
Line 946 ⟶ 1,201:
foreign uiDiv(op1, op2) // divides a uint by an Mpf object
foreign div2(op1, op2) // divides op1 (Mpf) by 2^op2 (uint)
addSi(op1, op2) { (op2 >= 0) ? addUi(op1, op2) : subUi(op1, -op2) } // op2 is a sint
subSi(op1, op2) { (op2 >= 0) ? subUi(op1, op2) : addUi(op1, -op2) } // op2 is a sint
mulSi(op1, op2) { (op2 >= 0) ? mulUi(op1, op2) : mulUi(op1, -op2).neg } // op2 is a sint
divSi(op1, op2) { (op2 >= 0) ? divUi(op1, op2) : divUi(op1, -op2).neg } // op2 is a sint
foreign neg(op) // sets to -op (Mpf)
Line 993 ⟶ 1,253:
foreign atan(op) // arc-tangent of op
foreign atan2(y, x) // arc-tangent 2 of y and x
foreign cosh(op) // hyperbolic cosine of op
foreign sinh(op) // hyperbolic sine of op
foreign tanh(op) // hyperbolic tangent of op
foreign gamma(op) // gamma function of op
foreign gammaInc(op, op2) // incomplete gamma function of op and op2
foreign lgamma(op) // natural logarithm of the absolute value of gamma(op)
foreign digamma(op) // digamma (or psi) function of op
foreign zeta(op) // zeta function of op (Mpf)
foreign zetaUi(op) // zeta function of op (uint)
/* As above methods where the first (or only) argument is 'this'.
Unless otherwise noted, any other argument must either be another Mpf object or a
Other Nums are converted by GMP to uint. */
add(op) { (op is Mpf) ? add(this, op) : ((op is Num) ?
sub(op) { (op is Mpf) ? sub(this, op) : ((op is Num) ?
mul(op) { (op is Mpf) ? mul(this, op) : ((op is Num) ?
mul2(op) { mul2(this, op) }
div(op) { (op is Mpf) ? div(this, op) : ((op is Num) ?
div2(op) { div2(this, op) }
Line 1,016 ⟶ 1,285:
pow(exp) { pow(this, exp) } // exp must be a uint
powz(exp) { powz(this, exp) } // exp must be an Mpz object
powf(exp) {
square { mul(this, this) }
cube { pow(this, 3) }
Line 1,042 ⟶ 1,311:
atan { atan(this) }
atan2(x) { atan2(this, x) }
cosh { cosh(this) }
sinh { sinh(this) }
tanh { tanh(this) }
gamma { gamma(this) }
gammaInc(op) { gammaInc(this, op) }
lgamma { lgamma(this) }
digamma { digamma(this) }
zeta { zeta(this) }
/* As above methods where the first (or only) argument is 'this'
Line 1,061 ⟶ 1,339:
max(op) { (op > this) ? set(op) : this } // maximum of this and op
clamp(min, max) { // clamps this to the interval [min, max]
if (
if (this < min) set(min) else if (this > max)
return this.copy()
}
Line 1,092 ⟶ 1,370:
isZero { cmpUi(0) == 0 } // returns 'true' if the current instance is zero
fraction {
var t = copy().trunc
return t.sub(this, t)
Line 1,111 ⟶ 1,389:
static min(sf) { sf.reduce { |acc, x| (x > acc) ? x : acc } }
static max(sf) { sf.reduce { |acc, x| (x < acc) ? x : acc } }
}</
==Source code (C)==
<
#include <stdio.h>
Line 1,981 ⟶ 2,260:
size_t ret = mpz_sizeinbase(*pop, base);
wrenSetSlotDouble(vm, 0, (double)ret);
}
void Mpz_digitsInBase(WrenVM* vm) {
const mpz_t *pop = (const mpz_t*)wrenGetSlotForeign(vm, 0);
int base = (int)wrenGetSlotDouble(vm, 1);
char *ret = mpz_get_str(NULL, base, *pop);
size_t digits = strlen(ret);
if (ret[0] == '-') --digits;
wrenSetSlotDouble(vm, 0, (double)digits);
free(ret);
}
Line 2,769 ⟶ 3,058:
mpfr_clear(ry);
mpfr_clear(rx);
mpfr_clear(res);
}
void Mpf_cosh(WrenVM* vm) {
mpfr_t op, res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_cosh(res, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(res);
}
void Mpf_sinh(WrenVM* vm) {
mpfr_t op, res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_sinh(res, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(res);
}
void Mpf_tanh(WrenVM* vm) {
mpfr_t op, res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_tanh(res, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(res);
}
void Mpf_gamma(WrenVM* vm) {
mpfr_t op, res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_gamma(res, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(res);
}
void Mpf_gammaInc(WrenVM* vm) {
mpfr_t op, op2, res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpf_t *pop2 = (mpf_t*)wrenGetSlotForeign(vm, 2);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_prec_t prec3 = (mpfr_prec_t)mpf_get_prec(*pop2);
mpfr_init2(op2, prec3);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_set_f(op2, *pop2, MPFR_RNDN);
mpfr_gamma_inc(res, op, op2, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(op2);
mpfr_clear(res);
}
void Mpf_lgamma(WrenVM* vm) {
mpfr_t op, res;
int signnp = 0;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_lgamma(res, &signnp, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(res);
}
void Mpf_digamma(WrenVM* vm) {
mpfr_t op, res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_digamma(res, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(res);
}
void Mpf_zeta(WrenVM* vm) {
mpfr_t op, res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
mpf_t *pop = (mpf_t*)wrenGetSlotForeign(vm, 1);
mpfr_prec_t prec1 = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec1);
mpfr_prec_t prec2 = (mpfr_prec_t)mpf_get_prec(*pop);
mpfr_init2(op, prec2);
mpfr_set_f(op, *pop, MPFR_RNDN);
mpfr_zeta(res, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(op);
mpfr_clear(res);
}
void Mpf_zetaUi(WrenVM* vm) {
mpfr_t res;
mpf_t *pf = (mpf_t*)wrenGetSlotForeign(vm, 0);
unsigned long op = (unsigned long)wrenGetSlotDouble(vm, 1);
mpfr_prec_t prec = (mpfr_prec_t)mpf_get_prec(*pf);
mpfr_init2(res, prec);
mpfr_zeta_ui(res, op, MPFR_RNDN);
mpfr_get_f(*pf, res, MPFR_RNDN);
mpfr_clear(res);
}
Line 2,908 ⟶ 3,334:
if(!isStatic && strcmp(signature, "tstBit(_)") == 0) return Mpz_tstBit;
if(!isStatic && strcmp(signature, "sizeInBase(_)") == 0) return Mpz_sizeInBase;
if(!isStatic && strcmp(signature, "digitsInBase(_)") == 0) return Mpz_digitsInBase;
} else if (strcmp(className, "Mpq") == 0) {
if(isStatic && strcmp(signature, "swap(_,_)") == 0) return Mpq_swap;
Line 3,003 ⟶ 3,430:
if(!isStatic && strcmp(signature, "atan(_)") == 0) return Mpf_atan;
if(!isStatic && strcmp(signature, "atan2(_,_)") == 0) return Mpf_atan2;
if(!isStatic && strcmp(signature, "cosh(_)") == 0) return Mpf_cosh;
if(!isStatic && strcmp(signature, "sinh(_)") == 0) return Mpf_sinh;
if(!isStatic && strcmp(signature, "tanh(_)") == 0) return Mpf_tanh;
if(!isStatic && strcmp(signature, "gamma(_)") == 0) return Mpf_gamma;
if(!isStatic && strcmp(signature, "gammaInc(_,_)") == 0) return Mpf_gammaInc;
if(!isStatic && strcmp(signature, "lgamma(_)") == 0) return Mpf_lgamma;
if(!isStatic && strcmp(signature, "digamma(_)") == 0) return Mpf_digamma;
if(!isStatic && strcmp(signature, "zeta(_)") == 0) return Mpf_zeta;
if(!isStatic && strcmp(signature, "zetaUi(_)") == 0) return Mpf_zetaUi;
}
}
Line 3,086 ⟶ 3,522:
mpfr_free_cache();
return 0;
}</
| |||