| /* |
| * QuickJS Javascript Calculator |
| * |
| * Copyright (c) 2017-2020 Fabrice Bellard |
| * |
| * Permission is hereby granted, free of charge, to any person obtaining a copy |
| * of this software and associated documentation files (the "Software"), to deal |
| * in the Software without restriction, including without limitation the rights |
| * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
| * copies of the Software, and to permit persons to whom the Software is |
| * furnished to do so, subject to the following conditions: |
| * |
| * The above copyright notice and this permission notice shall be included in |
| * all copies or substantial portions of the Software. |
| * |
| * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
| * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
| * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
| * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
| * THE SOFTWARE. |
| */ |
| "use strict"; |
| "use math"; |
| |
| var Integer, Float, Fraction, Complex, Mod, Polynomial, PolyMod, RationalFunction, Series, Matrix; |
| |
| (function(global) { |
| global.Integer = global.BigInt; |
| global.Float = global.BigFloat; |
| global.algebraicMode = true; |
| |
| /* add non enumerable properties */ |
| function add_props(obj, props) { |
| var i, val, prop, tab, desc; |
| tab = Reflect.ownKeys(props); |
| for(i = 0; i < tab.length; i++) { |
| prop = tab[i]; |
| desc = Object.getOwnPropertyDescriptor(props, prop); |
| desc.enumerable = false; |
| if ("value" in desc) { |
| if (typeof desc.value !== "function") { |
| desc.writable = false; |
| desc.configurable = false; |
| } |
| } else { |
| /* getter/setter */ |
| desc.configurable = false; |
| } |
| Object.defineProperty(obj, prop, desc); |
| } |
| } |
| |
| /* same as proto[Symbol.operatorSet] = Operators.create(..op_list) |
| but allow shortcuts: left: [], right: [] or both |
| */ |
| function operators_set(proto, ...op_list) |
| { |
| var new_op_list, i, a, j, b, k, obj, tab; |
| var fields = [ "left", "right" ]; |
| new_op_list = []; |
| for(i = 0; i < op_list.length; i++) { |
| a = op_list[i]; |
| if (a.left || a.right) { |
| tab = [ a.left, a.right ]; |
| delete a.left; |
| delete a.right; |
| for(k = 0; k < 2; k++) { |
| obj = tab[k]; |
| if (obj) { |
| if (!Array.isArray(obj)) { |
| obj = [ obj ]; |
| } |
| for(j = 0; j < obj.length; j++) { |
| b = {}; |
| Object.assign(b, a); |
| b[fields[k]] = obj[j]; |
| new_op_list.push(b); |
| } |
| } |
| } |
| } else { |
| new_op_list.push(a); |
| } |
| } |
| proto[Symbol.operatorSet] = |
| Operators.create.call(null, ...new_op_list); |
| } |
| |
| /* Integer */ |
| |
| function generic_pow(a, b) { |
| var r, is_neg, i; |
| if (!Integer.isInteger(b)) { |
| return exp(log(a) * b); |
| } |
| if (Array.isArray(a) && !(a instanceof Polynomial || |
| a instanceof Series)) { |
| r = idn(Matrix.check_square(a)); |
| } else { |
| r = 1; |
| } |
| if (b == 0) |
| return r; |
| is_neg = false; |
| if (b < 0) { |
| is_neg = true; |
| b = -b; |
| } |
| r = a; |
| for(i = Integer.floorLog2(b) - 1; i >= 0; i--) { |
| r *= r; |
| if ((b >> i) & 1) |
| r *= a; |
| } |
| if (is_neg) { |
| if (typeof r.inverse != "function") |
| throw "negative powers are not supported for this type"; |
| r = r.inverse(); |
| } |
| return r; |
| } |
| |
| var small_primes = [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499 ]; |
| |
| function miller_rabin_test(n, t) { |
| var d, r, s, i, j, a; |
| d = n - 1; |
| s = 0; |
| while ((d & 1) == 0) { |
| d >>= 1; |
| s++; |
| } |
| if (small_primes.length < t) |
| t = small_primes.length; |
| loop: for(j = 0; j < t; j++) { |
| a = small_primes[j]; |
| r = Integer.pmod(a, d, n); |
| if (r == 1 || r == (n - 1)) |
| continue; |
| for(i = 1; i < s; i++) { |
| r = (r * r) % n; |
| if (r == 1) |
| return false; |
| if (r == (n - 1)) |
| continue loop; |
| } |
| return false; /* n is composite */ |
| } |
| return true; /* n is probably prime with probability (1-0.5^t) */ |
| } |
| |
| function fact_rec(a, b) { /* assumes a <= b */ |
| var i, r; |
| if ((b - a) <= 5) { |
| r = a; |
| for(i = a + 1; i <= b; i++) |
| r *= i; |
| return r; |
| } else { |
| /* to avoid a quadratic running time it is better to |
| multiply numbers of similar size */ |
| i = (a + b) >> 1; |
| return fact_rec(a, i) * fact_rec(i + 1, b); |
| } |
| } |
| |
| /* math mode specific quirk to overload the integer division and power */ |
| Operators.updateBigIntOperators( |
| { |
| "/"(a, b) { |
| if (algebraicMode) { |
| return Fraction.toFraction(a, b); |
| } else { |
| return Float(a) / Float(b); |
| } |
| }, |
| "**"(a, b) { |
| if (algebraicMode) { |
| return generic_pow(a, b); |
| } else { |
| return Float(a) ** Float(b); |
| } |
| } |
| }); |
| |
| add_props(Integer, { |
| isInteger(a) { |
| /* integers are represented either as bigint or as number */ |
| return typeof a === "bigint" || |
| (typeof a === "number" && Number.isSafeInteger(a)); |
| }, |
| gcd(a, b) { |
| var r; |
| while (b != 0) { |
| r = a % b; |
| a = b; |
| b = r; |
| } |
| return a; |
| }, |
| fact(n) { |
| return n <= 0 ? 1 : fact_rec(1, n); |
| }, |
| /* binomial coefficient */ |
| comb(n, k) { |
| if (k < 0 || k > n) |
| return 0; |
| if (k > n - k) |
| k = n - k; |
| if (k == 0) |
| return 1; |
| return Integer.tdiv(fact_rec(n - k + 1, n), fact_rec(1, k)); |
| }, |
| /* inverse of x modulo y */ |
| invmod(x, y) { |
| var q, u, v, a, c, t; |
| u = x; |
| v = y; |
| c = 1; |
| a = 0; |
| while (u != 0) { |
| t = Integer.fdivrem(v, u); |
| q = t[0]; |
| v = u; |
| u = t[1]; |
| t = c; |
| c = a - q * c; |
| a = t; |
| } |
| /* v = gcd(x, y) */ |
| if (v != 1) |
| throw RangeError("not invertible"); |
| return a % y; |
| }, |
| /* return a ^ b modulo m */ |
| pmod(a, b, m) { |
| var r; |
| if (b == 0) |
| return 1; |
| if (b < 0) { |
| a = Integer.invmod(a, m); |
| b = -b; |
| } |
| r = 1; |
| for(;;) { |
| if (b & 1) { |
| r = (r * a) % m; |
| } |
| b >>= 1; |
| if (b == 0) |
| break; |
| a = (a * a) % m; |
| } |
| return r; |
| }, |
| |
| /* return true if n is prime (or probably prime with |
| probability 1-0.5^t) */ |
| isPrime(n, t) { |
| var i, d, n1; |
| if (!Integer.isInteger(n)) |
| throw TypeError("invalid type"); |
| if (n <= 1) |
| return false; |
| n1 = small_primes.length; |
| /* XXX: need Integer.sqrt() */ |
| for(i = 0; i < n1; i++) { |
| d = small_primes[i]; |
| if (d == n) |
| return true; |
| if (d > n) |
| return false; |
| if ((n % d) == 0) |
| return false; |
| } |
| if (n < d * d) |
| return true; |
| if (typeof t == "undefined") |
| t = 64; |
| return miller_rabin_test(n, t); |
| }, |
| nextPrime(n) { |
| if (!Integer.isInteger(n)) |
| throw TypeError("invalid type"); |
| if (n < 1) |
| n = 1; |
| for(;;) { |
| n++; |
| if (Integer.isPrime(n)) |
| return n; |
| } |
| }, |
| factor(n) { |
| var r, d; |
| if (!Integer.isInteger(n)) |
| throw TypeError("invalid type"); |
| r = []; |
| if (abs(n) <= 1) { |
| r.push(n); |
| return r; |
| } |
| if (n < 0) { |
| r.push(-1); |
| n = -n; |
| } |
| |
| while ((n % 2) == 0) { |
| n >>= 1; |
| r.push(2); |
| } |
| |
| d = 3; |
| while (n != 1) { |
| if (Integer.isPrime(n)) { |
| r.push(n); |
| break; |
| } |
| /* we are sure there is at least one divisor, so one test */ |
| for(;;) { |
| if ((n % d) == 0) |
| break; |
| d += 2; |
| } |
| for(;;) { |
| r.push(d); |
| n = Integer.tdiv(n, d); |
| if ((n % d) != 0) |
| break; |
| } |
| } |
| return r; |
| }, |
| }); |
| |
| add_props(Integer.prototype, { |
| inverse() { |
| return 1 / this; |
| }, |
| norm2() { |
| return this * this; |
| }, |
| abs() { |
| var v = this; |
| if (v < 0) |
| v = -v; |
| return v; |
| }, |
| conj() { |
| return this; |
| }, |
| arg() { |
| if (this >= 0) |
| return 0; |
| else |
| return Float.PI; |
| }, |
| exp() { |
| if (this == 0) |
| return 1; |
| else |
| return Float.exp(this); |
| }, |
| log() { |
| if (this == 1) |
| return 0; |
| else |
| return Float(this).log(); |
| }, |
| }); |
| |
| /* Fraction */ |
| |
| Fraction = function Fraction(a, b) |
| { |
| var d, r, obj; |
| |
| if (new.target) |
| throw TypeError("not a constructor"); |
| if (a instanceof Fraction) |
| return a; |
| if (!Integer.isInteger(a)) |
| throw TypeError("integer expected"); |
| if (typeof b === "undefined") { |
| b = 1; |
| } else { |
| if (!Integer.isInteger(b)) |
| throw TypeError("integer expected"); |
| if (b == 0) |
| throw RangeError("division by zero"); |
| d = Integer.gcd(a, b); |
| if (d != 1) { |
| a = Integer.tdiv(a, d); |
| b = Integer.tdiv(b, d); |
| } |
| |
| /* the fractions are normalized with den > 0 */ |
| if (b < 0) { |
| a = -a; |
| b = -b; |
| } |
| } |
| obj = Object.create(Fraction.prototype); |
| obj.num = a; |
| obj.den = b; |
| return obj; |
| } |
| |
| function fraction_add(a, b) { |
| a = Fraction(a); |
| b = Fraction(b); |
| return Fraction.toFraction(a.num * b.den + a.den * b.num, a.den * b.den); |
| } |
| function fraction_sub(a, b) { |
| a = Fraction(a); |
| b = Fraction(b); |
| return Fraction.toFraction(a.num * b.den - a.den * b.num, a.den * b.den); |
| } |
| function fraction_mul(a, b) { |
| a = Fraction(a); |
| b = Fraction(b); |
| return Fraction.toFraction(a.num * b.num, a.den * b.den); |
| } |
| function fraction_div(a, b) { |
| a = Fraction(a); |
| b = Fraction(b); |
| return Fraction.toFraction(a.num * b.den, a.den * b.num); |
| } |
| function fraction_mod(a, b) { |
| var a1 = Fraction(a); |
| var b1 = Fraction(b); |
| return a - Integer.ediv(a1.num * b1.den, a1.den * b1.num) * b; |
| } |
| function fraction_eq(a, b) { |
| a = Fraction(a); |
| b = Fraction(b); |
| /* we assume the fractions are normalized */ |
| return (a.num == b.num && a.den == b.den); |
| } |
| function fraction_lt(a, b) { |
| a = Fraction(a); |
| b = Fraction(b); |
| return (a.num * b.den < b.num * a.den); |
| } |
| |
| /* operators are needed for fractions */ |
| function float_add(a, b) { |
| return Float(a) + Float(b); |
| } |
| function float_sub(a, b) { |
| return Float(a) - Float(b); |
| } |
| function float_mul(a, b) { |
| return Float(a) * Float(b); |
| } |
| function float_div(a, b) { |
| return Float(a) / Float(b); |
| } |
| function float_mod(a, b) { |
| return Float(a) % Float(b); |
| } |
| function float_pow(a, b) { |
| return Float(a) ** Float(b); |
| } |
| function float_eq(a, b) { |
| /* XXX: may be better to use infinite precision for the comparison */ |
| return Float(a) === Float(b); |
| } |
| function float_lt(a, b) { |
| a = Float(a); |
| b = Float(b); |
| /* XXX: may be better to use infinite precision for the comparison */ |
| if (Float.isNaN(a) || Float.isNaN(b)) |
| return undefined; |
| else |
| return a < b; |
| } |
| |
| operators_set(Fraction.prototype, |
| { |
| "+": fraction_add, |
| "-": fraction_sub, |
| "*": fraction_mul, |
| "/": fraction_div, |
| "%": fraction_mod, |
| "**": generic_pow, |
| "==": fraction_eq, |
| "<": fraction_lt, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| return Fraction(-a.num, a.den); |
| }, |
| }, |
| { |
| left: [Number, BigInt], |
| right: [Number, BigInt], |
| "+": fraction_add, |
| "-": fraction_sub, |
| "*": fraction_mul, |
| "/": fraction_div, |
| "%": fraction_mod, |
| "**": generic_pow, |
| "==": fraction_eq, |
| "<": fraction_lt, |
| }, |
| { |
| left: Float, |
| right: Float, |
| "+": float_add, |
| "-": float_sub, |
| "*": float_mul, |
| "/": float_div, |
| "%": float_mod, |
| "**": float_pow, |
| "==": float_eq, |
| "<": float_lt, |
| }); |
| |
| add_props(Fraction, { |
| /* (internal use) simplify 'a' to an integer when possible */ |
| toFraction(a, b) { |
| var r = Fraction(a, b); |
| if (algebraicMode && r.den == 1) |
| return r.num; |
| else |
| return r; |
| }, |
| }); |
| |
| add_props(Fraction.prototype, { |
| [Symbol.toPrimitive](hint) { |
| if (hint === "string") { |
| return this.toString(); |
| } else { |
| return Float(this.num) / this.den; |
| } |
| }, |
| inverse() { |
| return Fraction(this.den, this.num); |
| }, |
| toString() { |
| return this.num + "/" + this.den; |
| }, |
| norm2() { |
| return this * this; |
| }, |
| abs() { |
| if (this.num < 0) |
| return -this; |
| else |
| return this; |
| }, |
| conj() { |
| return this; |
| }, |
| arg() { |
| if (this.num >= 0) |
| return 0; |
| else |
| return Float.PI; |
| }, |
| exp() { |
| return Float.exp(Float(this)); |
| }, |
| log() { |
| return Float(this).log(); |
| }, |
| }); |
| |
| /* Number (Float64) */ |
| |
| add_props(Number.prototype, { |
| inverse() { |
| return 1 / this; |
| }, |
| norm2() { |
| return this * this; |
| }, |
| abs() { |
| return Math.abs(this); |
| }, |
| conj() { |
| return this; |
| }, |
| arg() { |
| if (this >= 0) |
| return 0; |
| else |
| return Float.PI; |
| }, |
| exp() { |
| return Float.exp(this); |
| }, |
| log() { |
| if (this < 0) { |
| return Complex(this).log(); |
| } else { |
| return Float.log(this); |
| } |
| }, |
| }); |
| |
| /* Float */ |
| |
| var const_tab = []; |
| |
| /* we cache the constants for small precisions */ |
| function get_const(n) { |
| var t, c, p; |
| t = const_tab[n]; |
| p = BigFloatEnv.prec; |
| if (t && t.prec == p) { |
| return t.val; |
| } else { |
| switch(n) { |
| case 0: c = Float.exp(1); break; |
| case 1: c = Float.log(10); break; |
| // case 2: c = Float.log(2); break; |
| case 3: c = 1/Float.log(2); break; |
| case 4: c = 1/Float.log(10); break; |
| // case 5: c = Float.atan(1) * 4; break; |
| case 6: c = Float.sqrt(0.5); break; |
| case 7: c = Float.sqrt(2); break; |
| } |
| if (p <= 1024) { |
| const_tab[n] = { prec: p, val: c }; |
| } |
| return c; |
| } |
| } |
| |
| add_props(Float, { |
| isFloat(a) { |
| return typeof a === "number" || typeof a === "bigfloat"; |
| }, |
| bestappr(u, b) { |
| var num1, num0, den1, den0, u, num, den, n; |
| |
| if (typeof b === "undefined") |
| throw TypeError("second argument expected"); |
| num1 = 1; |
| num0 = 0; |
| den1 = 0; |
| den0 = 1; |
| for(;;) { |
| n = Integer(Float.floor(u)); |
| num = n * num1 + num0; |
| den = n * den1 + den0; |
| if (den > b) |
| break; |
| u = 1.0 / (u - n); |
| num0 = num1; |
| num1 = num; |
| den0 = den1; |
| den1 = den; |
| } |
| return Fraction(num1, den1); |
| }, |
| /* similar constants as Math.x */ |
| get E() { return get_const(0); }, |
| get LN10() { return get_const(1); }, |
| // get LN2() { return get_const(2); }, |
| get LOG2E() { return get_const(3); }, |
| get LOG10E() { return get_const(4); }, |
| // get PI() { return get_const(5); }, |
| get SQRT1_2() { return get_const(6); }, |
| get SQRT2() { return get_const(7); }, |
| }); |
| |
| add_props(Float.prototype, { |
| inverse() { |
| return 1.0 / this; |
| }, |
| norm2() { |
| return this * this; |
| }, |
| abs() { |
| return Float.abs(this); |
| }, |
| conj() { |
| return this; |
| }, |
| arg() { |
| if (this >= 0) |
| return 0; |
| else |
| return Float.PI; |
| }, |
| exp() { |
| return Float.exp(this); |
| }, |
| log() { |
| if (this < 0) { |
| return Complex(this).log(); |
| } else { |
| return Float.log(this); |
| } |
| }, |
| }); |
| |
| /* Complex */ |
| |
| Complex = function Complex(re, im) |
| { |
| var obj; |
| if (new.target) |
| throw TypeError("not a constructor"); |
| if (re instanceof Complex) |
| return re; |
| if (typeof im === "undefined") { |
| im = 0; |
| } |
| obj = Object.create(Complex.prototype); |
| obj.re = re; |
| obj.im = im; |
| return obj; |
| } |
| |
| |
| function complex_add(a, b) { |
| a = Complex(a); |
| b = Complex(b); |
| return Complex.toComplex(a.re + b.re, a.im + b.im); |
| } |
| function complex_sub(a, b) { |
| a = Complex(a); |
| b = Complex(b); |
| return Complex.toComplex(a.re - b.re, a.im - b.im); |
| } |
| function complex_mul(a, b) { |
| a = Complex(a); |
| b = Complex(b); |
| return Complex.toComplex(a.re * b.re - a.im * b.im, |
| a.re * b.im + a.im * b.re); |
| } |
| function complex_div(a, b) { |
| a = Complex(a); |
| b = Complex(b); |
| return a * b.inverse(); |
| } |
| function complex_eq(a, b) { |
| a = Complex(a); |
| b = Complex(b); |
| return a.re == b.re && a.im == b.im; |
| } |
| |
| operators_set(Complex.prototype, |
| { |
| "+": complex_add, |
| "-": complex_sub, |
| "*": complex_mul, |
| "/": complex_div, |
| "**": generic_pow, |
| "==": complex_eq, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| return Complex(-a.re, -a.im); |
| } |
| }, |
| { |
| left: [Number, BigInt, Float, Fraction], |
| right: [Number, BigInt, Float, Fraction], |
| "+": complex_add, |
| "-": complex_sub, |
| "*": complex_mul, |
| "/": complex_div, |
| "**": generic_pow, |
| "==": complex_eq, |
| }); |
| |
| add_props(Complex, { |
| /* simplify to real number when possible */ |
| toComplex(re, im) { |
| if (algebraicMode && im == 0) |
| return re; |
| else |
| return Complex(re, im); |
| }, |
| }); |
| |
| add_props(Complex.prototype, { |
| inverse() { |
| var c = this.norm2(); |
| return Complex(this.re / c, -this.im / c); |
| }, |
| toString() { |
| var v, s = "", a = this; |
| if (a.re != 0) |
| s += a.re.toString(); |
| if (a.im == 1) { |
| if (s != "") |
| s += "+"; |
| s += "I"; |
| } else if (a.im == -1) { |
| s += "-I"; |
| } else { |
| v = a.im.toString(); |
| if (v[0] != "-" && s != "") |
| s += "+"; |
| s += v + "*I"; |
| } |
| return s; |
| }, |
| norm2() { |
| return this.re * this.re + this.im * this.im; |
| }, |
| abs() { |
| return Float.sqrt(norm2(this)); |
| }, |
| conj() { |
| return Complex(this.re, -this.im); |
| }, |
| arg() { |
| return Float.atan2(this.im, this.re); |
| }, |
| exp() { |
| var arg = this.im, r = this.re.exp(); |
| return Complex(r * cos(arg), r * sin(arg)); |
| }, |
| log() { |
| return Complex(abs(this).log(), atan2(this.im, this.re)); |
| }, |
| }); |
| |
| /* Mod */ |
| |
| Mod = function Mod(a, m) { |
| var obj, t; |
| if (new.target) |
| throw TypeError("not a constructor"); |
| obj = Object.create(Mod.prototype); |
| if (Integer.isInteger(m)) { |
| if (m <= 0) |
| throw RangeError("the modulo cannot be <= 0"); |
| if (Integer.isInteger(a)) { |
| a %= m; |
| } else if (a instanceof Fraction) { |
| return Mod(a.num, m) / a.den; |
| } else { |
| throw TypeError("invalid types"); |
| } |
| } else { |
| throw TypeError("invalid types"); |
| } |
| obj.res = a; |
| obj.mod = m; |
| return obj; |
| }; |
| |
| function mod_add(a, b) { |
| if (!(a instanceof Mod)) { |
| return Mod(a + b.res, b.mod); |
| } else if (!(b instanceof Mod)) { |
| return Mod(a.res + b, a.mod); |
| } else { |
| if (a.mod != b.mod) |
| throw TypeError("different modulo for binary operator"); |
| return Mod(a.res + b.res, a.mod); |
| } |
| } |
| function mod_sub(a, b) { |
| if (!(a instanceof Mod)) { |
| return Mod(a - b.res, b.mod); |
| } else if (!(b instanceof Mod)) { |
| return Mod(a.res - b, a.mod); |
| } else { |
| if (a.mod != b.mod) |
| throw TypeError("different modulo for binary operator"); |
| return Mod(a.res - b.res, a.mod); |
| } |
| } |
| function mod_mul(a, b) { |
| if (!(a instanceof Mod)) { |
| return Mod(a * b.res, b.mod); |
| } else if (!(b instanceof Mod)) { |
| return Mod(a.res * b, a.mod); |
| } else { |
| if (a.mod != b.mod) |
| throw TypeError("different modulo for binary operator"); |
| return Mod(a.res * b.res, a.mod); |
| } |
| } |
| function mod_div(a, b) { |
| if (!(b instanceof Mod)) |
| b = Mod(b, a.mod); |
| return mod_mul(a, b.inverse()); |
| } |
| function mod_eq(a, b) { |
| return (a.mod == b.mod && a.res == b.res); |
| } |
| |
| operators_set(Mod.prototype, |
| { |
| "+": mod_add, |
| "-": mod_sub, |
| "*": mod_mul, |
| "/": mod_div, |
| "**": generic_pow, |
| "==": mod_eq, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| return Mod(-a.res, a.mod); |
| } |
| }, |
| { |
| left: [Number, BigInt, Float, Fraction], |
| right: [Number, BigInt, Float, Fraction], |
| "+": mod_add, |
| "-": mod_sub, |
| "*": mod_mul, |
| "/": mod_div, |
| "**": generic_pow, |
| }); |
| |
| add_props(Mod.prototype, { |
| inverse() { |
| var a = this, m = a.mod; |
| if (Integer.isInteger(m)) { |
| return Mod(Integer.invmod(a.res, m), m); |
| } else { |
| throw TypeError("unsupported type"); |
| } |
| }, |
| toString() { |
| return "Mod(" + this.res + "," + this.mod + ")"; |
| }, |
| }); |
| |
| /* Polynomial */ |
| |
| function polynomial_is_scalar(a) |
| { |
| if (typeof a === "number" || |
| typeof a === "bigint" || |
| typeof a === "bigfloat") |
| return true; |
| if (a instanceof Fraction || |
| a instanceof Complex || |
| a instanceof Mod) |
| return true; |
| return false; |
| } |
| |
| Polynomial = function Polynomial(a) |
| { |
| if (new.target) |
| throw TypeError("not a constructor"); |
| if (a instanceof Polynomial) { |
| return a; |
| } else if (Array.isArray(a)) { |
| if (a.length == 0) |
| a = [ 0 ]; |
| Object.setPrototypeOf(a, Polynomial.prototype); |
| return a.trim(); |
| } else if (polynomial_is_scalar(a)) { |
| a = [a]; |
| Object.setPrototypeOf(a, Polynomial.prototype); |
| return a; |
| } else { |
| throw TypeError("invalid type"); |
| } |
| } |
| |
| function number_need_paren(c) |
| { |
| return !(Integer.isInteger(c) || |
| Float.isFloat(c) || |
| c instanceof Fraction || |
| (c instanceof Complex && c.re == 0)); |
| } |
| |
| /* string for c*X^i */ |
| function monomial_toString(c, i) |
| { |
| var str1; |
| if (i == 0) { |
| str1 = c.toString(); |
| } else { |
| if (c == 1) { |
| str1 = ""; |
| } else if (c == -1) { |
| str1 = "-"; |
| } else { |
| if (number_need_paren(c)) { |
| str1 = "(" + c + ")"; |
| } else { |
| str1 = String(c); |
| } |
| str1 += "*"; |
| } |
| str1 += "X"; |
| if (i != 1) { |
| str1 += "^" + i; |
| } |
| } |
| return str1; |
| } |
| |
| /* find one complex root of 'p' starting from z at precision eps using |
| at most max_it iterations. Return null if could not find root. */ |
| function poly_root_laguerre1(p, z, max_it) |
| { |
| var p1, p2, i, z0, z1, z2, d, t0, t1, d1, d2, e, el, zl; |
| |
| d = p.deg(); |
| if (d == 1) { |
| /* monomial case */ |
| return -p[0] / p[1]; |
| } |
| /* trivial zero */ |
| if (p[0] == 0) |
| return 0.0; |
| |
| p1 = p.deriv(); |
| p2 = p1.deriv(); |
| el = 0.0; |
| zl = 0.0; |
| for(i = 0; i < max_it; i++) { |
| z0 = p.apply(z); |
| if (z0 == 0) |
| return z; /* simple exit case */ |
| |
| /* Ward stopping criteria */ |
| e = abs(z - zl); |
| // print("e", i, e); |
| if (i >= 2 && e >= el) { |
| if (abs(zl) < 1e-4) { |
| if (e < 1e-7) |
| return zl; |
| } else { |
| if (e < abs(zl) * 1e-3) |
| return zl; |
| } |
| } |
| el = e; |
| zl = z; |
| |
| z1 = p1.apply(z); |
| z2 = p2.apply(z); |
| t0 = (d - 1) * z1; |
| t0 = t0 * t0; |
| t1 = d * (d - 1) * z0 * z2; |
| t0 = sqrt(t0 - t1); |
| d1 = z1 + t0; |
| d2 = z1 - t0; |
| if (norm2(d2) > norm2(d1)) |
| d1 = d2; |
| if (d1 == 0) |
| return null; |
| z = z - d * z0 / d1; |
| } |
| return null; |
| } |
| |
| function poly_roots(p) |
| { |
| var d, i, roots, j, z, eps; |
| var start_points = [ 0.1, -1.4, 1.7 ]; |
| |
| if (!(p instanceof Polynomial)) |
| throw TypeError("polynomial expected"); |
| d = p.deg(); |
| if (d <= 0) |
| return []; |
| eps = 2.0 ^ (-BigFloatEnv.prec); |
| roots = []; |
| for(i = 0; i < d; i++) { |
| /* XXX: should select another start point if error */ |
| for(j = 0; j < 3; j++) { |
| z = poly_root_laguerre1(p, start_points[j], 100); |
| if (z !== null) |
| break; |
| } |
| if (j == 3) |
| throw RangeError("error in root finding algorithm"); |
| roots[i] = z; |
| p = Polynomial.divrem(p, X - z)[0]; |
| } |
| return roots; |
| } |
| |
| add_props(Polynomial.prototype, { |
| trim() { |
| var a = this, i; |
| i = a.length; |
| while (i > 1 && a[i - 1] == 0) |
| i--; |
| a.length = i; |
| return a; |
| }, |
| conj() { |
| var r, i, n, a; |
| a = this; |
| n = a.length; |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = a[i].conj(); |
| return Polynomial(r); |
| }, |
| inverse() { |
| return RationalFunction(Polynomial([1]), this); |
| }, |
| toString() { |
| var i, str, str1, c, a = this; |
| if (a.length == 1) { |
| return a[0].toString(); |
| } |
| str=""; |
| for(i = a.length - 1; i >= 0; i--) { |
| c = a[i]; |
| if (c == 0 || |
| (c instanceof Mod) && c.res == 0) |
| continue; |
| str1 = monomial_toString(c, i); |
| if (str1[0] != "-") { |
| if (str != "") |
| str += "+"; |
| } |
| str += str1; |
| } |
| return str; |
| }, |
| deg() { |
| if (this.length == 1 && this[0] == 0) |
| return -Infinity; |
| else |
| return this.length - 1; |
| }, |
| apply(b) { |
| var i, n, r, a = this; |
| n = a.length - 1; |
| r = a[n]; |
| while (n > 0) { |
| n--; |
| r = r * b + a[n]; |
| } |
| return r; |
| }, |
| deriv() { |
| var a = this, n, r, i; |
| n = a.length; |
| if (n == 1) { |
| return Polynomial(0); |
| } else { |
| r = []; |
| for(i = 1; i < n; i++) { |
| r[i - 1] = i * a[i]; |
| } |
| return Polynomial(r); |
| } |
| }, |
| integ() { |
| var a = this, n, r, i; |
| n = a.length; |
| r = [0]; |
| for(i = 0; i < n; i++) { |
| r[i + 1] = a[i] / (i + 1); |
| } |
| return Polynomial(r); |
| }, |
| norm2() { |
| var a = this, n, r, i; |
| n = a.length; |
| r = 0; |
| for(i = 0; i < n; i++) { |
| r += a[i].norm2(); |
| } |
| return r; |
| }, |
| }); |
| |
| |
| function polynomial_add(a, b) { |
| var tmp, r, i, n1, n2; |
| a = Polynomial(a); |
| b = Polynomial(b); |
| if (a.length < b.length) { |
| tmp = a; |
| a = b; |
| b = tmp; |
| } |
| n1 = b.length; |
| n2 = a.length; |
| r = []; |
| for(i = 0; i < n1; i++) |
| r[i] = a[i] + b[i]; |
| for(i = n1; i < n2; i++) |
| r[i] = a[i]; |
| return Polynomial(r); |
| } |
| function polynomial_sub(a, b) { |
| return polynomial_add(a, -b); |
| } |
| function polynomial_mul(a, b) { |
| var i, j, n1, n2, n, r; |
| a = Polynomial(a); |
| b = Polynomial(b); |
| n1 = a.length; |
| n2 = b.length; |
| n = n1 + n2 - 1; |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = 0; |
| for(i = 0; i < n1; i++) { |
| for(j = 0; j < n2; j++) { |
| r[i + j] += a[i] * b[j]; |
| } |
| } |
| return Polynomial(r); |
| } |
| function polynomial_div_scalar(a, b) { |
| return a * (1 / b); |
| } |
| function polynomial_div(a, b) |
| { |
| return RationalFunction(Polynomial(a), |
| Polynomial(b)); |
| } |
| function polynomial_mod(a, b) { |
| return Polynomial.divrem(a, b)[1]; |
| } |
| function polynomial_eq(a, b) { |
| var n, i; |
| n = a.length; |
| if (n != b.length) |
| return false; |
| for(i = 0; i < n; i++) { |
| if (a[i] != b[i]) |
| return false; |
| } |
| return true; |
| } |
| |
| operators_set(Polynomial.prototype, |
| { |
| "+": polynomial_add, |
| "-": polynomial_sub, |
| "*": polynomial_mul, |
| "/": polynomial_div, |
| "**": generic_pow, |
| "==": polynomial_eq, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| var r, i, n, a; |
| n = a.length; |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = -a[i]; |
| return Polynomial(r); |
| }, |
| }, |
| { |
| left: [Number, BigInt, Float, Fraction, Complex, Mod], |
| "+": polynomial_add, |
| "-": polynomial_sub, |
| "*": polynomial_mul, |
| "/": polynomial_div, |
| "**": generic_pow, /* XXX: only for integer */ |
| }, |
| { |
| right: [Number, BigInt, Float, Fraction, Complex, Mod], |
| "+": polynomial_add, |
| "-": polynomial_sub, |
| "*": polynomial_mul, |
| "/": polynomial_div_scalar, |
| "**": generic_pow, /* XXX: only for integer */ |
| }); |
| |
| add_props(Polynomial, { |
| divrem(a, b) { |
| var n1, n2, i, j, q, r, n, c; |
| if (b.deg() < 0) |
| throw RangeError("division by zero"); |
| n1 = a.length; |
| n2 = b.length; |
| if (n1 < n2) |
| return [Polynomial([0]), a]; |
| r = Array.prototype.dup.call(a); |
| q = []; |
| n2--; |
| n = n1 - n2; |
| for(i = 0; i < n; i++) |
| q[i] = 0; |
| for(i = n - 1; i >= 0; i--) { |
| c = r[i + n2]; |
| if (c != 0) { |
| c = c / b[n2]; |
| r[i + n2] = 0; |
| for(j = 0; j < n2; j++) { |
| r[i + j] -= b[j] * c; |
| } |
| q[i] = c; |
| } |
| } |
| return [Polynomial(q), Polynomial(r)]; |
| }, |
| gcd(a, b) { |
| var t; |
| while (b.deg() >= 0) { |
| t = Polynomial.divrem(a, b); |
| a = b; |
| b = t[1]; |
| } |
| /* convert to monic form */ |
| return a / a[a.length - 1]; |
| }, |
| invmod(x, y) { |
| var q, u, v, a, c, t; |
| u = x; |
| v = y; |
| c = Polynomial([1]); |
| a = Polynomial([0]); |
| while (u.deg() >= 0) { |
| t = Polynomial.divrem(v, u); |
| q = t[0]; |
| v = u; |
| u = t[1]; |
| t = c; |
| c = a - q * c; |
| a = t; |
| } |
| /* v = gcd(x, y) */ |
| if (v.deg() > 0) |
| throw RangeError("not invertible"); |
| return Polynomial.divrem(a, y)[1]; |
| }, |
| roots(p) { |
| return poly_roots(p); |
| } |
| }); |
| |
| /* Polynomial Modulo Q */ |
| |
| PolyMod = function PolyMod(a, m) { |
| var obj, t; |
| if (new.target) |
| throw TypeError("not a constructor"); |
| obj = Object.create(PolyMod.prototype); |
| if (m instanceof Polynomial) { |
| if (m.deg() <= 0) |
| throw RangeError("the modulo cannot have a degree <= 0"); |
| if (a instanceof RationalFunction) { |
| return PolyMod(a.num, m) / a.den; |
| } else { |
| a = Polynomial(a); |
| t = Polynomial.divrem(a, m); |
| a = t[1]; |
| } |
| } else { |
| throw TypeError("invalid types"); |
| } |
| obj.res = a; |
| obj.mod = m; |
| return obj; |
| }; |
| |
| function polymod_add(a, b) { |
| if (!(a instanceof PolyMod)) { |
| return PolyMod(a + b.res, b.mod); |
| } else if (!(b instanceof PolyMod)) { |
| return PolyMod(a.res + b, a.mod); |
| } else { |
| if (a.mod != b.mod) |
| throw TypeError("different modulo for binary operator"); |
| return PolyMod(a.res + b.res, a.mod); |
| } |
| } |
| function polymod_sub(a, b) { |
| return polymod_add(a, -b); |
| } |
| function polymod_mul(a, b) { |
| if (!(a instanceof PolyMod)) { |
| return PolyMod(a * b.res, b.mod); |
| } else if (!(b instanceof PolyMod)) { |
| return PolyMod(a.res * b, a.mod); |
| } else { |
| if (a.mod != b.mod) |
| throw TypeError("different modulo for binary operator"); |
| return PolyMod(a.res * b.res, a.mod); |
| } |
| } |
| function polymod_div(a, b) { |
| if (!(b instanceof PolyMod)) |
| b = PolyMod(b, a.mod); |
| return polymod_mul(a, b.inverse()); |
| } |
| function polymod_eq(a, b) { |
| return (a.mod == b.mod && a.res == b.res); |
| } |
| |
| operators_set(PolyMod.prototype, |
| { |
| "+": polymod_add, |
| "-": polymod_sub, |
| "*": polymod_mul, |
| "/": polymod_div, |
| "**": generic_pow, |
| "==": polymod_eq, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| return PolyMod(-a.res, a.mod); |
| }, |
| }, |
| { |
| left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial], |
| right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial], |
| "+": polymod_add, |
| "-": polymod_sub, |
| "*": polymod_mul, |
| "/": polymod_div, |
| "**": generic_pow, /* XXX: only for integer */ |
| }); |
| |
| add_props(PolyMod.prototype, { |
| inverse() { |
| var a = this, m = a.mod; |
| if (m instanceof Polynomial) { |
| return PolyMod(Polynomial.invmod(a.res, m), m); |
| } else { |
| throw TypeError("unsupported type"); |
| } |
| }, |
| toString() { |
| return "PolyMod(" + this.res + "," + this.mod + ")"; |
| }, |
| }); |
| |
| /* Rational function */ |
| |
| RationalFunction = function RationalFunction(a, b) |
| { |
| var t, r, d, obj; |
| if (new.target) |
| throw TypeError("not a constructor"); |
| if (!(a instanceof Polynomial) || |
| !(b instanceof Polynomial)) |
| throw TypeError("polynomial expected"); |
| t = Polynomial.divrem(a, b); |
| r = t[1]; |
| if (r.deg() < 0) |
| return t[0]; /* no need for a fraction */ |
| d = Polynomial.gcd(b, r); |
| if (d.deg() > 0) { |
| a = Polynomial.divrem(a, d)[0]; |
| b = Polynomial.divrem(b, d)[0]; |
| } |
| obj = Object.create(RationalFunction.prototype); |
| obj.num = a; |
| obj.den = b; |
| return obj; |
| } |
| |
| add_props(RationalFunction.prototype, { |
| inverse() { |
| return RationalFunction(this.den, this.num); |
| }, |
| conj() { |
| return RationalFunction(this.num.conj(), this.den.conj()); |
| }, |
| toString() { |
| var str; |
| if (this.num.deg() <= 0 && |
| !number_need_paren(this.num[0])) |
| str = this.num.toString(); |
| else |
| str = "(" + this.num.toString() + ")"; |
| str += "/(" + this.den.toString() + ")" |
| return str; |
| }, |
| apply(b) { |
| return this.num.apply(b) / this.den.apply(b); |
| }, |
| deriv() { |
| var n = this.num, d = this.den; |
| return RationalFunction(n.deriv() * d - n * d.deriv(), d * d); |
| }, |
| }); |
| |
| function ratfunc_add(a, b) { |
| a = RationalFunction.toRationalFunction(a); |
| b = RationalFunction.toRationalFunction(b); |
| return RationalFunction(a.num * b.den + a.den * b.num, a.den * b.den); |
| } |
| function ratfunc_sub(a, b) { |
| a = RationalFunction.toRationalFunction(a); |
| b = RationalFunction.toRationalFunction(b); |
| return RationalFunction(a.num * b.den - a.den * b.num, a.den * b.den); |
| } |
| function ratfunc_mul(a, b) { |
| a = RationalFunction.toRationalFunction(a); |
| b = RationalFunction.toRationalFunction(b); |
| return RationalFunction(a.num * b.num, a.den * b.den); |
| } |
| function ratfunc_div(a, b) { |
| a = RationalFunction.toRationalFunction(a); |
| b = RationalFunction.toRationalFunction(b); |
| return RationalFunction(a.num * b.den, a.den * b.num); |
| } |
| function ratfunc_eq(a, b) { |
| a = RationalFunction.toRationalFunction(a); |
| b = RationalFunction.toRationalFunction(b); |
| /* we assume the fractions are normalized */ |
| return (a.num == b.num && a.den == b.den); |
| } |
| |
| operators_set(RationalFunction.prototype, |
| { |
| "+": ratfunc_add, |
| "-": ratfunc_sub, |
| "*": ratfunc_mul, |
| "/": ratfunc_div, |
| "**": generic_pow, |
| "==": ratfunc_eq, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| return RationalFunction(-this.num, this.den); |
| }, |
| }, |
| { |
| left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial], |
| right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial], |
| "+": ratfunc_add, |
| "-": ratfunc_sub, |
| "*": ratfunc_mul, |
| "/": ratfunc_div, |
| "**": generic_pow, /* should only be used with integers */ |
| }); |
| |
| add_props(RationalFunction, { |
| /* This function always return a RationalFunction object even |
| if it could simplified to a polynomial, so it is not |
| equivalent to RationalFunction(a) */ |
| toRationalFunction(a) { |
| var obj; |
| if (a instanceof RationalFunction) { |
| return a; |
| } else { |
| obj = Object.create(RationalFunction.prototype); |
| obj.num = Polynomial(a); |
| obj.den = Polynomial(1); |
| return obj; |
| } |
| }, |
| }); |
| |
| /* Power series */ |
| |
| /* 'a' is an array */ |
| function get_emin(a) { |
| var i, n; |
| n = a.length; |
| for(i = 0; i < n; i++) { |
| if (a[i] != 0) |
| return i; |
| } |
| return n; |
| }; |
| |
| function series_is_scalar_or_polynomial(a) |
| { |
| return polynomial_is_scalar(a) || |
| (a instanceof Polynomial); |
| } |
| |
| /* n is the maximum number of terms if 'a' is not a serie */ |
| Series = function Series(a, n) { |
| var emin, r, i; |
| |
| if (a instanceof Series) { |
| return a; |
| } else if (series_is_scalar_or_polynomial(a)) { |
| if (n <= 0) { |
| /* XXX: should still use the polynomial degree */ |
| return Series.zero(0, 0); |
| } else { |
| a = Polynomial(a); |
| emin = get_emin(a); |
| r = Series.zero(n, emin); |
| n = Math.min(a.length - emin, n); |
| for(i = 0; i < n; i++) |
| r[i] = a[i + emin]; |
| return r; |
| } |
| } else if (a instanceof RationalFunction) { |
| return Series(a.num, n) / a.den; |
| } else { |
| throw TypeError("invalid type"); |
| } |
| }; |
| |
| function series_add(v1, v2) { |
| var tmp, d, emin, n, r, i, j, v2_emin, c1, c2; |
| if (!(v1 instanceof Series)) { |
| tmp = v1; |
| v1 = v2; |
| v2 = tmp; |
| } |
| d = v1.emin + v1.length; |
| if (series_is_scalar_or_polynomial(v2)) { |
| v2 = Polynomial(v2); |
| if (d <= 0) |
| return v1; |
| v2_emin = 0; |
| } else if (v2 instanceof RationalFunction) { |
| /* compute the emin of the rational fonction */ |
| i = get_emin(v2.num) - get_emin(v2.den); |
| if (d <= i) |
| return v1; |
| /* compute the serie with the required terms */ |
| v2 = Series(v2, d - i); |
| v2_emin = v2.emin; |
| } else { |
| v2_emin = v2.emin; |
| d = Math.min(d, v2_emin + v2.length); |
| } |
| emin = Math.min(v1.emin, v2_emin); |
| n = d - emin; |
| r = Series.zero(n, emin); |
| /* XXX: slow */ |
| for(i = emin; i < d; i++) { |
| j = i - v1.emin; |
| if (j >= 0 && j < v1.length) |
| c1 = v1[j]; |
| else |
| c1 = 0; |
| j = i - v2_emin; |
| if (j >= 0 && j < v2.length) |
| c2 = v2[j]; |
| else |
| c2 = 0; |
| r[i - emin] = c1 + c2; |
| } |
| return r.trim(); |
| } |
| function series_sub(a, b) { |
| return series_add(a, -b); |
| } |
| function series_mul(v1, v2) { |
| var n, i, j, r, n, emin, n1, n2, k; |
| if (!(v1 instanceof Series)) |
| v1 = Series(v1, v2.length); |
| else if (!(v2 instanceof Series)) |
| v2 = Series(v2, v1.length); |
| emin = v1.emin + v2.emin; |
| n = Math.min(v1.length, v2.length); |
| n1 = v1.length; |
| n2 = v2.length; |
| r = Series.zero(n, emin); |
| for(i = 0; i < n1; i++) { |
| k = Math.min(n2, n - i); |
| for(j = 0; j < k; j++) { |
| r[i + j] += v1[i] * v2[j]; |
| } |
| } |
| return r.trim(); |
| } |
| function series_div(v1, v2) { |
| if (!(v2 instanceof Series)) |
| v2 = Series(v2, v1.length); |
| return series_mul(v1, v2.inverse()); |
| } |
| function series_pow(a, b) { |
| if (Integer.isInteger(b)) { |
| return generic_pow(a, b); |
| } else { |
| if (!(a instanceof Series)) |
| a = Series(a, b.length); |
| return exp(log(a) * b); |
| } |
| } |
| function series_eq(a, b) { |
| var n, i; |
| if (a.emin != b.emin) |
| return false; |
| n = a.length; |
| if (n != b.length) |
| return false; |
| for(i = 0; i < n; i++) { |
| if (a[i] != b[i]) |
| return false; |
| } |
| return true; |
| } |
| |
| operators_set(Series.prototype, |
| { |
| "+": series_add, |
| "-": series_sub, |
| "*": series_mul, |
| "/": series_div, |
| "**": series_pow, |
| "==": series_eq, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| var obj, n, i; |
| n = a.length; |
| obj = Series.zero(a.length, a.emin); |
| for(i = 0; i < n; i++) { |
| obj[i] = -a[i]; |
| } |
| return obj; |
| }, |
| }, |
| { |
| left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial], |
| right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial], |
| "+": series_add, |
| "-": series_sub, |
| "*": series_mul, |
| "/": series_div, |
| "**": series_pow, |
| }); |
| |
| add_props(Series.prototype, { |
| conj() { |
| var obj, n, i; |
| n = this.length; |
| obj = Series.zero(this.length, this.emin); |
| for(i = 0; i < n; i++) { |
| obj[i] = this[i].conj(); |
| } |
| return obj; |
| }, |
| inverse() { |
| var r, n, i, j, sum, v1 = this; |
| n = v1.length; |
| if (n == 0) |
| throw RangeError("division by zero"); |
| r = Series.zero(n, -v1.emin); |
| r[0] = 1 / v1[0]; |
| for(i = 1; i < n; i++) { |
| sum = 0; |
| for(j = 1; j <= i; j++) { |
| sum += v1[j] * r[i - j]; |
| } |
| r[i] = -sum * r[0]; |
| } |
| return r; |
| }, |
| /* remove leading zero terms */ |
| trim() { |
| var i, j, n, r, v1 = this; |
| n = v1.length; |
| i = 0; |
| while (i < n && v1[i] == 0) |
| i++; |
| if (i == 0) |
| return v1; |
| for(j = i; j < n; j++) |
| v1[j - i] = v1[j]; |
| v1.length = n - i; |
| v1.__proto__.emin += i; |
| return v1; |
| }, |
| toString() { |
| var i, j, str, str1, c, a = this, emin, n; |
| str=""; |
| emin = this.emin; |
| n = this.length; |
| for(j = 0; j < n; j++) { |
| i = j + emin; |
| c = a[j]; |
| if (c != 0) { |
| str1 = monomial_toString(c, i); |
| if (str1[0] != "-") { |
| if (str != "") |
| str += "+"; |
| } |
| str += str1; |
| } |
| } |
| if (str != "") |
| str += "+"; |
| str += "O(" + monomial_toString(1, n + emin) + ")"; |
| return str; |
| }, |
| apply(b) { |
| var i, n, r, a = this; |
| n = a.length; |
| if (n == 0) |
| return 0; |
| r = a[--n]; |
| while (n > 0) { |
| n--; |
| r = r * b + a[n]; |
| } |
| if (a.emin != 0) |
| r *= b ^ a.emin; |
| return r; |
| }, |
| deriv() { |
| var a = this, n = a.length, emin = a.emin, r, i, j; |
| if (n == 0 && emin == 0) { |
| return Series.zero(0, 0); |
| } else { |
| r = Series.zero(n, emin - 1); |
| for(i = 0; i < n; i++) { |
| j = emin + i; |
| if (j == 0) |
| r[i] = 0; |
| else |
| r[i] = j * a[i]; |
| } |
| return r.trim(); |
| } |
| }, |
| integ() { |
| var a = this, n = a.length, emin = a.emin, i, j, r; |
| r = Series.zero(n, emin + 1); |
| for(i = 0; i < n; i++) { |
| j = emin + i; |
| if (j == -1) { |
| if (a[i] != 0) |
| throw RangeError("cannot represent integ(1/X)"); |
| } else { |
| r[i] = a[i] / (j + 1); |
| } |
| } |
| return r.trim(); |
| }, |
| exp() { |
| var c, i, r, n, a = this; |
| if (a.emin < 0) |
| throw RangeError("negative exponent in exp"); |
| n = a.emin + a.length; |
| if (a.emin > 0 || a[0] == 0) { |
| c = 1; |
| } else { |
| c = global.exp(a[0]); |
| a -= a[0]; |
| } |
| r = Series.zero(n, 0); |
| for(i = 0; i < n; i++) { |
| r[i] = c / fact(i); |
| } |
| return r.apply(a); |
| }, |
| log() { |
| var a = this, r; |
| if (a.emin != 0) |
| throw RangeError("log argument must have a non zero constant term"); |
| r = integ(deriv(a) / a); |
| /* add the constant term */ |
| r += global.log(a[0]); |
| return r; |
| }, |
| }); |
| |
| add_props(Series, { |
| /* new series of length n and first exponent emin */ |
| zero(n, emin) { |
| var r, i, obj; |
| |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = 0; |
| /* we return an array and store emin in its prototype */ |
| obj = Object.create(Series.prototype); |
| obj.emin = emin; |
| Object.setPrototypeOf(r, obj); |
| return r; |
| }, |
| O(a) { |
| function ErrorO() { |
| return TypeError("invalid O() argument"); |
| } |
| var n; |
| if (series_is_scalar_or_polynomial(a)) { |
| a = Polynomial(a); |
| n = a.deg(); |
| if (n < 0) |
| throw ErrorO(); |
| } else if (a instanceof RationalFunction) { |
| if (a.num.deg() != 0) |
| throw ErrorO(); |
| n = a.den.deg(); |
| if (n < 0) |
| throw ErrorO(); |
| n = -n; |
| } else |
| throw ErrorO(); |
| return Series.zero(0, n); |
| }, |
| }); |
| |
| /* Array (Matrix) */ |
| |
| Matrix = function Matrix(h, w) { |
| var i, j, r, rl; |
| if (typeof w === "undefined") |
| w = h; |
| r = []; |
| for(i = 0; i < h; i++) { |
| rl = []; |
| for(j = 0; j < w; j++) |
| rl[j] = 0; |
| r[i] = rl; |
| } |
| return r; |
| }; |
| |
| add_props(Matrix, { |
| idn(n) { |
| var r, i; |
| r = Matrix(n, n); |
| for(i = 0; i < n; i++) |
| r[i][i] = 1; |
| return r; |
| }, |
| diag(a) { |
| var r, i, n; |
| n = a.length; |
| r = Matrix(n, n); |
| for(i = 0; i < n; i++) |
| r[i][i] = a[i]; |
| return r; |
| }, |
| hilbert(n) { |
| var i, j, r; |
| r = Matrix(n); |
| for(i = 0; i < n; i++) { |
| for(j = 0; j < n; j++) { |
| r[i][j] = 1 / (1 + i + j); |
| } |
| } |
| return r; |
| }, |
| trans(a) { |
| var h, w, r, i, j; |
| if (!Array.isArray(a)) |
| throw TypeError("matrix expected"); |
| h = a.length; |
| if (!Array.isArray(a[0])) { |
| w = 1; |
| r = Matrix(w, h); |
| for(i = 0; i < h; i++) { |
| r[0][i] = a[i]; |
| } |
| } else { |
| w = a[0].length; |
| r = Matrix(w, h); |
| for(i = 0; i < h; i++) { |
| for(j = 0; j < w; j++) { |
| r[j][i] = a[i][j]; |
| } |
| } |
| } |
| return r; |
| }, |
| check_square(a) { |
| var a, n; |
| if (!Array.isArray(a)) |
| throw TypeError("array expected"); |
| n = a.length; |
| if (!Array.isArray(a[0]) || n != a[0].length) |
| throw TypeError("square matrix expected"); |
| return n; |
| }, |
| trace(a) { |
| var n, r, i; |
| n = Matrix.check_square(a); |
| r = a[0][0]; |
| for(i = 1; i < n; i++) { |
| r += a[i][i]; |
| } |
| return r; |
| }, |
| charpoly(a) { |
| var n, p, c, i, j, coef; |
| n = Matrix.check_square(a); |
| p = []; |
| for(i = 0; i < n + 1; i++) |
| p[i] = 0; |
| p[n] = 1; |
| c = Matrix.idn(n); |
| for(i = 0; i < n; i++) { |
| c = c * a; |
| coef = -trace(c) / (i + 1); |
| p[n - i - 1] = coef; |
| for(j = 0; j < n; j++) |
| c[j][j] += coef; |
| } |
| return Polynomial(p); |
| }, |
| eigenvals(a) { |
| return Polynomial.roots(Matrix.charpoly(a)); |
| }, |
| det(a) { |
| var n, i, j, k, s, src, v, c; |
| |
| n = Matrix.check_square(a); |
| s = 1; |
| src = a.dup(); |
| for(i=0;i<n;i++) { |
| for(j = i; j < n; j++) { |
| if (src[j][i] != 0) |
| break; |
| } |
| if (j == n) |
| return 0; |
| if (j != i) { |
| for(k = 0;k < n; k++) { |
| v = src[j][k]; |
| src[j][k] = src[i][k]; |
| src[i][k] = v; |
| } |
| s = -s; |
| } |
| c = src[i][i].inverse(); |
| for(j = i + 1; j < n; j++) { |
| v = c * src[j][i]; |
| for(k = 0;k < n; k++) { |
| src[j][k] -= src[i][k] * v; |
| } |
| } |
| } |
| c = s; |
| for(i=0;i<n;i++) |
| c *= src[i][i]; |
| return c; |
| }, |
| inverse(a) { |
| var n, dst, src, i, j, k, n2, r, c, v; |
| n = Matrix.check_square(a); |
| src = a.dup(); |
| dst = Matrix.idn(n); |
| for(i=0;i<n;i++) { |
| for(j = i; j < n; j++) { |
| if (src[j][i] != 0) |
| break; |
| } |
| if (j == n) |
| throw RangeError("matrix is not invertible"); |
| if (j != i) { |
| /* swap lines in src and dst */ |
| v = src[j]; |
| src[j] = src[i]; |
| src[i] = v; |
| v = dst[j]; |
| dst[j] = dst[i]; |
| dst[i] = v; |
| } |
| |
| c = src[i][i].inverse(); |
| for(k = 0; k < n; k++) { |
| src[i][k] *= c; |
| dst[i][k] *= c; |
| } |
| |
| for(j = 0; j < n; j++) { |
| if (j != i) { |
| c = src[j][i]; |
| for(k = i; k < n; k++) { |
| src[j][k] -= src[i][k] * c; |
| } |
| for(k = 0; k < n; k++) { |
| dst[j][k] -= dst[i][k] * c; |
| } |
| } |
| } |
| } |
| return dst; |
| }, |
| rank(a) { |
| var src, i, j, k, w, h, l, c; |
| |
| if (!Array.isArray(a) || |
| !Array.isArray(a[0])) |
| throw TypeError("matrix expected"); |
| h = a.length; |
| w = a[0].length; |
| src = a.dup(); |
| l = 0; |
| for(i=0;i<w;i++) { |
| for(j = l; j < h; j++) { |
| if (src[j][i] != 0) |
| break; |
| } |
| if (j == h) |
| continue; |
| if (j != l) { |
| /* swap lines */ |
| for(k = 0; k < w; k++) { |
| v = src[j][k]; |
| src[j][k] = src[l][k]; |
| src[l][k] = v; |
| } |
| } |
| |
| c = src[l][i].inverse(); |
| for(k = 0; k < w; k++) { |
| src[l][k] *= c; |
| } |
| |
| for(j = l + 1; j < h; j++) { |
| c = src[j][i]; |
| for(k = i; k < w; k++) { |
| src[j][k] -= src[l][k] * c; |
| } |
| } |
| l++; |
| } |
| return l; |
| }, |
| ker(a) { |
| var src, i, j, k, w, h, l, m, r, im_cols, ker_dim, c; |
| |
| if (!Array.isArray(a) || |
| !Array.isArray(a[0])) |
| throw TypeError("matrix expected"); |
| h = a.length; |
| w = a[0].length; |
| src = a.dup(); |
| im_cols = []; |
| l = 0; |
| for(i=0;i<w;i++) { |
| im_cols[i] = false; |
| for(j = l; j < h; j++) { |
| if (src[j][i] != 0) |
| break; |
| } |
| if (j == h) |
| continue; |
| im_cols[i] = true; |
| if (j != l) { |
| /* swap lines */ |
| for(k = 0; k < w; k++) { |
| v = src[j][k]; |
| src[j][k] = src[l][k]; |
| src[l][k] = v; |
| } |
| } |
| |
| c = src[l][i].inverse(); |
| for(k = 0; k < w; k++) { |
| src[l][k] *= c; |
| } |
| |
| for(j = 0; j < h; j++) { |
| if (j != l) { |
| c = src[j][i]; |
| for(k = i; k < w; k++) { |
| src[j][k] -= src[l][k] * c; |
| } |
| } |
| } |
| l++; |
| // log_str("m=" + cval_toString(v1) + "\n"); |
| } |
| // log_str("im cols="+im_cols+"\n"); |
| |
| /* build the kernel vectors */ |
| ker_dim = w - l; |
| r = Matrix(w, ker_dim); |
| k = 0; |
| for(i = 0; i < w; i++) { |
| if (!im_cols[i]) { |
| /* select this column from the matrix */ |
| l = 0; |
| m = 0; |
| for(j = 0; j < w; j++) { |
| if (im_cols[j]) { |
| r[j][k] = -src[m][i]; |
| m++; |
| } else { |
| if (l == k) { |
| r[j][k] = 1; |
| } else { |
| r[j][k] = 0; |
| } |
| l++; |
| } |
| } |
| k++; |
| } |
| } |
| return r; |
| }, |
| dp(a, b) { |
| var i, n, r; |
| n = a.length; |
| if (n != b.length) |
| throw TypeError("incompatible array length"); |
| /* XXX: could do complex product */ |
| r = 0; |
| for(i = 0; i < n; i++) { |
| r += a[i] * b[i]; |
| } |
| return r; |
| }, |
| /* cross product */ |
| cp(v1, v2) { |
| var r; |
| if (v1.length != 3 || v2.length != 3) |
| throw TypeError("vectors must have 3 elements"); |
| r = []; |
| r[0] = v1[1] * v2[2] - v1[2] * v2[1]; |
| r[1] = v1[2] * v2[0] - v1[0] * v2[2]; |
| r[2] = v1[0] * v2[1] - v1[1] * v2[0]; |
| return r; |
| }, |
| }); |
| |
| function array_add(a, b) { |
| var r, i, n; |
| n = a.length; |
| if (n != b.length) |
| throw TypeError("incompatible array size"); |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = a[i] + b[i]; |
| return r; |
| } |
| function array_sub(a, b) { |
| var r, i, n; |
| n = a.length; |
| if (n != b.length) |
| throw TypeError("incompatible array size"); |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = a[i] - b[i]; |
| return r; |
| } |
| function array_scalar_mul(a, b) { |
| var r, i, n; |
| n = a.length; |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = a[i] * b; |
| return r; |
| } |
| function array_mul(a, b) { |
| var h, w, l, i, j, k, r, rl, sum, a_mat, b_mat; |
| h = a.length; |
| a_mat = Array.isArray(a[0]); |
| if (a_mat) { |
| l = a[0].length; |
| } else { |
| l = 1; |
| } |
| if (l != b.length) |
| throw RangeError("incompatible matrix size"); |
| b_mat = Array.isArray(b[0]); |
| if (b_mat) |
| w = b[0].length; |
| else |
| w = 1; |
| r = []; |
| if (a_mat && b_mat) { |
| for(i = 0; i < h; i++) { |
| rl = []; |
| for(j = 0; j < w; j++) { |
| sum = 0; |
| for(k = 0; k < l; k++) { |
| sum += a[i][k] * b[k][j]; |
| } |
| rl[j] = sum; |
| } |
| r[i] = rl; |
| } |
| } else if (a_mat && !b_mat) { |
| for(i = 0; i < h; i++) { |
| sum = 0; |
| for(k = 0; k < l; k++) { |
| sum += a[i][k] * b[k]; |
| } |
| r[i] = sum; |
| } |
| } else if (!a_mat && b_mat) { |
| for(i = 0; i < h; i++) { |
| rl = []; |
| for(j = 0; j < w; j++) { |
| rl[j] = a[i] * b[0][j]; |
| } |
| r[i] = rl; |
| } |
| } else { |
| for(i = 0; i < h; i++) { |
| r[i] = a[i] * b[0]; |
| } |
| } |
| return r; |
| } |
| function array_div(a, b) { |
| return array_mul(a, b.inverse()); |
| } |
| function array_scalar_div(a, b) { |
| return a * b.inverse(); |
| } |
| function array_eq(a, b) { |
| var n, i; |
| n = a.length; |
| if (n != b.length) |
| return false; |
| for(i = 0; i < n; i++) { |
| if (a[i] != b[i]) |
| return false; |
| } |
| return true; |
| } |
| |
| operators_set(Array.prototype, |
| { |
| "+": array_add, |
| "-": array_sub, |
| "*": array_mul, |
| "/": array_div, |
| "==": array_eq, |
| "pos"(a) { |
| return a; |
| }, |
| "neg"(a) { |
| var i, n, r; |
| n = a.length; |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = -a[i]; |
| return r; |
| } |
| }, |
| { |
| right: [Number, BigInt, Float, Fraction, Complex, Mod, |
| Polynomial, PolyMod, RationalFunction, Series], |
| "*": array_scalar_mul, |
| "/": array_scalar_div, |
| "**": generic_pow, /* XXX: only for integer */ |
| }, |
| { |
| left: [Number, BigInt, Float, Fraction, Complex, Mod, |
| Polynomial, PolyMod, RationalFunction, Series], |
| "*"(a, b) { return array_scalar_mul(b, a); }, |
| "/"(a, b) { return array_scalar_div(b, a); }, |
| }); |
| |
| add_props(Array.prototype, { |
| conj() { |
| var i, n, r; |
| n = this.length; |
| r = []; |
| for(i = 0; i < n; i++) |
| r[i] = this[i].conj(); |
| return r; |
| }, |
| dup() { |
| var r, i, n, el, a = this; |
| r = []; |
| n = a.length; |
| for(i = 0; i < n; i++) { |
| el = a[i]; |
| if (Array.isArray(el)) |
| el = el.dup(); |
| r[i] = el; |
| } |
| return r; |
| }, |
| inverse() { |
| return Matrix.inverse(this); |
| }, |
| norm2: Polynomial.prototype.norm2, |
| }); |
| |
| })(this); |
| |
| /* global definitions */ |
| var I = Complex(0, 1); |
| var X = Polynomial([0, 1]); |
| var O = Series.O; |
| |
| Object.defineProperty(this, "PI", { get: function () { return Float.PI } }); |
| |
| /* put frequently used functions in the global context */ |
| var gcd = Integer.gcd; |
| var fact = Integer.fact; |
| var comb = Integer.comb; |
| var pmod = Integer.pmod; |
| var invmod = Integer.invmod; |
| var factor = Integer.factor; |
| var isprime = Integer.isPrime; |
| var nextprime = Integer.nextPrime; |
| |
| function deriv(a) |
| { |
| return a.deriv(); |
| } |
| |
| function integ(a) |
| { |
| return a.integ(); |
| } |
| |
| function norm2(a) |
| { |
| return a.norm2(); |
| } |
| |
| function abs(a) |
| { |
| return a.abs(); |
| } |
| |
| function conj(a) |
| { |
| return a.conj(); |
| } |
| |
| function arg(a) |
| { |
| return a.arg(); |
| } |
| |
| function inverse(a) |
| { |
| return a.inverse(); |
| } |
| |
| function trunc(a) |
| { |
| if (Integer.isInteger(a)) { |
| return a; |
| } else if (a instanceof Fraction) { |
| return Integer.tdiv(a.num, a.den); |
| } else if (a instanceof Polynomial) { |
| return a; |
| } else if (a instanceof RationalFunction) { |
| return Polynomial.divrem(a.num, a.den)[0]; |
| } else { |
| return Float.ceil(a); |
| } |
| } |
| |
| function floor(a) |
| { |
| if (Integer.isInteger(a)) { |
| return a; |
| } else if (a instanceof Fraction) { |
| return Integer.fdiv(a.num, a.den); |
| } else { |
| return Float.floor(a); |
| } |
| } |
| |
| function ceil(a) |
| { |
| if (Integer.isInteger(a)) { |
| return a; |
| } else if (a instanceof Fraction) { |
| return Integer.cdiv(a.num, a.den); |
| } else { |
| return Float.ceil(a); |
| } |
| } |
| |
| function sqrt(a) |
| { |
| var t, u, re, im; |
| if (a instanceof Series) { |
| return a ^ (1/2); |
| } else if (a instanceof Complex) { |
| t = abs(a); |
| u = a.re; |
| re = sqrt((t + u) / 2); |
| im = sqrt((t - u) / 2); |
| if (a.im < 0) |
| im = -im; |
| return Complex.toComplex(re, im); |
| } else { |
| a = Float(a); |
| if (a < 0) { |
| return Complex(0, Float.sqrt(-a)); |
| } else { |
| return Float.sqrt(a); |
| } |
| } |
| } |
| |
| function exp(a) |
| { |
| return a.exp(); |
| } |
| |
| function log(a) |
| { |
| return a.log(); |
| } |
| |
| function log2(a) |
| { |
| return log(a) * Float.LOG2E; |
| } |
| |
| function log10(a) |
| { |
| return log(a) * Float.LOG10E; |
| } |
| |
| function todb(a) |
| { |
| return log10(a) * 10; |
| } |
| |
| function fromdb(a) |
| { |
| return 10 ^ (a / 10); |
| } |
| |
| function sin(a) |
| { |
| var t; |
| if (a instanceof Complex || a instanceof Series) { |
| t = exp(a * I); |
| return (t - 1/t) / (2 * I); |
| } else { |
| return Float.sin(Float(a)); |
| } |
| } |
| |
| function cos(a) |
| { |
| var t; |
| if (a instanceof Complex || a instanceof Series) { |
| t = exp(a * I); |
| return (t + 1/t) / 2; |
| } else { |
| return Float.cos(Float(a)); |
| } |
| } |
| |
| function tan(a) |
| { |
| if (a instanceof Complex || a instanceof Series) { |
| return sin(a) / cos(a); |
| } else { |
| return Float.tan(Float(a)); |
| } |
| } |
| |
| function asin(a) |
| { |
| return Float.asin(Float(a)); |
| } |
| |
| function acos(a) |
| { |
| return Float.acos(Float(a)); |
| } |
| |
| function atan(a) |
| { |
| return Float.atan(Float(a)); |
| } |
| |
| function atan2(a, b) |
| { |
| return Float.atan2(Float(a), Float(b)); |
| } |
| |
| function sinc(a) |
| { |
| if (a == 0) { |
| return 1; |
| } else { |
| a *= Float.PI; |
| return sin(a) / a; |
| } |
| } |
| |
| function todeg(a) |
| { |
| return a * 180 / Float.PI; |
| } |
| |
| function fromdeg(a) |
| { |
| return a * Float.PI / 180; |
| } |
| |
| function sinh(a) |
| { |
| var e = Float.exp(Float(a)); |
| return (e - 1/e) * 0.5; |
| } |
| |
| function cosh(a) |
| { |
| var e = Float.exp(Float(a)); |
| return (e + 1/e) * 0.5; |
| } |
| |
| function tanh(a) |
| { |
| var e = Float.exp(Float(a) * 2); |
| return (e - 1) / (e + 1); |
| } |
| |
| function asinh(a) |
| { |
| var x = Float(a); |
| return log(sqrt(x * x + 1) + x); |
| } |
| |
| function acosh(a) |
| { |
| var x = Float(a); |
| return log(sqrt(x * x - 1) + x); |
| } |
| |
| function atanh(a) |
| { |
| var x = Float(a); |
| return 0.5 * log((1 + x) / (1 - x)); |
| } |
| |
| var idn = Matrix.idn; |
| var diag = Matrix.diag; |
| var trans = Matrix.trans; |
| var trace = Matrix.trace; |
| var charpoly = Matrix.charpoly; |
| var eigenvals = Matrix.eigenvals; |
| var det = Matrix.det; |
| var rank = Matrix.rank; |
| var ker = Matrix.ker; |
| var cp = Matrix.cp; |
| var dp = Matrix.dp; |
| |
| var polroots = Polynomial.roots; |
| var bestappr = Float.bestappr; |
