quickjs-tart

quickjs-based runtime for wallet-core logic
Log | Files | Refs | README | LICENSE

qjscalc.js (72792B)

      1 /*
      2  * QuickJS Javascript Calculator
      3  *
      4  * Copyright (c) 2017-2020 Fabrice Bellard
      5  *
      6  * Permission is hereby granted, free of charge, to any person obtaining a copy
      7  * of this software and associated documentation files (the "Software"), to deal
      8  * in the Software without restriction, including without limitation the rights
      9  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
     10  * copies of the Software, and to permit persons to whom the Software is
     11  * furnished to do so, subject to the following conditions:
     12  *
     13  * The above copyright notice and this permission notice shall be included in
     14  * all copies or substantial portions of the Software.
     15  *
     16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
     17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
     18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
     19  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
     20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
     21  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
     22  * THE SOFTWARE.
     23  */
     24 "use strict";
     25 "use math";
     26 
     27 var Integer, Float, Fraction, Complex, Mod, Polynomial, PolyMod, RationalFunction, Series, Matrix;
     28 
     29 (function(global) {
     30     global.Integer = global.BigInt;
     31     global.Float = global.BigFloat;
     32     global.algebraicMode = true;
     33 
     34     /* add non enumerable properties */
     35     function add_props(obj, props) {
     36         var i, val, prop, tab, desc;
     37         tab = Reflect.ownKeys(props);
     38         for(i = 0; i < tab.length; i++) {
     39             prop = tab[i];
     40             desc = Object.getOwnPropertyDescriptor(props, prop);
     41             desc.enumerable = false;
     42             if ("value" in desc) {
     43                 if (typeof desc.value !== "function") {
     44                     desc.writable = false;
     45                     desc.configurable = false;
     46                 }
     47             } else {
     48                 /* getter/setter */
     49                 desc.configurable = false;
     50             }
     51             Object.defineProperty(obj, prop, desc);
     52         }
     53     }
     54 
     55     /* same as proto[Symbol.operatorSet] = Operators.create(..op_list)
     56        but allow shortcuts: left: [], right: [] or both
     57     */
     58     function operators_set(proto, ...op_list)
     59     {
     60         var new_op_list, i, a, j, b, k, obj, tab;
     61         var fields = [ "left", "right" ];
     62         new_op_list = [];
     63         for(i = 0; i < op_list.length; i++) {
     64             a = op_list[i];
     65             if (a.left || a.right) {
     66                 tab = [ a.left, a.right ];
     67                 delete a.left;
     68                 delete a.right;
     69                 for(k = 0; k < 2; k++) {
     70                     obj = tab[k];
     71                     if (obj) {
     72                         if (!Array.isArray(obj)) {
     73                             obj = [ obj ];
     74                         }
     75                         for(j = 0; j < obj.length; j++) {
     76                             b = {};
     77                             Object.assign(b, a);
     78                             b[fields[k]] = obj[j];
     79                             new_op_list.push(b);
     80                         }
     81                     }
     82                 }
     83             } else {
     84                 new_op_list.push(a);
     85             }
     86         }
     87         proto[Symbol.operatorSet] =
     88             Operators.create.call(null, ...new_op_list);
     89     }
     90 
     91     /* Integer */
     92 
     93     function generic_pow(a, b) {
     94         var r, is_neg, i;
     95         if (!Integer.isInteger(b)) {
     96             return exp(log(a) * b);
     97         }
     98         if (Array.isArray(a) && !(a instanceof Polynomial ||
     99                                   a instanceof Series)) {
    100             r = idn(Matrix.check_square(a));
    101         } else {
    102             r = 1;
    103         }
    104         if (b == 0)
    105             return r;
    106         is_neg = false;
    107         if (b < 0) {
    108             is_neg = true;
    109             b = -b;
    110         }
    111         r = a;
    112         for(i = Integer.floorLog2(b) - 1; i >= 0; i--) {
    113             r *= r;
    114             if ((b >> i) & 1)
    115                 r *= a;
    116         }
    117         if (is_neg) {
    118             if (typeof r.inverse != "function")
    119                 throw "negative powers are not supported for this type";
    120             r = r.inverse();
    121         }
    122         return r;
    123     }
    124 
    125     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 ];
    126 
    127     function miller_rabin_test(n, t) {
    128         var d, r, s, i, j, a;
    129         d = n - 1;
    130         s = 0;
    131         while ((d & 1) == 0) {
    132             d >>= 1;
    133             s++;
    134         }
    135         if (small_primes.length < t)
    136             t = small_primes.length;
    137         loop: for(j = 0; j < t; j++) {
    138             a = small_primes[j];
    139             r = Integer.pmod(a, d, n);
    140             if (r == 1 || r == (n - 1))
    141                 continue;
    142             for(i = 1; i < s; i++) {
    143                 r = (r * r) % n;
    144                 if (r == 1)
    145                     return false;
    146                 if (r == (n - 1))
    147                     continue loop;
    148             }
    149             return false; /* n is composite */
    150         }
    151         return true; /* n is probably prime with probability (1-0.5^t) */
    152     }
    153 
    154     function fact_rec(a, b) {  /* assumes a <= b */
    155         var i, r;
    156         if ((b - a) <= 5) {
    157             r = a;
    158             for(i = a + 1; i <= b; i++)
    159                 r *= i;
    160             return r;
    161         } else {
    162             /* to avoid a quadratic running time it is better to
    163                multiply numbers of similar size */
    164             i = (a + b) >> 1;
    165             return fact_rec(a, i) * fact_rec(i + 1, b);
    166         }
    167     }
    168 
    169     /* math mode specific quirk to overload the integer division and power */
    170     Operators.updateBigIntOperators(
    171         {
    172             "/"(a, b) {
    173                 if (algebraicMode) {
    174                     return Fraction.toFraction(a, b);
    175                 } else {
    176                     return Float(a) / Float(b);
    177                 }
    178             },
    179             "**"(a, b) {
    180                 if (algebraicMode) {
    181                     return generic_pow(a, b);
    182                 } else {
    183                     return Float(a) ** Float(b);
    184                 }
    185             }
    186         });
    187 
    188     add_props(Integer, {
    189         isInteger(a) {
    190             /* integers are represented either as bigint or as number */
    191             return typeof a === "bigint" ||
    192                 (typeof a === "number" && Number.isSafeInteger(a));
    193         },
    194         gcd(a, b) {
    195             var r;
    196             while (b != 0) {
    197                 r = a % b;
    198                 a = b;
    199                 b = r;
    200             }
    201             return a;
    202         },
    203         fact(n) {
    204             return n <= 0 ? 1 : fact_rec(1, n);
    205         },
    206         /* binomial coefficient */
    207         comb(n, k) {
    208             if (k < 0 || k > n)
    209                 return 0;
    210             if (k > n - k)
    211                 k = n - k;
    212             if (k == 0)
    213                 return 1;
    214             return Integer.tdiv(fact_rec(n - k + 1, n), fact_rec(1, k));
    215         },
    216         /* inverse of x modulo y */
    217         invmod(x, y) {
    218             var q, u, v, a, c, t;
    219             u = x;
    220             v = y;
    221             c = 1;
    222             a = 0;
    223             while (u != 0) {
    224                 t = Integer.fdivrem(v, u);
    225                 q = t[0];
    226                 v = u;
    227                 u = t[1];
    228                 t = c;
    229                 c = a - q * c;
    230                 a = t;
    231             }
    232             /* v = gcd(x, y) */
    233             if (v != 1)
    234                 throw RangeError("not invertible");
    235             return a % y;
    236         },
    237         /* return a ^ b modulo m */
    238         pmod(a, b, m) {
    239             var r;
    240             if (b == 0)
    241                 return 1;
    242             if (b < 0) {
    243                 a = Integer.invmod(a, m);
    244                 b = -b;
    245             }
    246             r = 1;
    247             for(;;) {
    248                 if (b & 1) {
    249                     r = (r * a) % m;
    250                 }
    251                 b >>= 1;
    252                 if (b == 0)
    253                     break;
    254                 a = (a * a) % m;
    255             }
    256             return r;
    257         },
    258 
    259         /* return true if n is prime (or probably prime with
    260            probability 1-0.5^t) */
    261         isPrime(n, t) {
    262             var i, d, n1;
    263             if (!Integer.isInteger(n))
    264                 throw TypeError("invalid type");
    265             if (n <= 1)
    266                 return false;
    267             n1 = small_primes.length;
    268             /* XXX: need Integer.sqrt() */
    269             for(i = 0; i < n1; i++) {
    270                 d = small_primes[i];
    271                 if (d == n)
    272                     return true;
    273                 if (d > n)
    274                     return false;
    275                 if ((n % d) == 0)
    276                     return false;
    277             }
    278             if (n < d * d)
    279                 return true;
    280             if (typeof t == "undefined")
    281                 t = 64;
    282             return miller_rabin_test(n, t);
    283         },
    284         nextPrime(n) {
    285             if (!Integer.isInteger(n))
    286                 throw TypeError("invalid type");
    287             if (n < 1)
    288                 n = 1;
    289             for(;;) {
    290                 n++;
    291                 if (Integer.isPrime(n))
    292                     return n;
    293             }
    294         },
    295         factor(n) {
    296             var r, d;
    297             if (!Integer.isInteger(n))
    298                 throw TypeError("invalid type");
    299             r = [];
    300             if (abs(n) <= 1) {
    301                 r.push(n);
    302                 return r;
    303             }
    304             if (n < 0) {
    305                 r.push(-1);
    306                 n = -n;
    307             }
    308 
    309             while ((n % 2) == 0) {
    310                 n >>= 1;
    311                 r.push(2);
    312             }
    313 
    314             d = 3;
    315             while (n != 1) {
    316                 if (Integer.isPrime(n)) {
    317                     r.push(n);
    318                     break;
    319                 }
    320                 /* we are sure there is at least one divisor, so one test */
    321                 for(;;) {
    322                     if ((n % d) == 0)
    323                         break;
    324                     d += 2;
    325                 }
    326                 for(;;) {
    327                     r.push(d);
    328                     n = Integer.tdiv(n, d);
    329                     if ((n % d) != 0)
    330                         break;
    331                 }
    332             }
    333             return r;
    334         },
    335     });
    336 
    337     add_props(Integer.prototype, {
    338         inverse() {
    339             return 1 / this;
    340         },
    341         norm2() {
    342             return this * this;
    343         },
    344         abs() {
    345             var v = this;
    346             if (v < 0)
    347                 v = -v;
    348             return v;
    349         },
    350         conj() {
    351             return this;
    352         },
    353         arg() {
    354             if (this >= 0)
    355                 return 0;
    356             else
    357                 return Float.PI;
    358         },
    359         exp() {
    360             if (this == 0)
    361                 return 1;
    362             else
    363                 return Float.exp(this);
    364         },
    365         log() {
    366             if (this == 1)
    367                 return 0;
    368             else
    369                 return Float(this).log();
    370         },
    371     });
    372 
    373     /* Fraction */
    374 
    375     Fraction = function Fraction(a, b)
    376     {
    377         var d, r, obj;
    378 
    379         if (new.target)
    380             throw TypeError("not a constructor");
    381         if (a instanceof Fraction)
    382             return a;
    383         if (!Integer.isInteger(a))
    384             throw TypeError("integer expected");
    385         if (typeof b === "undefined") {
    386             b = 1;
    387         } else {
    388             if (!Integer.isInteger(b))
    389                 throw TypeError("integer expected");
    390             if (b == 0)
    391                 throw RangeError("division by zero");
    392             d = Integer.gcd(a, b);
    393             if (d != 1) {
    394                 a = Integer.tdiv(a, d);
    395                 b = Integer.tdiv(b, d);
    396             }
    397 
    398             /* the fractions are normalized with den > 0 */
    399             if (b < 0) {
    400                 a = -a;
    401                 b = -b;
    402             }
    403         }
    404         obj = Object.create(Fraction.prototype);
    405         obj.num = a;
    406         obj.den = b;
    407         return obj;
    408     }
    409 
    410     function fraction_add(a, b) {
    411         a = Fraction(a);
    412         b = Fraction(b);
    413         return Fraction.toFraction(a.num * b.den + a.den * b.num, a.den * b.den);
    414     }
    415     function fraction_sub(a, b) {
    416         a = Fraction(a);
    417         b = Fraction(b);
    418         return Fraction.toFraction(a.num * b.den - a.den * b.num, a.den * b.den);
    419     }
    420     function fraction_mul(a, b) {
    421         a = Fraction(a);
    422         b = Fraction(b);
    423         return Fraction.toFraction(a.num * b.num, a.den * b.den);
    424     }
    425     function fraction_div(a, b) {
    426         a = Fraction(a);
    427         b = Fraction(b);
    428         return Fraction.toFraction(a.num * b.den, a.den * b.num);
    429     }
    430     function fraction_mod(a, b) {
    431         var a1 = Fraction(a);
    432         var b1 = Fraction(b);
    433         return a - Integer.ediv(a1.num * b1.den, a1.den * b1.num) * b;
    434     }
    435     function fraction_eq(a, b) {
    436         a = Fraction(a);
    437         b = Fraction(b);
    438         /* we assume the fractions are normalized */
    439         return (a.num == b.num && a.den == b.den);
    440     }
    441     function fraction_lt(a, b) {
    442         a = Fraction(a);
    443         b = Fraction(b);
    444         return (a.num * b.den < b.num * a.den);
    445     }
    446 
    447     /* operators are needed for fractions */
    448     function float_add(a, b) {
    449         return Float(a) + Float(b);
    450     }
    451     function float_sub(a, b) {
    452         return Float(a) - Float(b);
    453     }
    454     function float_mul(a, b) {
    455         return Float(a) * Float(b);
    456     }
    457     function float_div(a, b) {
    458         return Float(a) / Float(b);
    459     }
    460     function float_mod(a, b) {
    461         return Float(a) % Float(b);
    462     }
    463     function float_pow(a, b) {
    464         return Float(a) ** Float(b);
    465     }
    466     function float_eq(a, b) {
    467         /* XXX: may be better to use infinite precision for the comparison */
    468         return Float(a) === Float(b);
    469     }
    470     function float_lt(a, b) {
    471         a = Float(a);
    472         b = Float(b);
    473         /* XXX: may be better to use infinite precision for the comparison */
    474         if (Float.isNaN(a) || Float.isNaN(b))
    475             return undefined;
    476         else
    477             return a < b;
    478     }
    479 
    480     operators_set(Fraction.prototype,
    481         {
    482             "+": fraction_add,
    483             "-": fraction_sub,
    484             "*": fraction_mul,
    485             "/": fraction_div,
    486             "%": fraction_mod,
    487             "**": generic_pow,
    488             "==": fraction_eq,
    489             "<": fraction_lt,
    490             "pos"(a) {
    491                 return a;
    492             },
    493             "neg"(a) {
    494                 return Fraction(-a.num, a.den);
    495             },
    496         },
    497         {
    498             left: [Number, BigInt],
    499             right: [Number, BigInt],
    500             "+": fraction_add,
    501             "-": fraction_sub,
    502             "*": fraction_mul,
    503             "/": fraction_div,
    504             "%": fraction_mod,
    505             "**": generic_pow,
    506             "==": fraction_eq,
    507             "<": fraction_lt,
    508         },
    509         {
    510             left: Float,
    511             right: Float,
    512             "+": float_add,
    513             "-": float_sub,
    514             "*": float_mul,
    515             "/": float_div,
    516             "%": float_mod,
    517             "**": float_pow,
    518             "==": float_eq,
    519             "<": float_lt,
    520         });
    521 
    522     add_props(Fraction, {
    523         /* (internal use) simplify 'a' to an integer when possible */
    524         toFraction(a, b) {
    525             var r = Fraction(a, b);
    526             if (algebraicMode && r.den == 1)
    527                 return r.num;
    528             else
    529                 return r;
    530         },
    531     });
    532 
    533     add_props(Fraction.prototype, {
    534         [Symbol.toPrimitive](hint) {
    535             if (hint === "string") {
    536                 return this.toString();
    537             } else {
    538                 return Float(this.num) / this.den;
    539             }
    540         },
    541         inverse() {
    542             return Fraction(this.den, this.num);
    543         },
    544         toString() {
    545             return this.num + "/" + this.den;
    546         },
    547         norm2() {
    548             return this * this;
    549         },
    550         abs() {
    551             if (this.num < 0)
    552                 return -this;
    553             else
    554                 return this;
    555         },
    556         conj() {
    557             return this;
    558         },
    559         arg() {
    560             if (this.num >= 0)
    561                 return 0;
    562             else
    563                 return Float.PI;
    564         },
    565         exp() {
    566             return Float.exp(Float(this));
    567         },
    568         log() {
    569             return Float(this).log();
    570         },
    571     });
    572 
    573     /* Number (Float64) */
    574 
    575     add_props(Number.prototype, {
    576         inverse() {
    577             return 1 / this;
    578         },
    579         norm2() {
    580             return this * this;
    581         },
    582         abs() {
    583             return Math.abs(this);
    584         },
    585         conj() {
    586             return this;
    587         },
    588         arg() {
    589             if (this >= 0)
    590                 return 0;
    591             else
    592                 return Float.PI;
    593         },
    594         exp() {
    595             return Float.exp(this);
    596         },
    597         log() {
    598             if (this < 0) {
    599                 return Complex(this).log();
    600             } else {
    601                 return Float.log(this);
    602             }
    603         },
    604     });
    605 
    606     /* Float */
    607 
    608     var const_tab = [];
    609 
    610     /* we cache the constants for small precisions */
    611     function get_const(n) {
    612         var t, c, p;
    613         t = const_tab[n];
    614         p = BigFloatEnv.prec;
    615         if (t && t.prec == p) {
    616             return t.val;
    617         } else {
    618             switch(n) {
    619             case 0: c = Float.exp(1); break;
    620             case 1: c = Float.log(10); break;
    621 //            case 2: c = Float.log(2); break;
    622             case 3: c = 1/Float.log(2); break;
    623             case 4: c = 1/Float.log(10); break;
    624 //            case 5: c = Float.atan(1) * 4; break;
    625             case 6: c = Float.sqrt(0.5); break;
    626             case 7: c = Float.sqrt(2); break;
    627             }
    628             if (p <= 1024) {
    629                 const_tab[n] = { prec: p, val: c };
    630             }
    631             return c;
    632         }
    633     }
    634 
    635     add_props(Float, {
    636         isFloat(a) {
    637             return typeof a === "number" || typeof a === "bigfloat";
    638         },
    639         bestappr(u, b) {
    640             var num1, num0, den1, den0, u, num, den, n;
    641 
    642             if (typeof b === "undefined")
    643                 throw TypeError("second argument expected");
    644             num1 = 1;
    645             num0 = 0;
    646             den1 = 0;
    647             den0 = 1;
    648             for(;;) {
    649                 n = Integer(Float.floor(u));
    650                 num = n * num1 + num0;
    651                 den = n * den1 + den0;
    652                 if (den > b)
    653                     break;
    654                 u = 1.0 / (u - n);
    655                 num0 = num1;
    656                 num1 = num;
    657                 den0 = den1;
    658                 den1 = den;
    659             }
    660             return Fraction(num1, den1);
    661         },
    662         /* similar constants as Math.x */
    663         get E() { return get_const(0); },
    664         get LN10() { return get_const(1); },
    665 //        get LN2() { return get_const(2); },
    666         get LOG2E() { return get_const(3); },
    667         get LOG10E() { return get_const(4); },
    668 //        get PI() { return get_const(5); },
    669         get SQRT1_2() { return get_const(6); },
    670         get SQRT2() { return get_const(7); },
    671     });
    672 
    673     add_props(Float.prototype, {
    674         inverse() {
    675             return 1.0 / this;
    676         },
    677         norm2() {
    678             return this * this;
    679         },
    680         abs() {
    681             return Float.abs(this);
    682         },
    683         conj() {
    684             return this;
    685         },
    686         arg() {
    687             if (this >= 0)
    688                 return 0;
    689             else
    690                 return Float.PI;
    691         },
    692         exp() {
    693             return Float.exp(this);
    694         },
    695         log() {
    696             if (this < 0) {
    697                 return Complex(this).log();
    698             } else {
    699                 return Float.log(this);
    700             }
    701         },
    702     });
    703 
    704     /* Complex */
    705 
    706     Complex = function Complex(re, im)
    707     {
    708         var obj;
    709         if (new.target)
    710             throw TypeError("not a constructor");
    711         if (re instanceof Complex)
    712             return re;
    713         if (typeof im === "undefined") {
    714             im = 0;
    715         }
    716         obj = Object.create(Complex.prototype);
    717         obj.re = re;
    718         obj.im = im;
    719         return obj;
    720     }
    721 
    722 
    723     function complex_add(a, b) {
    724         a = Complex(a);
    725         b = Complex(b);
    726         return Complex.toComplex(a.re + b.re, a.im + b.im);
    727     }
    728     function complex_sub(a, b) {
    729         a = Complex(a);
    730         b = Complex(b);
    731         return Complex.toComplex(a.re - b.re, a.im - b.im);
    732     }
    733     function complex_mul(a, b) {
    734         a = Complex(a);
    735         b = Complex(b);
    736         return Complex.toComplex(a.re * b.re - a.im * b.im,
    737                                  a.re * b.im + a.im * b.re);
    738     }
    739     function complex_div(a, b) {
    740         a = Complex(a);
    741         b = Complex(b);
    742         return a * b.inverse();
    743     }
    744     function complex_eq(a, b) {
    745         a = Complex(a);
    746         b = Complex(b);
    747         return a.re == b.re && a.im == b.im;
    748     }
    749 
    750     operators_set(Complex.prototype,
    751         {
    752             "+": complex_add,
    753             "-": complex_sub,
    754             "*": complex_mul,
    755             "/": complex_div,
    756             "**": generic_pow,
    757             "==": complex_eq,
    758             "pos"(a) {
    759                 return a;
    760             },
    761             "neg"(a) {
    762                 return Complex(-a.re, -a.im);
    763             }
    764         },
    765         {
    766             left: [Number, BigInt, Float, Fraction],
    767             right: [Number, BigInt, Float, Fraction],
    768             "+": complex_add,
    769             "-": complex_sub,
    770             "*": complex_mul,
    771             "/": complex_div,
    772             "**": generic_pow,
    773             "==": complex_eq,
    774         });
    775 
    776     add_props(Complex, {
    777         /* simplify to real number when possible */
    778         toComplex(re, im) {
    779             if (algebraicMode && im == 0)
    780                 return re;
    781             else
    782                 return Complex(re, im);
    783         },
    784     });
    785 
    786     add_props(Complex.prototype, {
    787         inverse() {
    788             var c = this.norm2();
    789             return Complex(this.re / c, -this.im / c);
    790         },
    791         toString() {
    792             var v, s = "", a = this;
    793             if (a.re != 0)
    794                 s += a.re.toString();
    795             if (a.im == 1) {
    796                 if (s != "")
    797                     s += "+";
    798                 s += "I";
    799             } else if (a.im == -1) {
    800                 s += "-I";
    801             } else {
    802                 v = a.im.toString();
    803                 if (v[0] != "-" && s != "")
    804                     s += "+";
    805                 s += v + "*I";
    806             }
    807             return s;
    808         },
    809         norm2() {
    810             return this.re * this.re + this.im * this.im;
    811         },
    812         abs() {
    813             return Float.sqrt(norm2(this));
    814         },
    815         conj() {
    816             return Complex(this.re, -this.im);
    817         },
    818         arg() {
    819             return Float.atan2(this.im, this.re);
    820         },
    821         exp() {
    822             var arg = this.im, r = this.re.exp();
    823             return Complex(r * cos(arg), r * sin(arg));
    824         },
    825         log() {
    826             return Complex(abs(this).log(), atan2(this.im, this.re));
    827         },
    828     });
    829 
    830     /* Mod */
    831 
    832     Mod = function Mod(a, m) {
    833         var obj, t;
    834         if (new.target)
    835             throw TypeError("not a constructor");
    836         obj = Object.create(Mod.prototype);
    837         if (Integer.isInteger(m)) {
    838             if (m <= 0)
    839                 throw RangeError("the modulo cannot be <= 0");
    840             if (Integer.isInteger(a)) {
    841                 a %= m;
    842             } else if (a instanceof Fraction) {
    843                 return Mod(a.num, m) / a.den;
    844             } else {
    845                 throw TypeError("invalid types");
    846             }
    847         } else {
    848             throw TypeError("invalid types");
    849         }
    850         obj.res = a;
    851         obj.mod = m;
    852         return obj;
    853     };
    854 
    855     function mod_add(a, b) {
    856         if (!(a instanceof Mod)) {
    857             return Mod(a + b.res, b.mod);
    858         } else if (!(b instanceof Mod)) {
    859             return Mod(a.res + b, a.mod);
    860         } else {
    861             if (a.mod != b.mod)
    862                 throw TypeError("different modulo for binary operator");
    863             return Mod(a.res + b.res, a.mod);
    864         }
    865     }
    866     function mod_sub(a, b) {
    867         if (!(a instanceof Mod)) {
    868             return Mod(a - b.res, b.mod);
    869         } else if (!(b instanceof Mod)) {
    870             return Mod(a.res - b, a.mod);
    871         } else {
    872             if (a.mod != b.mod)
    873                 throw TypeError("different modulo for binary operator");
    874             return Mod(a.res - b.res, a.mod);
    875         }
    876     }
    877     function mod_mul(a, b) {
    878         if (!(a instanceof Mod)) {
    879             return Mod(a * b.res, b.mod);
    880         } else if (!(b instanceof Mod)) {
    881             return Mod(a.res * b, a.mod);
    882         } else {
    883             if (a.mod != b.mod)
    884                 throw TypeError("different modulo for binary operator");
    885             return Mod(a.res * b.res, a.mod);
    886         }
    887     }
    888     function mod_div(a, b) {
    889         if (!(b instanceof Mod))
    890             b = Mod(b, a.mod);
    891         return mod_mul(a, b.inverse());
    892     }
    893     function mod_eq(a, b) {
    894         return (a.mod == b.mod && a.res == b.res);
    895     }
    896 
    897     operators_set(Mod.prototype,
    898         {
    899             "+": mod_add,
    900             "-": mod_sub,
    901             "*": mod_mul,
    902             "/": mod_div,
    903             "**": generic_pow,
    904             "==": mod_eq,
    905             "pos"(a) {
    906                 return a;
    907             },
    908             "neg"(a) {
    909                 return Mod(-a.res, a.mod);
    910             }
    911         },
    912         {
    913             left: [Number, BigInt, Float, Fraction],
    914             right: [Number, BigInt, Float, Fraction],
    915             "+": mod_add,
    916             "-": mod_sub,
    917             "*": mod_mul,
    918             "/": mod_div,
    919             "**": generic_pow,
    920         });
    921 
    922     add_props(Mod.prototype, {
    923         inverse() {
    924             var a = this, m = a.mod;
    925             if (Integer.isInteger(m)) {
    926                 return Mod(Integer.invmod(a.res, m), m);
    927             } else {
    928                 throw TypeError("unsupported type");
    929             }
    930         },
    931         toString() {
    932             return "Mod(" + this.res + "," + this.mod + ")";
    933         },
    934     });
    935 
    936     /* Polynomial */
    937 
    938     function polynomial_is_scalar(a)
    939     {
    940         if (typeof a === "number" ||
    941             typeof a === "bigint" ||
    942             typeof a === "bigfloat")
    943             return true;
    944         if (a instanceof Fraction ||
    945             a instanceof Complex ||
    946             a instanceof Mod)
    947             return true;
    948         return false;
    949     }
    950 
    951     Polynomial = function Polynomial(a)
    952     {
    953         if (new.target)
    954             throw TypeError("not a constructor");
    955         if (a instanceof Polynomial) {
    956             return a;
    957         } else if (Array.isArray(a)) {
    958             if (a.length == 0)
    959                 a = [ 0 ];
    960             Object.setPrototypeOf(a, Polynomial.prototype);
    961             return a.trim();
    962         } else if (polynomial_is_scalar(a)) {
    963             a = [a];
    964             Object.setPrototypeOf(a, Polynomial.prototype);
    965             return a;
    966         } else {
    967             throw TypeError("invalid type");
    968         }
    969     }
    970 
    971     function number_need_paren(c)
    972     {
    973         return !(Integer.isInteger(c) ||
    974                  Float.isFloat(c) ||
    975                  c instanceof Fraction ||
    976                  (c instanceof Complex && c.re == 0));
    977     }
    978 
    979     /* string for c*X^i */
    980     function monomial_toString(c, i)
    981     {
    982         var str1;
    983         if (i == 0) {
    984             str1 = c.toString();
    985         } else {
    986             if (c == 1) {
    987                 str1 = "";
    988             } else if (c == -1) {
    989                 str1 = "-";
    990             } else {
    991                 if (number_need_paren(c)) {
    992                     str1 = "(" + c + ")";
    993                 } else {
    994                     str1 = String(c);
    995                 }
    996                 str1 += "*";
    997             }
    998             str1 += "X";
    999             if (i != 1) {
   1000                 str1 += "^" + i;
   1001             }
   1002         }
   1003         return str1;
   1004     }
   1005 
   1006     /* find one complex root of 'p' starting from z at precision eps using
   1007        at most max_it iterations. Return null if could not find root. */
   1008     function poly_root_laguerre1(p, z, max_it)
   1009     {
   1010         var p1, p2, i, z0, z1, z2, d, t0, t1, d1, d2, e, el, zl;
   1011 
   1012         d = p.deg();
   1013         if (d == 1) {
   1014             /* monomial case */
   1015             return -p[0] / p[1];
   1016         }
   1017         /* trivial zero */
   1018         if (p[0] == 0)
   1019             return 0.0;
   1020 
   1021         p1 = p.deriv();
   1022         p2 = p1.deriv();
   1023         el = 0.0;
   1024         zl = 0.0;
   1025         for(i = 0; i < max_it; i++) {
   1026             z0 = p.apply(z);
   1027             if (z0 == 0)
   1028                 return z; /* simple exit case */
   1029 
   1030             /* Ward stopping criteria */
   1031             e = abs(z - zl);
   1032 //            print("e", i, e);
   1033             if (i >= 2 && e >= el) {
   1034                 if (abs(zl) < 1e-4) {
   1035                     if (e < 1e-7)
   1036                         return zl;
   1037                 } else {
   1038                     if (e < abs(zl) * 1e-3)
   1039                         return zl;
   1040                 }
   1041             }
   1042             el = e;
   1043             zl = z;
   1044 
   1045             z1 = p1.apply(z);
   1046             z2 = p2.apply(z);
   1047             t0 = (d - 1) * z1;
   1048             t0 = t0 * t0;
   1049             t1 = d * (d - 1) * z0 * z2;
   1050             t0 = sqrt(t0 - t1);
   1051             d1 = z1 + t0;
   1052             d2 = z1 - t0;
   1053             if (norm2(d2) > norm2(d1))
   1054                 d1 = d2;
   1055             if (d1 == 0)
   1056                 return null;
   1057             z = z - d * z0 / d1;
   1058         }
   1059         return null;
   1060     }
   1061 
   1062     function poly_roots(p)
   1063     {
   1064         var d, i, roots, j, z, eps;
   1065         var start_points = [ 0.1, -1.4, 1.7 ];
   1066 
   1067         if (!(p instanceof Polynomial))
   1068             throw TypeError("polynomial expected");
   1069         d = p.deg();
   1070         if (d <= 0)
   1071             return [];
   1072         eps = 2.0 ^ (-BigFloatEnv.prec);
   1073         roots = [];
   1074         for(i = 0; i < d; i++) {
   1075             /* XXX: should select another start point if error */
   1076             for(j = 0; j < 3; j++) {
   1077                 z = poly_root_laguerre1(p, start_points[j], 100);
   1078                 if (z !== null)
   1079                     break;
   1080             }
   1081             if (j == 3)
   1082                 throw RangeError("error in root finding algorithm");
   1083             roots[i] = z;
   1084             p = Polynomial.divrem(p, X - z)[0];
   1085         }
   1086         return roots;
   1087     }
   1088 
   1089     add_props(Polynomial.prototype, {
   1090         trim() {
   1091             var a = this, i;
   1092             i = a.length;
   1093             while (i > 1 && a[i - 1] == 0)
   1094                 i--;
   1095             a.length = i;
   1096             return a;
   1097         },
   1098         conj() {
   1099             var r, i, n, a;
   1100             a = this;
   1101             n = a.length;
   1102             r = [];
   1103             for(i = 0; i < n; i++)
   1104                 r[i] = a[i].conj();
   1105             return Polynomial(r);
   1106         },
   1107         inverse() {
   1108             return RationalFunction(Polynomial([1]), this);
   1109         },
   1110         toString() {
   1111             var i, str, str1, c, a = this;
   1112             if (a.length == 1) {
   1113                 return a[0].toString();
   1114             }
   1115             str="";
   1116             for(i = a.length - 1; i >= 0; i--) {
   1117                 c = a[i];
   1118                 if (c == 0 ||
   1119                     (c instanceof Mod) && c.res == 0)
   1120                     continue;
   1121                 str1 = monomial_toString(c, i);
   1122                 if (str1[0] != "-") {
   1123                     if (str != "")
   1124                         str += "+";
   1125                 }
   1126                 str += str1;
   1127             }
   1128             return str;
   1129         },
   1130         deg() {
   1131             if (this.length == 1 && this[0] == 0)
   1132                 return -Infinity;
   1133             else
   1134                 return this.length - 1;
   1135         },
   1136         apply(b) {
   1137             var i, n, r, a = this;
   1138             n = a.length - 1;
   1139             r = a[n];
   1140             while (n > 0) {
   1141                 n--;
   1142                 r = r * b + a[n];
   1143             }
   1144             return r;
   1145         },
   1146         deriv() {
   1147             var a = this, n, r, i;
   1148             n = a.length;
   1149             if (n == 1) {
   1150                 return Polynomial(0);
   1151             } else {
   1152                 r = [];
   1153                 for(i = 1; i < n; i++) {
   1154                     r[i - 1] = i * a[i];
   1155                 }
   1156                 return Polynomial(r);
   1157             }
   1158         },
   1159         integ() {
   1160             var a = this, n, r, i;
   1161             n = a.length;
   1162             r = [0];
   1163             for(i = 0; i < n; i++) {
   1164                 r[i + 1] = a[i] / (i + 1);
   1165             }
   1166             return Polynomial(r);
   1167         },
   1168         norm2() {
   1169             var a = this, n, r, i;
   1170             n = a.length;
   1171             r = 0;
   1172             for(i = 0; i < n; i++) {
   1173                 r += a[i].norm2();
   1174             }
   1175             return r;
   1176         },
   1177     });
   1178 
   1179 
   1180     function polynomial_add(a, b) {
   1181         var tmp, r, i, n1, n2;
   1182         a = Polynomial(a);
   1183         b = Polynomial(b);
   1184         if (a.length < b.length) {
   1185             tmp = a;
   1186             a = b;
   1187             b = tmp;
   1188         }
   1189         n1 = b.length;
   1190         n2 = a.length;
   1191         r = [];
   1192         for(i = 0; i < n1; i++)
   1193             r[i] = a[i] + b[i];
   1194         for(i = n1; i < n2; i++)
   1195             r[i] = a[i];
   1196         return Polynomial(r);
   1197     }
   1198     function polynomial_sub(a, b) {
   1199         return polynomial_add(a, -b);
   1200     }
   1201     function polynomial_mul(a, b) {
   1202         var i, j, n1, n2, n, r;
   1203         a = Polynomial(a);
   1204         b = Polynomial(b);
   1205         n1 = a.length;
   1206         n2 = b.length;
   1207         n = n1 + n2 - 1;
   1208         r = [];
   1209         for(i = 0; i < n; i++)
   1210             r[i] = 0;
   1211         for(i = 0; i < n1; i++) {
   1212             for(j = 0; j < n2; j++) {
   1213                 r[i + j] += a[i] * b[j];
   1214             }
   1215         }
   1216         return Polynomial(r);
   1217     }
   1218     function polynomial_div_scalar(a, b) {
   1219         return a * (1 / b);
   1220     }
   1221     function polynomial_div(a, b)
   1222     {
   1223         return RationalFunction(Polynomial(a),
   1224                                 Polynomial(b));
   1225     }
   1226     function polynomial_mod(a, b) {
   1227         return Polynomial.divrem(a, b)[1];
   1228     }
   1229     function polynomial_eq(a, b) {
   1230         var n, i;
   1231         n = a.length;
   1232         if (n != b.length)
   1233             return false;
   1234         for(i = 0; i < n; i++) {
   1235             if (a[i] != b[i])
   1236                 return false;
   1237         }
   1238         return true;
   1239     }
   1240 
   1241     operators_set(Polynomial.prototype,
   1242         {
   1243             "+": polynomial_add,
   1244             "-": polynomial_sub,
   1245             "*": polynomial_mul,
   1246             "/": polynomial_div,
   1247             "**": generic_pow,
   1248             "==": polynomial_eq,
   1249             "pos"(a) {
   1250                 return a;
   1251             },
   1252             "neg"(a) {
   1253                 var r, i, n, a;
   1254                 n = a.length;
   1255                 r = [];
   1256                 for(i = 0; i < n; i++)
   1257                 r[i] = -a[i];
   1258                 return Polynomial(r);
   1259             },
   1260         },
   1261         {
   1262             left: [Number, BigInt, Float, Fraction, Complex, Mod],
   1263             "+": polynomial_add,
   1264             "-": polynomial_sub,
   1265             "*": polynomial_mul,
   1266             "/": polynomial_div,
   1267             "**": generic_pow, /* XXX: only for integer */
   1268         },
   1269         {
   1270             right: [Number, BigInt, Float, Fraction, Complex, Mod],
   1271             "+": polynomial_add,
   1272             "-": polynomial_sub,
   1273             "*": polynomial_mul,
   1274             "/": polynomial_div_scalar,
   1275             "**": generic_pow, /* XXX: only for integer */
   1276         });
   1277 
   1278     add_props(Polynomial, {
   1279         divrem(a, b) {
   1280             var n1, n2, i, j, q, r, n, c;
   1281             if (b.deg() < 0)
   1282                 throw RangeError("division by zero");
   1283             n1 = a.length;
   1284             n2 = b.length;
   1285             if (n1 < n2)
   1286                 return [Polynomial([0]), a];
   1287             r = Array.prototype.dup.call(a);
   1288             q = [];
   1289             n2--;
   1290             n = n1 - n2;
   1291             for(i = 0; i < n; i++)
   1292                 q[i] = 0;
   1293             for(i = n - 1; i >= 0; i--) {
   1294                 c = r[i + n2];
   1295                 if (c != 0) {
   1296                     c = c / b[n2];
   1297                     r[i + n2] = 0;
   1298                     for(j = 0; j < n2; j++) {
   1299                         r[i + j] -= b[j] * c;
   1300                     }
   1301                     q[i] = c;
   1302                 }
   1303             }
   1304             return [Polynomial(q), Polynomial(r)];
   1305         },
   1306         gcd(a, b) {
   1307             var t;
   1308             while (b.deg() >= 0) {
   1309                 t = Polynomial.divrem(a, b);
   1310                 a = b;
   1311                 b = t[1];
   1312             }
   1313             /* convert to monic form */
   1314             return a / a[a.length - 1];
   1315         },
   1316         invmod(x, y) {
   1317             var q, u, v, a, c, t;
   1318             u = x;
   1319             v = y;
   1320             c = Polynomial([1]);
   1321             a = Polynomial([0]);
   1322             while (u.deg() >= 0) {
   1323                 t = Polynomial.divrem(v, u);
   1324                 q = t[0];
   1325                 v = u;
   1326                 u = t[1];
   1327                 t = c;
   1328                 c = a - q * c;
   1329                 a = t;
   1330             }
   1331             /* v = gcd(x, y) */
   1332             if (v.deg() > 0)
   1333                 throw RangeError("not invertible");
   1334             return Polynomial.divrem(a, y)[1];
   1335         },
   1336         roots(p) {
   1337             return poly_roots(p);
   1338         }
   1339     });
   1340 
   1341     /* Polynomial Modulo Q */
   1342 
   1343     PolyMod = function PolyMod(a, m) {
   1344         var obj, t;
   1345         if (new.target)
   1346             throw TypeError("not a constructor");
   1347         obj = Object.create(PolyMod.prototype);
   1348         if (m instanceof Polynomial) {
   1349             if (m.deg() <= 0)
   1350                 throw RangeError("the modulo cannot have a degree <= 0");
   1351             if (a instanceof RationalFunction) {
   1352                 return PolyMod(a.num, m) / a.den;
   1353             } else {
   1354                 a = Polynomial(a);
   1355                 t = Polynomial.divrem(a, m);
   1356                 a = t[1];
   1357             }
   1358         } else {
   1359             throw TypeError("invalid types");
   1360         }
   1361         obj.res = a;
   1362         obj.mod = m;
   1363         return obj;
   1364     };
   1365 
   1366     function polymod_add(a, b) {
   1367         if (!(a instanceof PolyMod)) {
   1368             return PolyMod(a + b.res, b.mod);
   1369         } else if (!(b instanceof PolyMod)) {
   1370             return PolyMod(a.res + b, a.mod);
   1371         } else {
   1372             if (a.mod != b.mod)
   1373                 throw TypeError("different modulo for binary operator");
   1374             return PolyMod(a.res + b.res, a.mod);
   1375         }
   1376     }
   1377     function polymod_sub(a, b) {
   1378         return polymod_add(a, -b);
   1379     }
   1380     function polymod_mul(a, b) {
   1381         if (!(a instanceof PolyMod)) {
   1382             return PolyMod(a * b.res, b.mod);
   1383         } else if (!(b instanceof PolyMod)) {
   1384             return PolyMod(a.res * b, a.mod);
   1385         } else {
   1386             if (a.mod != b.mod)
   1387                     throw TypeError("different modulo for binary operator");
   1388             return PolyMod(a.res * b.res, a.mod);
   1389         }
   1390     }
   1391     function polymod_div(a, b) {
   1392         if (!(b instanceof PolyMod))
   1393             b = PolyMod(b, a.mod);
   1394         return polymod_mul(a, b.inverse());
   1395     }
   1396     function polymod_eq(a, b) {
   1397         return (a.mod == b.mod && a.res == b.res);
   1398     }
   1399 
   1400     operators_set(PolyMod.prototype,
   1401         {
   1402             "+": polymod_add,
   1403             "-": polymod_sub,
   1404             "*": polymod_mul,
   1405             "/": polymod_div,
   1406             "**": generic_pow,
   1407             "==": polymod_eq,
   1408             "pos"(a) {
   1409                 return a;
   1410             },
   1411             "neg"(a) {
   1412                 return PolyMod(-a.res, a.mod);
   1413             },
   1414         },
   1415         {
   1416             left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
   1417             right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
   1418             "+": polymod_add,
   1419             "-": polymod_sub,
   1420             "*": polymod_mul,
   1421             "/": polymod_div,
   1422             "**": generic_pow, /* XXX: only for integer */
   1423         });
   1424 
   1425     add_props(PolyMod.prototype, {
   1426         inverse() {
   1427             var a = this, m = a.mod;
   1428             if (m instanceof Polynomial) {
   1429                 return PolyMod(Polynomial.invmod(a.res, m), m);
   1430             } else {
   1431                 throw TypeError("unsupported type");
   1432             }
   1433         },
   1434         toString() {
   1435             return "PolyMod(" + this.res + "," + this.mod + ")";
   1436         },
   1437     });
   1438 
   1439     /* Rational function */
   1440 
   1441     RationalFunction = function RationalFunction(a, b)
   1442     {
   1443         var t, r, d, obj;
   1444         if (new.target)
   1445             throw TypeError("not a constructor");
   1446         if (!(a instanceof Polynomial) ||
   1447             !(b instanceof Polynomial))
   1448             throw TypeError("polynomial expected");
   1449         t = Polynomial.divrem(a, b);
   1450         r = t[1];
   1451         if (r.deg() < 0)
   1452             return t[0]; /* no need for a fraction */
   1453         d = Polynomial.gcd(b, r);
   1454         if (d.deg() > 0) {
   1455             a = Polynomial.divrem(a, d)[0];
   1456             b = Polynomial.divrem(b, d)[0];
   1457         }
   1458         obj = Object.create(RationalFunction.prototype);
   1459         obj.num = a;
   1460         obj.den = b;
   1461         return obj;
   1462     }
   1463 
   1464     add_props(RationalFunction.prototype, {
   1465         inverse() {
   1466             return RationalFunction(this.den, this.num);
   1467         },
   1468         conj() {
   1469             return RationalFunction(this.num.conj(), this.den.conj());
   1470         },
   1471         toString() {
   1472             var str;
   1473             if (this.num.deg() <= 0 &&
   1474                 !number_need_paren(this.num[0]))
   1475                 str = this.num.toString();
   1476             else
   1477                 str = "(" + this.num.toString() + ")";
   1478             str += "/(" + this.den.toString() + ")"
   1479             return str;
   1480         },
   1481         apply(b) {
   1482             return this.num.apply(b) / this.den.apply(b);
   1483         },
   1484         deriv() {
   1485             var n = this.num, d = this.den;
   1486             return RationalFunction(n.deriv() * d - n * d.deriv(), d * d);
   1487         },
   1488     });
   1489 
   1490     function ratfunc_add(a, b) {
   1491         a = RationalFunction.toRationalFunction(a);
   1492         b = RationalFunction.toRationalFunction(b);
   1493         return RationalFunction(a.num * b.den + a.den * b.num, a.den * b.den);
   1494     }
   1495     function ratfunc_sub(a, b) {
   1496         a = RationalFunction.toRationalFunction(a);
   1497         b = RationalFunction.toRationalFunction(b);
   1498         return RationalFunction(a.num * b.den - a.den * b.num, a.den * b.den);
   1499     }
   1500     function ratfunc_mul(a, b) {
   1501         a = RationalFunction.toRationalFunction(a);
   1502         b = RationalFunction.toRationalFunction(b);
   1503         return RationalFunction(a.num * b.num, a.den * b.den);
   1504     }
   1505     function ratfunc_div(a, b) {
   1506         a = RationalFunction.toRationalFunction(a);
   1507         b = RationalFunction.toRationalFunction(b);
   1508         return RationalFunction(a.num * b.den, a.den * b.num);
   1509     }
   1510     function ratfunc_eq(a, b) {
   1511         a = RationalFunction.toRationalFunction(a);
   1512         b = RationalFunction.toRationalFunction(b);
   1513         /* we assume the fractions are normalized */
   1514         return (a.num == b.num && a.den == b.den);
   1515     }
   1516 
   1517     operators_set(RationalFunction.prototype,
   1518         {
   1519             "+": ratfunc_add,
   1520             "-": ratfunc_sub,
   1521             "*": ratfunc_mul,
   1522             "/": ratfunc_div,
   1523             "**": generic_pow,
   1524             "==": ratfunc_eq,
   1525             "pos"(a) {
   1526                 return a;
   1527             },
   1528             "neg"(a) {
   1529                 return RationalFunction(-this.num, this.den);
   1530             },
   1531         },
   1532         {
   1533             left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
   1534             right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
   1535             "+": ratfunc_add,
   1536             "-": ratfunc_sub,
   1537             "*": ratfunc_mul,
   1538             "/": ratfunc_div,
   1539             "**": generic_pow, /* should only be used with integers */
   1540         });
   1541 
   1542     add_props(RationalFunction, {
   1543         /* This function always return a RationalFunction object even
   1544            if it could simplified to a polynomial, so it is not
   1545            equivalent to RationalFunction(a) */
   1546         toRationalFunction(a) {
   1547             var obj;
   1548             if (a instanceof RationalFunction) {
   1549                 return a;
   1550             } else {
   1551                 obj = Object.create(RationalFunction.prototype);
   1552                 obj.num = Polynomial(a);
   1553                 obj.den = Polynomial(1);
   1554                 return obj;
   1555             }
   1556         },
   1557     });
   1558 
   1559     /* Power series */
   1560 
   1561     /* 'a' is an array */
   1562     function get_emin(a) {
   1563         var i, n;
   1564         n = a.length;
   1565         for(i = 0; i < n; i++) {
   1566             if (a[i] != 0)
   1567                 return i;
   1568         }
   1569         return n;
   1570     };
   1571 
   1572     function series_is_scalar_or_polynomial(a)
   1573     {
   1574         return polynomial_is_scalar(a) ||
   1575             (a instanceof Polynomial);
   1576     }
   1577 
   1578     /* n is the maximum number of terms if 'a' is not a serie */
   1579     Series = function Series(a, n) {
   1580         var emin, r, i;
   1581 
   1582         if (a instanceof Series) {
   1583             return a;
   1584         } else if (series_is_scalar_or_polynomial(a)) {
   1585             if (n <= 0) {
   1586                 /* XXX: should still use the polynomial degree */
   1587                 return Series.zero(0, 0);
   1588             } else {
   1589                 a = Polynomial(a);
   1590                 emin = get_emin(a);
   1591                 r = Series.zero(n, emin);
   1592                 n = Math.min(a.length - emin, n);
   1593                 for(i = 0; i < n; i++)
   1594                     r[i] = a[i + emin];
   1595                 return r;
   1596             }
   1597         } else if (a instanceof RationalFunction) {
   1598             return Series(a.num, n) / a.den;
   1599         } else {
   1600             throw TypeError("invalid type");
   1601         }
   1602     };
   1603 
   1604     function series_add(v1, v2) {
   1605         var tmp, d, emin, n, r, i, j, v2_emin, c1, c2;
   1606         if (!(v1 instanceof Series)) {
   1607             tmp = v1;
   1608             v1 = v2;
   1609             v2 = tmp;
   1610         }
   1611         d = v1.emin + v1.length;
   1612         if (series_is_scalar_or_polynomial(v2)) {
   1613             v2 = Polynomial(v2);
   1614             if (d <= 0)
   1615                 return v1;
   1616             v2_emin = 0;
   1617         } else if (v2 instanceof RationalFunction) {
   1618             /* compute the emin of the rational fonction */
   1619             i = get_emin(v2.num) - get_emin(v2.den);
   1620             if (d <= i)
   1621                 return v1;
   1622             /* compute the serie with the required terms */
   1623             v2 = Series(v2, d - i);
   1624             v2_emin = v2.emin;
   1625         } else {
   1626             v2_emin = v2.emin;
   1627             d = Math.min(d, v2_emin + v2.length);
   1628         }
   1629         emin = Math.min(v1.emin, v2_emin);
   1630         n = d - emin;
   1631         r = Series.zero(n, emin);
   1632         /* XXX: slow */
   1633         for(i = emin; i < d; i++) {
   1634             j = i - v1.emin;
   1635             if (j >= 0 && j < v1.length)
   1636                 c1 = v1[j];
   1637             else
   1638                 c1 = 0;
   1639             j = i - v2_emin;
   1640             if (j >= 0 && j < v2.length)
   1641                 c2 = v2[j];
   1642             else
   1643                 c2 = 0;
   1644             r[i - emin] = c1 + c2;
   1645         }
   1646         return r.trim();
   1647     }
   1648     function series_sub(a, b) {
   1649         return series_add(a, -b);
   1650     }
   1651     function series_mul(v1, v2) {
   1652         var n, i, j, r, n, emin, n1, n2, k;
   1653         if (!(v1 instanceof Series))
   1654             v1 = Series(v1, v2.length);
   1655         else if (!(v2 instanceof Series))
   1656             v2 = Series(v2, v1.length);
   1657         emin = v1.emin + v2.emin;
   1658         n = Math.min(v1.length, v2.length);
   1659         n1 = v1.length;
   1660         n2 = v2.length;
   1661         r = Series.zero(n, emin);
   1662         for(i = 0; i < n1; i++) {
   1663             k = Math.min(n2, n - i);
   1664             for(j = 0; j < k; j++) {
   1665                 r[i + j] += v1[i] * v2[j];
   1666             }
   1667         }
   1668         return r.trim();
   1669     }
   1670     function series_div(v1, v2) {
   1671         if (!(v2 instanceof Series))
   1672             v2 = Series(v2, v1.length);
   1673         return series_mul(v1, v2.inverse());
   1674     }
   1675     function series_pow(a, b) {
   1676         if (Integer.isInteger(b)) {
   1677             return generic_pow(a, b);
   1678         } else {
   1679             if (!(a instanceof Series))
   1680                 a = Series(a, b.length);
   1681             return exp(log(a) * b);
   1682         }
   1683     }
   1684     function series_eq(a, b) {
   1685         var n, i;
   1686         if (a.emin != b.emin)
   1687             return false;
   1688         n = a.length;
   1689         if (n != b.length)
   1690             return false;
   1691         for(i = 0; i < n; i++) {
   1692             if (a[i] != b[i])
   1693                 return false;
   1694         }
   1695         return true;
   1696     }
   1697 
   1698     operators_set(Series.prototype,
   1699         {
   1700             "+": series_add,
   1701             "-": series_sub,
   1702             "*": series_mul,
   1703             "/": series_div,
   1704             "**": series_pow,
   1705             "==": series_eq,
   1706             "pos"(a) {
   1707                 return a;
   1708             },
   1709             "neg"(a) {
   1710                 var obj, n, i;
   1711                 n = a.length;
   1712                 obj = Series.zero(a.length, a.emin);
   1713                 for(i = 0; i < n; i++) {
   1714                     obj[i] = -a[i];
   1715                 }
   1716                 return obj;
   1717             },
   1718         },
   1719         {
   1720             left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
   1721             right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
   1722             "+": series_add,
   1723             "-": series_sub,
   1724             "*": series_mul,
   1725             "/": series_div,
   1726             "**": series_pow,
   1727         });
   1728 
   1729     add_props(Series.prototype, {
   1730         conj() {
   1731             var obj, n, i;
   1732             n = this.length;
   1733             obj = Series.zero(this.length, this.emin);
   1734             for(i = 0; i < n; i++) {
   1735                 obj[i] = this[i].conj();
   1736             }
   1737             return obj;
   1738         },
   1739         inverse() {
   1740             var r, n, i, j, sum, v1 = this;
   1741             n = v1.length;
   1742             if (n == 0)
   1743                 throw RangeError("division by zero");
   1744             r = Series.zero(n, -v1.emin);
   1745             r[0] = 1 / v1[0];
   1746             for(i = 1; i < n; i++) {
   1747                 sum = 0;
   1748                 for(j = 1; j <= i; j++) {
   1749                     sum += v1[j] * r[i - j];
   1750                 }
   1751                 r[i] = -sum * r[0];
   1752             }
   1753             return r;
   1754         },
   1755         /* remove leading zero terms */
   1756         trim() {
   1757             var i, j, n, r, v1 = this;
   1758             n = v1.length;
   1759             i = 0;
   1760             while (i < n && v1[i] == 0)
   1761                 i++;
   1762             if (i == 0)
   1763                 return v1;
   1764             for(j = i; j < n; j++)
   1765                 v1[j - i] = v1[j];
   1766             v1.length = n - i;
   1767             v1.__proto__.emin += i;
   1768             return v1;
   1769         },
   1770         toString() {
   1771             var i, j, str, str1, c, a = this, emin, n;
   1772             str="";
   1773             emin = this.emin;
   1774             n = this.length;
   1775             for(j = 0; j < n; j++) {
   1776                 i = j + emin;
   1777                 c = a[j];
   1778                 if (c != 0) {
   1779                     str1 = monomial_toString(c, i);
   1780                     if (str1[0] != "-") {
   1781                         if (str != "")
   1782                             str += "+";
   1783                     }
   1784                     str += str1;
   1785                 }
   1786             }
   1787             if (str != "")
   1788                 str += "+";
   1789             str += "O(" + monomial_toString(1, n + emin) + ")";
   1790             return str;
   1791         },
   1792         apply(b) {
   1793             var i, n, r, a = this;
   1794             n = a.length;
   1795             if (n == 0)
   1796                 return 0;
   1797             r = a[--n];
   1798             while (n > 0) {
   1799                 n--;
   1800                 r = r * b + a[n];
   1801             }
   1802             if (a.emin != 0)
   1803                 r *= b ^ a.emin;
   1804             return r;
   1805         },
   1806         deriv() {
   1807             var a = this, n = a.length, emin = a.emin, r, i, j;
   1808             if (n == 0 && emin == 0) {
   1809                 return Series.zero(0, 0);
   1810             } else {
   1811                 r = Series.zero(n, emin - 1);
   1812                 for(i = 0; i < n; i++) {
   1813                     j = emin + i;
   1814                     if (j == 0)
   1815                         r[i] = 0;
   1816                     else
   1817                         r[i] = j * a[i];
   1818                 }
   1819                 return r.trim();
   1820             }
   1821         },
   1822         integ() {
   1823             var a = this, n = a.length, emin = a.emin, i, j, r;
   1824             r = Series.zero(n, emin + 1);
   1825             for(i = 0; i < n; i++) {
   1826                 j = emin + i;
   1827                 if (j == -1) {
   1828                     if (a[i] != 0)
   1829                         throw RangeError("cannot represent integ(1/X)");
   1830                 } else {
   1831                     r[i] = a[i] / (j + 1);
   1832                 }
   1833             }
   1834             return r.trim();
   1835         },
   1836         exp() {
   1837             var c, i, r, n, a = this;
   1838             if (a.emin < 0)
   1839                 throw RangeError("negative exponent in exp");
   1840             n = a.emin + a.length;
   1841             if (a.emin > 0 || a[0] == 0) {
   1842                 c = 1;
   1843             } else {
   1844                 c = global.exp(a[0]);
   1845                 a -= a[0];
   1846             }
   1847             r = Series.zero(n, 0);
   1848             for(i = 0; i < n; i++) {
   1849                 r[i] = c / fact(i);
   1850             }
   1851             return r.apply(a);
   1852         },
   1853         log() {
   1854             var a = this, r;
   1855             if (a.emin != 0)
   1856                 throw RangeError("log argument must have a non zero constant term");
   1857             r = integ(deriv(a) / a);
   1858             /* add the constant term */
   1859             r += global.log(a[0]);
   1860             return r;
   1861         },
   1862     });
   1863 
   1864     add_props(Series, {
   1865         /* new series of length n and first exponent emin */
   1866         zero(n, emin) {
   1867             var r, i, obj;
   1868 
   1869             r = [];
   1870             for(i = 0; i < n; i++)
   1871                 r[i] = 0;
   1872             /* we return an array and store emin in its prototype */
   1873             obj = Object.create(Series.prototype);
   1874             obj.emin = emin;
   1875             Object.setPrototypeOf(r, obj);
   1876             return r;
   1877         },
   1878         O(a) {
   1879             function ErrorO() {
   1880                 return TypeError("invalid O() argument");
   1881             }
   1882             var n;
   1883             if (series_is_scalar_or_polynomial(a)) {
   1884                 a = Polynomial(a);
   1885                 n = a.deg();
   1886                 if (n < 0)
   1887                     throw ErrorO();
   1888             } else if (a instanceof RationalFunction) {
   1889                 if (a.num.deg() != 0)
   1890                     throw ErrorO();
   1891                 n = a.den.deg();
   1892                 if (n < 0)
   1893                     throw ErrorO();
   1894                 n = -n;
   1895             } else
   1896                 throw ErrorO();
   1897             return Series.zero(0, n);
   1898         },
   1899     });
   1900 
   1901     /* Array (Matrix) */
   1902 
   1903     Matrix = function Matrix(h, w) {
   1904         var i, j, r, rl;
   1905         if (typeof w === "undefined")
   1906             w = h;
   1907         r = [];
   1908         for(i = 0; i < h; i++) {
   1909             rl = [];
   1910             for(j = 0; j < w; j++)
   1911                 rl[j] = 0;
   1912             r[i] = rl;
   1913         }
   1914         return r;
   1915     };
   1916 
   1917     add_props(Matrix, {
   1918         idn(n) {
   1919             var r, i;
   1920             r = Matrix(n, n);
   1921             for(i = 0; i < n; i++)
   1922                 r[i][i] = 1;
   1923             return r;
   1924         },
   1925         diag(a) {
   1926             var r, i, n;
   1927             n = a.length;
   1928             r = Matrix(n, n);
   1929             for(i = 0; i < n; i++)
   1930                 r[i][i] = a[i];
   1931             return r;
   1932         },
   1933         hilbert(n) {
   1934             var i, j, r;
   1935             r = Matrix(n);
   1936             for(i = 0; i < n; i++) {
   1937                 for(j = 0; j < n; j++) {
   1938                     r[i][j] = 1 / (1 + i + j);
   1939                 }
   1940             }
   1941             return r;
   1942         },
   1943         trans(a) {
   1944             var h, w, r, i, j;
   1945             if (!Array.isArray(a))
   1946                 throw TypeError("matrix expected");
   1947             h = a.length;
   1948             if (!Array.isArray(a[0])) {
   1949                 w = 1;
   1950                 r = Matrix(w, h);
   1951                 for(i = 0; i < h; i++) {
   1952                     r[0][i] = a[i];
   1953                 }
   1954             } else {
   1955                 w = a[0].length;
   1956                 r = Matrix(w, h);
   1957                 for(i = 0; i < h; i++) {
   1958                     for(j = 0; j < w; j++) {
   1959                         r[j][i] = a[i][j];
   1960                     }
   1961                 }
   1962             }
   1963             return r;
   1964         },
   1965         check_square(a) {
   1966             var a, n;
   1967             if (!Array.isArray(a))
   1968                 throw TypeError("array expected");
   1969             n = a.length;
   1970             if (!Array.isArray(a[0]) || n != a[0].length)
   1971                 throw TypeError("square matrix expected");
   1972             return n;
   1973         },
   1974         trace(a) {
   1975             var n, r, i;
   1976             n = Matrix.check_square(a);
   1977             r = a[0][0];
   1978             for(i = 1; i < n; i++) {
   1979                 r += a[i][i];
   1980             }
   1981             return r;
   1982         },
   1983         charpoly(a) {
   1984             var n, p, c, i, j, coef;
   1985             n = Matrix.check_square(a);
   1986             p = [];
   1987             for(i = 0; i < n + 1; i++)
   1988                 p[i] = 0;
   1989             p[n] = 1;
   1990             c = Matrix.idn(n);
   1991             for(i = 0; i < n; i++) {
   1992                 c = c * a;
   1993                 coef = -trace(c) / (i + 1);
   1994                 p[n - i - 1] = coef;
   1995                 for(j = 0; j < n; j++)
   1996                     c[j][j] += coef;
   1997             }
   1998             return Polynomial(p);
   1999         },
   2000         eigenvals(a) {
   2001             return Polynomial.roots(Matrix.charpoly(a));
   2002         },
   2003         det(a) {
   2004             var n, i, j, k, s, src, v, c;
   2005 
   2006             n = Matrix.check_square(a);
   2007             s = 1;
   2008             src = a.dup();
   2009             for(i=0;i<n;i++) {
   2010                 for(j = i; j < n; j++) {
   2011                     if (src[j][i] != 0)
   2012                         break;
   2013                 }
   2014                 if (j == n)
   2015                     return 0;
   2016                 if (j != i) {
   2017                     for(k = 0;k < n; k++) {
   2018                         v = src[j][k];
   2019                         src[j][k] = src[i][k];
   2020                         src[i][k] = v;
   2021                     }
   2022                     s = -s;
   2023                 }
   2024                 c = src[i][i].inverse();
   2025                 for(j = i + 1; j < n; j++) {
   2026                     v = c * src[j][i];
   2027                     for(k = 0;k < n; k++) {
   2028                         src[j][k] -= src[i][k] * v;
   2029                     }
   2030                 }
   2031             }
   2032             c = s;
   2033             for(i=0;i<n;i++)
   2034                 c *= src[i][i];
   2035             return c;
   2036         },
   2037         inverse(a) {
   2038             var n, dst, src, i, j, k, n2, r, c, v;
   2039             n = Matrix.check_square(a);
   2040             src = a.dup();
   2041             dst = Matrix.idn(n);
   2042             for(i=0;i<n;i++) {
   2043                 for(j = i; j < n; j++) {
   2044                     if (src[j][i] != 0)
   2045                         break;
   2046                 }
   2047                 if (j == n)
   2048                     throw RangeError("matrix is not invertible");
   2049                 if (j != i) {
   2050                     /* swap lines in src and dst */
   2051                     v = src[j];
   2052                     src[j] = src[i];
   2053                     src[i] = v;
   2054                     v = dst[j];
   2055                     dst[j] = dst[i];
   2056                     dst[i] = v;
   2057                 }
   2058 
   2059                 c = src[i][i].inverse();
   2060                 for(k = 0; k < n; k++) {
   2061                     src[i][k] *= c;
   2062                     dst[i][k] *= c;
   2063                 }
   2064 
   2065                 for(j = 0; j < n; j++) {
   2066                     if (j != i) {
   2067                         c = src[j][i];
   2068                         for(k = i; k < n; k++) {
   2069                             src[j][k] -= src[i][k] * c;
   2070                         }
   2071                         for(k = 0; k < n; k++) {
   2072                             dst[j][k] -= dst[i][k] * c;
   2073                         }
   2074                     }
   2075                 }
   2076             }
   2077             return dst;
   2078         },
   2079         rank(a) {
   2080             var src, i, j, k, w, h, l, c;
   2081 
   2082             if (!Array.isArray(a) ||
   2083                 !Array.isArray(a[0]))
   2084                 throw TypeError("matrix expected");
   2085             h = a.length;
   2086             w = a[0].length;
   2087             src = a.dup();
   2088             l = 0;
   2089             for(i=0;i<w;i++) {
   2090                 for(j = l; j < h; j++) {
   2091                     if (src[j][i] != 0)
   2092                         break;
   2093                 }
   2094                 if (j == h)
   2095                     continue;
   2096                 if (j != l) {
   2097                     /* swap lines */
   2098                     for(k = 0; k < w; k++) {
   2099                         v = src[j][k];
   2100                         src[j][k] = src[l][k];
   2101                         src[l][k] = v;
   2102                     }
   2103                 }
   2104 
   2105                 c = src[l][i].inverse();
   2106                 for(k = 0; k < w; k++) {
   2107                     src[l][k] *= c;
   2108                 }
   2109 
   2110                 for(j = l + 1; j < h; j++) {
   2111                     c = src[j][i];
   2112                     for(k = i; k < w; k++) {
   2113                         src[j][k] -= src[l][k] * c;
   2114                     }
   2115                 }
   2116                 l++;
   2117             }
   2118             return l;
   2119         },
   2120         ker(a) {
   2121             var src, i, j, k, w, h, l, m, r, im_cols, ker_dim, c;
   2122 
   2123             if (!Array.isArray(a) ||
   2124                 !Array.isArray(a[0]))
   2125                 throw TypeError("matrix expected");
   2126             h = a.length;
   2127             w = a[0].length;
   2128             src = a.dup();
   2129             im_cols = [];
   2130             l = 0;
   2131             for(i=0;i<w;i++) {
   2132                 im_cols[i] = false;
   2133                 for(j = l; j < h; j++) {
   2134                     if (src[j][i] != 0)
   2135                         break;
   2136                 }
   2137                 if (j == h)
   2138                     continue;
   2139                 im_cols[i] = true;
   2140                 if (j != l) {
   2141                     /* swap lines */
   2142                     for(k = 0; k < w; k++) {
   2143                         v = src[j][k];
   2144                         src[j][k] = src[l][k];
   2145                         src[l][k] = v;
   2146                     }
   2147                 }
   2148 
   2149                 c = src[l][i].inverse();
   2150                 for(k = 0; k < w; k++) {
   2151                     src[l][k] *= c;
   2152                 }
   2153 
   2154                 for(j = 0; j < h; j++) {
   2155                     if (j != l) {
   2156                         c = src[j][i];
   2157                         for(k = i; k < w; k++) {
   2158                             src[j][k] -= src[l][k] * c;
   2159                         }
   2160                     }
   2161                 }
   2162                 l++;
   2163                 //        log_str("m=" + cval_toString(v1) + "\n");
   2164             }
   2165             //    log_str("im cols="+im_cols+"\n");
   2166 
   2167             /* build the kernel vectors */
   2168             ker_dim = w - l;
   2169             r = Matrix(w, ker_dim);
   2170             k = 0;
   2171             for(i = 0; i < w; i++) {
   2172                 if (!im_cols[i]) {
   2173                     /* select this column from the matrix */
   2174                     l = 0;
   2175                     m = 0;
   2176                     for(j = 0; j < w; j++) {
   2177                         if (im_cols[j]) {
   2178                             r[j][k] = -src[m][i];
   2179                             m++;
   2180                         } else {
   2181                             if (l == k) {
   2182                                 r[j][k] = 1;
   2183                             } else {
   2184                                 r[j][k] = 0;
   2185                             }
   2186                             l++;
   2187                         }
   2188                     }
   2189                     k++;
   2190                 }
   2191             }
   2192             return r;
   2193         },
   2194         dp(a, b) {
   2195             var i, n, r;
   2196             n = a.length;
   2197             if (n != b.length)
   2198                 throw TypeError("incompatible array length");
   2199             /* XXX: could do complex product */
   2200             r = 0;
   2201             for(i = 0; i < n; i++) {
   2202                 r += a[i] * b[i];
   2203             }
   2204             return r;
   2205         },
   2206         /* cross product */
   2207         cp(v1, v2) {
   2208             var r;
   2209             if (v1.length != 3 || v2.length != 3)
   2210                 throw TypeError("vectors must have 3 elements");
   2211             r = [];
   2212             r[0] = v1[1] * v2[2] - v1[2] * v2[1];
   2213             r[1] = v1[2] * v2[0] - v1[0] * v2[2];
   2214             r[2] = v1[0] * v2[1] - v1[1] * v2[0];
   2215             return r;
   2216         },
   2217     });
   2218 
   2219     function array_add(a, b) {
   2220         var r, i, n;
   2221         n = a.length;
   2222         if (n != b.length)
   2223             throw TypeError("incompatible array size");
   2224         r = [];
   2225         for(i = 0; i < n; i++)
   2226             r[i] = a[i] + b[i];
   2227         return r;
   2228     }
   2229     function array_sub(a, b) {
   2230         var r, i, n;
   2231         n = a.length;
   2232         if (n != b.length)
   2233             throw TypeError("incompatible array size");
   2234         r = [];
   2235         for(i = 0; i < n; i++)
   2236             r[i] = a[i] - b[i];
   2237         return r;
   2238     }
   2239     function array_scalar_mul(a, b) {
   2240         var r, i, n;
   2241         n = a.length;
   2242         r = [];
   2243         for(i = 0; i < n; i++)
   2244             r[i] = a[i] * b;
   2245         return r;
   2246     }
   2247     function array_mul(a, b) {
   2248         var h, w, l, i, j, k, r, rl, sum, a_mat, b_mat;
   2249         h = a.length;
   2250         a_mat = Array.isArray(a[0]);
   2251         if (a_mat) {
   2252             l = a[0].length;
   2253         } else {
   2254             l = 1;
   2255         }
   2256         if (l != b.length)
   2257             throw RangeError("incompatible matrix size");
   2258         b_mat = Array.isArray(b[0]);
   2259         if (b_mat)
   2260             w = b[0].length;
   2261         else
   2262             w = 1;
   2263         r = [];
   2264         if (a_mat && b_mat) {
   2265             for(i = 0; i < h; i++) {
   2266                 rl = [];
   2267                 for(j = 0; j < w; j++) {
   2268                     sum = 0;
   2269                     for(k = 0; k < l; k++) {
   2270                         sum += a[i][k] * b[k][j];
   2271                     }
   2272                     rl[j] = sum;
   2273                 }
   2274                 r[i] = rl;
   2275             }
   2276         } else if (a_mat && !b_mat) {
   2277             for(i = 0; i < h; i++) {
   2278                 sum = 0;
   2279                 for(k = 0; k < l; k++) {
   2280                     sum += a[i][k] * b[k];
   2281                 }
   2282                 r[i] = sum;
   2283             }
   2284         } else if (!a_mat && b_mat) {
   2285             for(i = 0; i < h; i++) {
   2286                 rl = [];
   2287                 for(j = 0; j < w; j++) {
   2288                     rl[j] = a[i] * b[0][j];
   2289                 }
   2290                 r[i] = rl;
   2291             }
   2292         } else {
   2293             for(i = 0; i < h; i++) {
   2294                 r[i] = a[i] * b[0];
   2295             }
   2296         }
   2297         return r;
   2298     }
   2299     function array_div(a, b) {
   2300         return array_mul(a, b.inverse());
   2301     }
   2302     function array_element_wise_inverse(a) {
   2303         var r, i, n;
   2304         n = a.length;
   2305         r = [];
   2306         for(i = 0; i < n; i++)
   2307             r[i] = a[i].inverse();
   2308         return r;
   2309     }
   2310     function array_eq(a, b) {
   2311         var n, i;
   2312         n = a.length;
   2313         if (n != b.length)
   2314             return false;
   2315         for(i = 0; i < n; i++) {
   2316             if (a[i] != b[i])
   2317                 return false;
   2318         }
   2319         return true;
   2320     }
   2321 
   2322     operators_set(Array.prototype,
   2323         {
   2324             "+": array_add,
   2325             "-": array_sub,
   2326             "*": array_mul,
   2327             "/": array_div,
   2328             "==": array_eq,
   2329             "pos"(a) {
   2330                 return a;
   2331             },
   2332             "neg"(a) {
   2333                 var i, n, r;
   2334                 n = a.length;
   2335                 r = [];
   2336                 for(i = 0; i < n; i++)
   2337                     r[i] = -a[i];
   2338                 return r;
   2339             }
   2340         },
   2341         {
   2342             right: [Number, BigInt, Float, Fraction, Complex, Mod,
   2343                     Polynomial, PolyMod, RationalFunction, Series],
   2344             "*": array_scalar_mul,
   2345             "/"(a, b) { return a * b.inverse(); },
   2346             "**": generic_pow, /* XXX: only for integer */
   2347         },
   2348         {
   2349             left: [Number, BigInt, Float, Fraction, Complex, Mod,
   2350                    Polynomial, PolyMod, RationalFunction, Series],
   2351             "*"(a, b) { return array_scalar_mul(b, a); },
   2352             "/"(a, b) { return a * array_element_wise_inverse(b); },
   2353         });
   2354 
   2355     add_props(Array.prototype, {
   2356         conj() {
   2357             var i, n, r;
   2358             n = this.length;
   2359             r = [];
   2360             for(i = 0; i < n; i++)
   2361                 r[i] = this[i].conj();
   2362             return r;
   2363         },
   2364         dup() {
   2365             var r, i, n, el, a = this;
   2366             r = [];
   2367             n = a.length;
   2368             for(i = 0; i < n; i++) {
   2369                 el = a[i];
   2370                 if (Array.isArray(el))
   2371                     el = el.dup();
   2372                 r[i] = el;
   2373             }
   2374             return r;
   2375         },
   2376         inverse() {
   2377             return Matrix.inverse(this);
   2378         },
   2379         norm2: Polynomial.prototype.norm2,
   2380     });
   2381 
   2382 })(this);
   2383 
   2384 /* global definitions */
   2385 var I = Complex(0, 1);
   2386 var X = Polynomial([0, 1]);
   2387 var O = Series.O;
   2388 
   2389 Object.defineProperty(this, "PI", { get: function () { return Float.PI } });
   2390 
   2391 /* put frequently used functions in the global context */
   2392 var gcd = Integer.gcd;
   2393 var fact = Integer.fact;
   2394 var comb = Integer.comb;
   2395 var pmod = Integer.pmod;
   2396 var invmod = Integer.invmod;
   2397 var factor = Integer.factor;
   2398 var isprime = Integer.isPrime;
   2399 var nextprime = Integer.nextPrime;
   2400 
   2401 function deriv(a)
   2402 {
   2403     return a.deriv();
   2404 }
   2405 
   2406 function integ(a)
   2407 {
   2408     return a.integ();
   2409 }
   2410 
   2411 function norm2(a)
   2412 {
   2413     return a.norm2();
   2414 }
   2415 
   2416 function abs(a)
   2417 {
   2418     return a.abs();
   2419 }
   2420 
   2421 function conj(a)
   2422 {
   2423     return a.conj();
   2424 }
   2425 
   2426 function arg(a)
   2427 {
   2428     return a.arg();
   2429 }
   2430 
   2431 function inverse(a)
   2432 {
   2433     return a.inverse();
   2434 }
   2435 
   2436 function trunc(a)
   2437 {
   2438     if (Integer.isInteger(a)) {
   2439         return a;
   2440     } else if (a instanceof Fraction) {
   2441         return Integer.tdiv(a.num, a.den);
   2442     } else if (a instanceof Polynomial) {
   2443         return a;
   2444     } else if (a instanceof RationalFunction) {
   2445         return Polynomial.divrem(a.num, a.den)[0];
   2446     } else {
   2447         return Float.ceil(a);
   2448     }
   2449 }
   2450 
   2451 function floor(a)
   2452 {
   2453     if (Integer.isInteger(a)) {
   2454         return a;
   2455     } else if (a instanceof Fraction) {
   2456         return Integer.fdiv(a.num, a.den);
   2457     } else {
   2458         return Float.floor(a);
   2459     }
   2460 }
   2461 
   2462 function ceil(a)
   2463 {
   2464     if (Integer.isInteger(a)) {
   2465         return a;
   2466     } else if (a instanceof Fraction) {
   2467         return Integer.cdiv(a.num, a.den);
   2468     } else {
   2469         return Float.ceil(a);
   2470     }
   2471 }
   2472 
   2473 function sqrt(a)
   2474 {
   2475     var t, u, re, im;
   2476     if (a instanceof Series) {
   2477         return a ^ (1/2);
   2478     } else if (a instanceof Complex) {
   2479         t = abs(a);
   2480         u = a.re;
   2481         re = sqrt((t + u) / 2);
   2482         im = sqrt((t - u) / 2);
   2483         if (a.im < 0)
   2484             im = -im;
   2485         return Complex.toComplex(re, im);
   2486     } else {
   2487         a = Float(a);
   2488         if (a < 0) {
   2489             return Complex(0, Float.sqrt(-a));
   2490         } else {
   2491             return Float.sqrt(a);
   2492         }
   2493     }
   2494 }
   2495 
   2496 function exp(a)
   2497 {
   2498     return a.exp();
   2499 }
   2500 
   2501 function log(a)
   2502 {
   2503     return a.log();
   2504 }
   2505 
   2506 function log2(a)
   2507 {
   2508     return log(a) * Float.LOG2E;
   2509 }
   2510 
   2511 function log10(a)
   2512 {
   2513     return log(a) * Float.LOG10E;
   2514 }
   2515 
   2516 function todb(a)
   2517 {
   2518     return log10(a) * 10;
   2519 }
   2520 
   2521 function fromdb(a)
   2522 {
   2523     return 10 ^ (a / 10);
   2524 }
   2525 
   2526 function sin(a)
   2527 {
   2528     var t;
   2529     if (a instanceof Complex || a instanceof Series) {
   2530         t = exp(a * I);
   2531         return (t - 1/t) / (2 * I);
   2532     } else {
   2533         return Float.sin(Float(a));
   2534     }
   2535 }
   2536 
   2537 function cos(a)
   2538 {
   2539     var t;
   2540     if (a instanceof Complex || a instanceof Series) {
   2541         t = exp(a * I);
   2542         return (t + 1/t) / 2;
   2543     } else {
   2544         return Float.cos(Float(a));
   2545     }
   2546 }
   2547 
   2548 function tan(a)
   2549 {
   2550     if (a instanceof Complex || a instanceof Series) {
   2551         return sin(a) / cos(a);
   2552     } else {
   2553         return Float.tan(Float(a));
   2554     }
   2555 }
   2556 
   2557 function asin(a)
   2558 {
   2559     return Float.asin(Float(a));
   2560 }
   2561 
   2562 function acos(a)
   2563 {
   2564     return Float.acos(Float(a));
   2565 }
   2566 
   2567 function atan(a)
   2568 {
   2569     return Float.atan(Float(a));
   2570 }
   2571 
   2572 function atan2(a, b)
   2573 {
   2574     return Float.atan2(Float(a), Float(b));
   2575 }
   2576 
   2577 function sinc(a)
   2578 {
   2579     if (a == 0) {
   2580         return 1;
   2581     } else {
   2582         a *= Float.PI;
   2583         return sin(a) / a;
   2584     }
   2585 }
   2586 
   2587 function todeg(a)
   2588 {
   2589     return a * 180 / Float.PI;
   2590 }
   2591 
   2592 function fromdeg(a)
   2593 {
   2594     return a * Float.PI / 180;
   2595 }
   2596 
   2597 function sinh(a)
   2598 {
   2599     var e = Float.exp(Float(a));
   2600     return (e - 1/e) * 0.5;
   2601 }
   2602 
   2603 function cosh(a)
   2604 {
   2605     var e = Float.exp(Float(a));
   2606     return (e + 1/e) * 0.5;
   2607 }
   2608 
   2609 function tanh(a)
   2610 {
   2611     var e = Float.exp(Float(a) * 2);
   2612     return (e - 1) / (e + 1);
   2613 }
   2614 
   2615 function asinh(a)
   2616 {
   2617     var x = Float(a);
   2618     return log(sqrt(x * x + 1) + x);
   2619 }
   2620 
   2621 function acosh(a)
   2622 {
   2623     var x = Float(a);
   2624     return log(sqrt(x * x - 1) + x);
   2625 }
   2626 
   2627 function atanh(a)
   2628 {
   2629     var x = Float(a);
   2630     return 0.5 * log((1 + x) / (1 - x));
   2631 }
   2632 
   2633 function sigmoid(x)
   2634 {
   2635     x = Float(x);
   2636     return 1 / (1 + exp(-x));
   2637 }
   2638 
   2639 function lerp(a, b, t)
   2640 {
   2641     return a + (b - a) * t;
   2642 }
   2643 
   2644 var idn = Matrix.idn;
   2645 var diag = Matrix.diag;
   2646 var trans = Matrix.trans;
   2647 var trace = Matrix.trace;
   2648 var charpoly = Matrix.charpoly;
   2649 var eigenvals = Matrix.eigenvals;
   2650 var det = Matrix.det;
   2651 var rank = Matrix.rank;
   2652 var ker = Matrix.ker;
   2653 var cp = Matrix.cp;
   2654 var dp = Matrix.dp;
   2655 
   2656 var polroots = Polynomial.roots;
   2657 var bestappr = Float.bestappr;