"use strict"
;
"use math"
;
var
Integer, Float, Fraction, Complex, Mod, Polynomial, PolyMod, RationalFunction, Series, Matrix;
(
function
(global) {
global.Integer = global.BigInt;
global.Float = global.BigFloat;
global.algebraicMode =
true
;
function
add_props(obj, props) {
var
i, val, prop, tab, desc;
tab = Reflect.ownKeys(props);
for
(i = 0; i < tab.length; i++) {
prop = tab[i];
desc = Object.getOwnPropertyDescriptor(props, prop);
desc.enumerable =
false
;
if
(
"value"
in
desc) {
if
(
typeof
desc.value !==
"function"
) {
desc.writable =
false
;
desc.configurable =
false
;
}
}
else
{
desc.configurable =
false
;
}
Object.defineProperty(obj, prop, desc);
}
}
function
operators_set(proto, ...op_list)
{
var
new_op_list, i, a, j, b, k, obj, tab;
var
fields = [
"left"
,
"right"
];
new_op_list = [];
for
(i = 0; i < op_list.length; i++) {
a = op_list[i];
if
(a.left || a.right) {
tab = [ a.left, a.right ];
delete
a.left;
delete
a.right;
for
(k = 0; k < 2; k++) {
obj = tab[k];
if
(obj) {
if
(!Array.isArray(obj)) {
obj = [ obj ];
}
for
(j = 0; j < obj.length; j++) {
b = {};
Object.assign(b, a);
b[fields[k]] = obj[j];
new_op_list.push(b);
}
}
}
}
else
{
new_op_list.push(a);
}
}
proto[Symbol.operatorSet] =
Operators.create.call(
null
, ...new_op_list);
}
function
generic_pow(a, b) {
var
r, is_neg, i;
if
(!Integer.isInteger(b)) {
return
exp(log(a) * b);
}
if
(Array.isArray(a) && !(a
instanceof
Polynomial ||
a
instanceof
Series)) {
r = idn(Matrix.check_square(a));
}
else
{
r = 1;
}
if
(b == 0)
return
r;
is_neg =
false
;
if
(b < 0) {
is_neg =
true
;
b = -b;
}
r = a;
for
(i = Integer.floorLog2(b) - 1; i >= 0; i--) {
r *= r;
if
((b >> i) & 1)
r *= a;
}
if
(is_neg) {
if
(
typeof
r.inverse !=
"function"
)
throw
"negative powers are not supported for this type"
;
r = r.inverse();
}
return
r;
}
var
small_primes = [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499 ];
function
miller_rabin_test(n, t) {
var
d, r, s, i, j, a;
d = n - 1;
s = 0;
while
((d & 1) == 0) {
d >>= 1;
s++;
}
if
(small_primes.length < t)
t = small_primes.length;
loop:
for
(j = 0; j < t; j++) {
a = small_primes[j];
r = Integer.pmod(a, d, n);
if
(r == 1 || r == (n - 1))
continue
;
for
(i = 1; i < s; i++) {
r = (r * r) % n;
if
(r == 1)
return
false
;
if
(r == (n - 1))
continue
loop;
}
return
false
;
}
return
true
;
}
function
fact_rec(a, b) {
var
i, r;
if
((b - a) <= 5) {
r = a;
for
(i = a + 1; i <= b; i++)
r *= i;
return
r;
}
else
{
i = (a + b) >> 1;
return
fact_rec(a, i) * fact_rec(i + 1, b);
}
}
Operators.updateBigIntOperators(
{
"/"
(a, b) {
if
(algebraicMode) {
return
Fraction.toFraction(a, b);
}
else
{
return
Float(a) / Float(b);
}
},
"**"
(a, b) {
if
(algebraicMode) {
return
generic_pow(a, b);
}
else
{
return
Float(a) ** Float(b);
}
}
});
add_props(Integer, {
isInteger(a) {
return
typeof
a ===
"bigint"
||
(
typeof
a ===
"number"
&& Number.isSafeInteger(a));
},
gcd(a, b) {
var
r;
while
(b != 0) {
r = a % b;
a = b;
b = r;
}
return
a;
},
fact(n) {
return
n <= 0 ? 1 : fact_rec(1, n);
},
comb(n, k) {
if
(k < 0 || k > n)
return
0;
if
(k > n - k)
k = n - k;
if
(k == 0)
return
1;
return
Integer.tdiv(fact_rec(n - k + 1, n), fact_rec(1, k));
},
invmod(x, y) {
var
q, u, v, a, c, t;
u = x;
v = y;
c = 1;
a = 0;
while
(u != 0) {
t = Integer.fdivrem(v, u);
q = t[0];
v = u;
u = t[1];
t = c;
c = a - q * c;
a = t;
}
if
(v != 1)
throw
RangeError(
"not invertible"
);
return
a % y;
},
pmod(a, b, m) {
var
r;
if
(b == 0)
return
1;
if
(b < 0) {
a = Integer.invmod(a, m);
b = -b;
}
r = 1;
for
(;;) {
if
(b & 1) {
r = (r * a) % m;
}
b >>= 1;
if
(b == 0)
break
;
a = (a * a) % m;
}
return
r;
},
isPrime(n, t) {
var
i, d, n1;
if
(!Integer.isInteger(n))
throw
TypeError(
"invalid type"
);
if
(n <= 1)
return
false
;
n1 = small_primes.length;
for
(i = 0; i < n1; i++) {
d = small_primes[i];
if
(d == n)
return
true
;
if
(d > n)
return
false
;
if
((n % d) == 0)
return
false
;
}
if
(n < d * d)
return
true
;
if
(
typeof
t ==
"undefined"
)
t = 64;
return
miller_rabin_test(n, t);
},
nextPrime(n) {
if
(!Integer.isInteger(n))
throw
TypeError(
"invalid type"
);
if
(n < 1)
n = 1;
for
(;;) {
n++;
if
(Integer.isPrime(n))
return
n;
}
},
factor(n) {
var
r, d;
if
(!Integer.isInteger(n))
throw
TypeError(
"invalid type"
);
r = [];
if
(abs(n) <= 1) {
r.push(n);
return
r;
}
if
(n < 0) {
r.push(-1);
n = -n;
}
while
((n % 2) == 0) {
n >>= 1;
r.push(2);
}
d = 3;
while
(n != 1) {
if
(Integer.isPrime(n)) {
r.push(n);
break
;
}
for
(;;) {
if
((n % d) == 0)
break
;
d += 2;
}
for
(;;) {
r.push(d);
n = Integer.tdiv(n, d);
if
((n % d) != 0)
break
;
}
}
return
r;
},
});
add_props(Integer.prototype, {
inverse() {
return
1 /
this
;
},
norm2() {
return
this
*
this
;
},
abs() {
var
v =
this
;
if
(v < 0)
v = -v;
return
v;
},
conj() {
return
this
;
},
arg() {
if
(
this
>= 0)
return
0;
else
return
Float.PI;
},
exp() {
if
(
this
== 0)
return
1;
else
return
Float.exp(
this
);
},
log() {
if
(
this
== 1)
return
0;
else
return
Float(
this
).log();
},
});
Fraction =
function
Fraction(a, b)
{
var
d, r, obj;
if
(
new
.target)
throw
TypeError(
"not a constructor"
);
if
(a
instanceof
Fraction)
return
a;
if
(!Integer.isInteger(a))
throw
TypeError(
"integer expected"
);
if
(
typeof
b ===
"undefined"
) {
b = 1;
}
else
{
if
(!Integer.isInteger(b))
throw
TypeError(
"integer expected"
);
if
(b == 0)
throw
RangeError(
"division by zero"
);
d = Integer.gcd(a, b);
if
(d != 1) {
a = Integer.tdiv(a, d);
b = Integer.tdiv(b, d);
}
if
(b < 0) {
a = -a;
b = -b;
}
}
obj = Object.create(Fraction.prototype);
obj.num = a;
obj.den = b;
return
obj;
}
function
fraction_add(a, b) {
a = Fraction(a);
b = Fraction(b);
return
Fraction.toFraction(a.num * b.den + a.den * b.num, a.den * b.den);
}
function
fraction_sub(a, b) {
a = Fraction(a);
b = Fraction(b);
return
Fraction.toFraction(a.num * b.den - a.den * b.num, a.den * b.den);
}
function
fraction_mul(a, b) {
a = Fraction(a);
b = Fraction(b);
return
Fraction.toFraction(a.num * b.num, a.den * b.den);
}
function
fraction_div(a, b) {
a = Fraction(a);
b = Fraction(b);
return
Fraction.toFraction(a.num * b.den, a.den * b.num);
}
function
fraction_mod(a, b) {
var
a1 = Fraction(a);
var
b1 = Fraction(b);
return
a - Integer.ediv(a1.num * b1.den, a1.den * b1.num) * b;
}
function
fraction_eq(a, b) {
a = Fraction(a);
b = Fraction(b);
return
(a.num == b.num && a.den == b.den);
}
function
fraction_lt(a, b) {
a = Fraction(a);
b = Fraction(b);
return
(a.num * b.den < b.num * a.den);
}
function
float_add(a, b) {
return
Float(a) + Float(b);
}
function
float_sub(a, b) {
return
Float(a) - Float(b);
}
function
float_mul(a, b) {
return
Float(a) * Float(b);
}
function
float_div(a, b) {
return
Float(a) / Float(b);
}
function
float_mod(a, b) {
return
Float(a) % Float(b);
}
function
float_pow(a, b) {
return
Float(a) ** Float(b);
}
function
float_eq(a, b) {
return
Float(a) === Float(b);
}
function
float_lt(a, b) {
a = Float(a);
b = Float(b);
if
(Float.isNaN(a) || Float.isNaN(b))
return
undefined;
else
return
a < b;
}
operators_set(Fraction.prototype,
{
"+"
: fraction_add,
"-"
: fraction_sub,
"*"
: fraction_mul,
"/"
: fraction_div,
"%"
: fraction_mod,
"**"
: generic_pow,
"=="
: fraction_eq,
"<"
: fraction_lt,
"pos"
(a) {
return
a;
},
"neg"
(a) {
return
Fraction(-a.num, a.den);
},
},
{
left: [Number, BigInt],
right: [Number, BigInt],
"+"
: fraction_add,
"-"
: fraction_sub,
"*"
: fraction_mul,
"/"
: fraction_div,
"%"
: fraction_mod,
"**"
: generic_pow,
"=="
: fraction_eq,
"<"
: fraction_lt,
},
{
left: Float,
right: Float,
"+"
: float_add,
"-"
: float_sub,
"*"
: float_mul,
"/"
: float_div,
"%"
: float_mod,
"**"
: float_pow,
"=="
: float_eq,
"<"
: float_lt,
});
add_props(Fraction, {
toFraction(a, b) {
var
r = Fraction(a, b);
if
(algebraicMode && r.den == 1)
return
r.num;
else
return
r;
},
});
add_props(Fraction.prototype, {
[Symbol.toPrimitive](hint) {
if
(hint ===
"string"
) {
return
this
.toString();
}
else
{
return
Float(
this
.num) /
this
.den;
}
},
inverse() {
return
Fraction(
this
.den,
this
.num);
},
toString() {
return
this
.num +
"/"
+
this
.den;
},
norm2() {
return
this
*
this
;
},
abs() {
if
(
this
.num < 0)
return
-
this
;
else
return
this
;
},
conj() {
return
this
;
},
arg() {
if
(
this
.num >= 0)
return
0;
else
return
Float.PI;
},
exp() {
return
Float.exp(Float(
this
));
},
log() {
return
Float(
this
).log();
},
});
add_props(Number.prototype, {
inverse() {
return
1 /
this
;
},
norm2() {
return
this
*
this
;
},
abs() {
return
Math.abs(
this
);
},
conj() {
return
this
;
},
arg() {
if
(
this
>= 0)
return
0;
else
return
Float.PI;
},
exp() {
return
Float.exp(
this
);
},
log() {
if
(
this
< 0) {
return
Complex(
this
).log();
}
else
{
return
Float.log(
this
);
}
},
});
var
const_tab = [];
function
get_const(n) {
var
t, c, p;
t = const_tab[n];
p = BigFloatEnv.prec;
if
(t && t.prec == p) {
return
t.val;
}
else
{
switch
(n) {
case
0: c = Float.exp(1);
break
;
case
1: c = Float.log(10);
break
;
case
3: c = 1/Float.log(2);
break
;
case
4: c = 1/Float.log(10);
break
;
case
6: c = Float.sqrt(0.5);
break
;
case
7: c = Float.sqrt(2);
break
;
}
if
(p <= 1024) {
const_tab[n] = { prec: p, val: c };
}
return
c;
}
}
add_props(Float, {
isFloat(a) {
return
typeof
a ===
"number"
||
typeof
a ===
"bigfloat"
;
},
bestappr(u, b) {
var
num1, num0, den1, den0, u, num, den, n;
if
(
typeof
b ===
"undefined"
)
throw
TypeError(
"second argument expected"
);
num1 = 1;
num0 = 0;
den1 = 0;
den0 = 1;
for
(;;) {
n = Integer(Float.floor(u));
num = n * num1 + num0;
den = n * den1 + den0;
if
(den > b)
break
;
u = 1.0 / (u - n);
num0 = num1;
num1 = num;
den0 = den1;
den1 = den;
}
return
Fraction(num1, den1);
},
get E() {
return
get_const(0); },
get LN10() {
return
get_const(1); },
get LOG2E() {
return
get_const(3); },
get LOG10E() {
return
get_const(4); },
get SQRT1_2() {
return
get_const(6); },
get SQRT2() {
return
get_const(7); },
});
add_props(Float.prototype, {
inverse() {
return
1.0 /
this
;
},
norm2() {
return
this
*
this
;
},
abs() {
return
Float.abs(
this
);
},
conj() {
return
this
;
},
arg() {
if
(
this
>= 0)
return
0;
else
return
Float.PI;
},
exp() {
return
Float.exp(
this
);
},
log() {
if
(
this
< 0) {
return
Complex(
this
).log();
}
else
{
return
Float.log(
this
);
}
},
});
Complex =
function
Complex(re, im)
{
var
obj;
if
(
new
.target)
throw
TypeError(
"not a constructor"
);
if
(re
instanceof
Complex)
return
re;
if
(
typeof
im ===
"undefined"
) {
im = 0;
}
obj = Object.create(Complex.prototype);
obj.re = re;
obj.im = im;
return
obj;
}
function
complex_add(a, b) {
a = Complex(a);
b = Complex(b);
return
Complex.toComplex(a.re + b.re, a.im + b.im);
}
function
complex_sub(a, b) {
a = Complex(a);
b = Complex(b);
return
Complex.toComplex(a.re - b.re, a.im - b.im);
}
function
complex_mul(a, b) {
a = Complex(a);
b = Complex(b);
return
Complex.toComplex(a.re * b.re - a.im * b.im,
a.re * b.im + a.im * b.re);
}
function
complex_div(a, b) {
a = Complex(a);
b = Complex(b);
return
a * b.inverse();
}
function
complex_eq(a, b) {
a = Complex(a);
b = Complex(b);
return
a.re == b.re && a.im == b.im;
}
operators_set(Complex.prototype,
{
"+"
: complex_add,
"-"
: complex_sub,
"*"
: complex_mul,
"/"
: complex_div,
"**"
: generic_pow,
"=="
: complex_eq,
"pos"
(a) {
return
a;
},
"neg"
(a) {
return
Complex(-a.re, -a.im);
}
},
{
left: [Number, BigInt, Float, Fraction],
right: [Number, BigInt, Float, Fraction],
"+"
: complex_add,
"-"
: complex_sub,
"*"
: complex_mul,
"/"
: complex_div,
"**"
: generic_pow,
"=="
: complex_eq,
});
add_props(Complex, {
toComplex(re, im) {
if
(algebraicMode && im == 0)
return
re;
else
return
Complex(re, im);
},
});
add_props(Complex.prototype, {
inverse() {
var
c =
this
.norm2();
return
Complex(
this
.re / c, -
this
.im / c);
},
toString() {
var
v, s =
""
, a =
this
;
if
(a.re != 0)
s += a.re.toString();
if
(a.im == 1) {
if
(s !=
""
)
s +=
"+"
;
s +=
"I"
;
}
else
if
(a.im == -1) {
s +=
"-I"
;
}
else
{
v = a.im.toString();
if
(v[0] !=
"-"
&& s !=
""
)
s +=
"+"
;
s += v +
"*I"
;
}
return
s;
},
norm2() {
return
this
.re *
this
.re +
this
.im *
this
.im;
},
abs() {
return
Float.sqrt(norm2(
this
));
},
conj() {
return
Complex(
this
.re, -
this
.im);
},
arg() {
return
Float.atan2(
this
.im,
this
.re);
},
exp() {
var
arg =
this
.im, r =
this
.re.exp();
return
Complex(r * cos(arg), r * sin(arg));
},
log() {
return
Complex(abs(
this
).log(), atan2(
this
.im,
this
.re));
},
});
Mod =
function
Mod(a, m) {
var
obj, t;
if
(
new
.target)
throw
TypeError(
"not a constructor"
);
obj = Object.create(Mod.prototype);
if
(Integer.isInteger(m)) {
if
(m <= 0)
throw
RangeError(
"the modulo cannot be <= 0"
);
if
(Integer.isInteger(a)) {
a %= m;
}
else
if
(a
instanceof
Fraction) {
return
Mod(a.num, m) / a.den;
}
else
{
throw
TypeError(
"invalid types"
);
}
}
else
{
throw
TypeError(
"invalid types"
);
}
obj.res = a;
obj.mod = m;
return
obj;
};
function
mod_add(a, b) {
if
(!(a
instanceof
Mod)) {
return
Mod(a + b.res, b.mod);
}
else
if
(!(b
instanceof
Mod)) {
return
Mod(a.res + b, a.mod);
}
else
{
if
(a.mod != b.mod)
throw
TypeError(
"different modulo for binary operator"
);
return
Mod(a.res + b.res, a.mod);
}
}
function
mod_sub(a, b) {
if
(!(a
instanceof
Mod)) {
return
Mod(a - b.res, b.mod);
}
else
if
(!(b
instanceof
Mod)) {
return
Mod(a.res - b, a.mod);
}
else
{
if
(a.mod != b.mod)
throw
TypeError(
"different modulo for binary operator"
);
return
Mod(a.res - b.res, a.mod);
}
}
function
mod_mul(a, b) {
if
(!(a
instanceof
Mod)) {
return
Mod(a * b.res, b.mod);
}
else
if
(!(b
instanceof
Mod)) {
return
Mod(a.res * b, a.mod);
}
else
{
if
(a.mod != b.mod)
throw
TypeError(
"different modulo for binary operator"
);
return
Mod(a.res * b.res, a.mod);
}
}
function
mod_div(a, b) {
if
(!(b
instanceof
Mod))
b = Mod(b, a.mod);
return
mod_mul(a, b.inverse());
}
function
mod_eq(a, b) {
return
(a.mod == b.mod && a.res == b.res);
}
operators_set(Mod.prototype,
{
"+"
: mod_add,
"-"
: mod_sub,
"*"
: mod_mul,
"/"
: mod_div,
"**"
: generic_pow,
"=="
: mod_eq,
"pos"
(a) {
return
a;
},
"neg"
(a) {
return
Mod(-a.res, a.mod);
}
},
{
left: [Number, BigInt, Float, Fraction],
right: [Number, BigInt, Float, Fraction],
"+"
: mod_add,
"-"
: mod_sub,
"*"
: mod_mul,
"/"
: mod_div,
"**"
: generic_pow,
});
add_props(Mod.prototype, {
inverse() {
var
a =
this
, m = a.mod;
if
(Integer.isInteger(m)) {
return
Mod(Integer.invmod(a.res, m), m);
}
else
{
throw
TypeError(
"unsupported type"
);
}
},
toString() {
return
"Mod("
+
this
.res +
","
+
this
.mod +
")"
;
},
});
function
polynomial_is_scalar(a)
{
if
(
typeof
a ===
"number"
||
typeof
a ===
"bigint"
||
typeof
a ===
"bigfloat"
)
return
true
;
if
(a
instanceof
Fraction ||
a
instanceof
Complex ||
a
instanceof
Mod)
return
true
;
return
false
;
}
Polynomial =
function
Polynomial(a)
{
if
(
new
.target)
throw
TypeError(
"not a constructor"
);
if
(a
instanceof
Polynomial) {
return
a;
}
else
if
(Array.isArray(a)) {
if
(a.length == 0)
a = [ 0 ];
Object.setPrototypeOf(a, Polynomial.prototype);
return
a.trim();
}
else
if
(polynomial_is_scalar(a)) {
a = [a];
Object.setPrototypeOf(a, Polynomial.prototype);
return
a;
}
else
{
throw
TypeError(
"invalid type"
);
}
}
function
number_need_paren(c)
{
return
!(Integer.isInteger(c) ||
Float.isFloat(c) ||
c
instanceof
Fraction ||
(c
instanceof
Complex && c.re == 0));
}
function
monomial_toString(c, i)
{
var
str1;
if
(i == 0) {
str1 = c.toString();
}
else
{
if
(c == 1) {
str1 =
""
;
}
else
if
(c == -1) {
str1 =
"-"
;
}
else
{
if
(number_need_paren(c)) {
str1 =
"("
+ c +
")"
;
}
else
{
str1 = String(c);
}
str1 +=
"*"
;
}
str1 +=
"X"
;
if
(i != 1) {
str1 +=
"^"
+ i;
}
}
return
str1;
}
function
poly_root_laguerre1(p, z, max_it)
{
var
p1, p2, i, z0, z1, z2, d, t0, t1, d1, d2, e, el, zl;
d = p.deg();
if
(d == 1) {
return
-p[0] / p[1];
}
if
(p[0] == 0)
return
0.0;
p1 = p.deriv();
p2 = p1.deriv();
el = 0.0;
zl = 0.0;
for
(i = 0; i < max_it; i++) {
z0 = p.apply(z);
if
(z0 == 0)
return
z;
e = abs(z - zl);
if
(i >= 2 && e >= el) {
if
(abs(zl) < 1e-4) {
if
(e < 1e-7)
return
zl;
}
else
{
if
(e < abs(zl) * 1e-3)
return
zl;
}
}
el = e;
zl = z;
z1 = p1.apply(z);
z2 = p2.apply(z);
t0 = (d - 1) * z1;
t0 = t0 * t0;
t1 = d * (d - 1) * z0 * z2;
t0 = sqrt(t0 - t1);
d1 = z1 + t0;
d2 = z1 - t0;
if
(norm2(d2) > norm2(d1))
d1 = d2;
if
(d1 == 0)
return
null
;
z = z - d * z0 / d1;
}
return
null
;
}
function
poly_roots(p)
{
var
d, i, roots, j, z, eps;
var
start_points = [ 0.1, -1.4, 1.7 ];
if
(!(p
instanceof
Polynomial))
throw
TypeError(
"polynomial expected"
);
d = p.deg();
if
(d <= 0)
return
[];
eps = 2.0 ^ (-BigFloatEnv.prec);
roots = [];
for
(i = 0; i < d; i++) {
for
(j = 0; j < 3; j++) {
z = poly_root_laguerre1(p, start_points[j], 100);
if
(z !==
null
)
break
;
}
if
(j == 3)
throw
RangeError(
"error in root finding algorithm"
);
roots[i] = z;
p = Polynomial.divrem(p, X - z)[0];
}
return
roots;
}
add_props(Polynomial.prototype, {
trim() {
var
a =
this
, i;
i = a.length;
while
(i > 1 && a[i - 1] == 0)
i--;
a.length = i;
return
a;
},
conj() {
var
r, i, n, a;
a =
this
;
n = a.length;
r = [];
for
(i = 0; i < n; i++)
r[i] = a[i].conj();
return
Polynomial(r);
},
inverse() {
return
RationalFunction(Polynomial([1]),
this
);
},
toString() {
var
i, str, str1, c, a =
this
;
if
(a.length == 1) {
return
a[0].toString();
}
str=
""
;
for
(i = a.length - 1; i >= 0; i--) {
c = a[i];
if
(c == 0 ||
(c
instanceof
Mod) && c.res == 0)
continue
;
str1 = monomial_toString(c, i);
if
(str1[0] !=
"-"
) {
if
(str !=
""
)
str +=
"+"
;
}
str += str1;
}
return
str;
},
deg() {
if
(
this
.length == 1 &&
this
[0] == 0)
return
-Infinity;
else
return
this
.length - 1;
},
apply(b) {
var
i, n, r, a =
this
;
n = a.length - 1;
r = a[n];
while
(n > 0) {
n--;
r = r * b + a[n];
}
return
r;
},
deriv() {
var
a =
this
, n, r, i;
n = a.length;
if
(n == 1) {
return
Polynomial(0);
}
else
{
r = [];
for
(i = 1; i < n; i++) {
r[i - 1] = i * a[i];
}
return
Polynomial(r);
}
},
integ() {
var
a =
this
, n, r, i;
n = a.length;
r = [0];
for
(i = 0; i < n; i++) {
r[i + 1] = a[i] / (i + 1);
}
return
Polynomial(r);
},
norm2() {
var
a =
this
, n, r, i;
n = a.length;
r = 0;
for
(i = 0; i < n; i++) {
r += a[i].norm2();
}
return
r;
},
});
function
polynomial_add(a, b) {
var
tmp, r, i, n1, n2;
a = Polynomial(a);
b = Polynomial(b);
if
(a.length < b.length) {
tmp = a;
a = b;
b = tmp;
}
n1 = b.length;
n2 = a.length;
r = [];
for
(i = 0; i < n1; i++)
r[i] = a[i] + b[i];
for
(i = n1; i < n2; i++)
r[i] = a[i];
return
Polynomial(r);
}
function
polynomial_sub(a, b) {
return
polynomial_add(a, -b);
}
function
polynomial_mul(a, b) {
var
i, j, n1, n2, n, r;
a = Polynomial(a);
b = Polynomial(b);
n1 = a.length;
n2 = b.length;
n = n1 + n2 - 1;
r = [];
for
(i = 0; i < n; i++)
r[i] = 0;
for
(i = 0; i < n1; i++) {
for
(j = 0; j < n2; j++) {
r[i + j] += a[i] * b[j];
}
}
return
Polynomial(r);
}
function
polynomial_div_scalar(a, b) {
return
a * (1 / b);
}
function
polynomial_div(a, b)
{
return
RationalFunction(Polynomial(a),
Polynomial(b));
}
function
polynomial_mod(a, b) {
return
Polynomial.divrem(a, b)[1];
}
function
polynomial_eq(a, b) {
var
n, i;
n = a.length;
if
(n != b.length)
return
false
;
for
(i = 0; i < n; i++) {
if
(a[i] != b[i])
return
false
;
}
return
true
;
}
operators_set(Polynomial.prototype,
{
"+"
: polynomial_add,
"-"
: polynomial_sub,
"*"
: polynomial_mul,
"/"
: polynomial_div,
"**"
: generic_pow,
"=="
: polynomial_eq,
"pos"
(a) {
return
a;
},
"neg"
(a) {
var
r, i, n, a;
n = a.length;
r = [];
for
(i = 0; i < n; i++)
r[i] = -a[i];
return
Polynomial(r);
},
},
{
left: [Number, BigInt, Float, Fraction, Complex, Mod],
"+"
: polynomial_add,
"-"
: polynomial_sub,
"*"
: polynomial_mul,
"/"
: polynomial_div,
"**"
: generic_pow,
},
{
right: [Number, BigInt, Float, Fraction, Complex, Mod],
"+"
: polynomial_add,
"-"
: polynomial_sub,
"*"
: polynomial_mul,
"/"
: polynomial_div_scalar,
"**"
: generic_pow,
});
add_props(Polynomial, {
divrem(a, b) {
var
n1, n2, i, j, q, r, n, c;
if
(b.deg() < 0)
throw
RangeError(
"division by zero"
);
n1 = a.length;
n2 = b.length;
if
(n1 < n2)
return
[Polynomial([0]), a];
r = Array.prototype.dup.call(a);
q = [];
n2--;
n = n1 - n2;
for
(i = 0; i < n; i++)
q[i] = 0;
for
(i = n - 1; i >= 0; i--) {
c = r[i + n2];
if
(c != 0) {
c = c / b[n2];
r[i + n2] = 0;
for
(j = 0; j < n2; j++) {
r[i + j] -= b[j] * c;
}
q[i] = c;
}
}
return
[Polynomial(q), Polynomial(r)];
},
gcd(a, b) {
var
t;
while
(b.deg() >= 0) {
t = Polynomial.divrem(a, b);
a = b;
b = t[1];
}
return
a / a[a.length - 1];
},
invmod(x, y) {
var
q, u, v, a, c, t;
u = x;
v = y;
c = Polynomial([1]);
a = Polynomial([0]);
while
(u.deg() >= 0) {
t = Polynomial.divrem(v, u);
q = t[0];
v = u;
u = t[1];
t = c;
c = a - q * c;
a = t;
}
if
(v.deg() > 0)
throw
RangeError(
"not invertible"
);
return
Polynomial.divrem(a, y)[1];
},
roots(p) {
return
poly_roots(p);
}
});
PolyMod =
function
PolyMod(a, m) {
var
obj, t;
if
(
new
.target)
throw
TypeError(
"not a constructor"
);
obj = Object.create(PolyMod.prototype);
if
(m
instanceof
Polynomial) {
if
(m.deg() <= 0)
throw
RangeError(
"the modulo cannot have a degree <= 0"
);
if
(a
instanceof
RationalFunction) {
return
PolyMod(a.num, m) / a.den;
}
else
{
a = Polynomial(a);
t = Polynomial.divrem(a, m);
a = t[1];
}
}
else
{
throw
TypeError(
"invalid types"
);
}
obj.res = a;
obj.mod = m;
return
obj;
};
function
polymod_add(a, b) {
if
(!(a
instanceof
PolyMod)) {
return
PolyMod(a + b.res, b.mod);
}
else
if
(!(b
instanceof
PolyMod)) {
return
PolyMod(a.res + b, a.mod);
}
else
{
if
(a.mod != b.mod)
throw
TypeError(
"different modulo for binary operator"
);
return
PolyMod(a.res + b.res, a.mod);
}
}
function
polymod_sub(a, b) {
return
polymod_add(a, -b);
}
function
polymod_mul(a, b) {
if
(!(a
instanceof
PolyMod)) {
return
PolyMod(a * b.res, b.mod);
}
else
if
(!(b
instanceof
PolyMod)) {
return
PolyMod(a.res * b, a.mod);
}
else
{
if
(a.mod != b.mod)
throw
TypeError(
"different modulo for binary operator"
);
return
PolyMod(a.res * b.res, a.mod);
}
}
function
polymod_div(a, b) {
if
(!(b
instanceof
PolyMod))
b = PolyMod(b, a.mod);
return
polymod_mul(a, b.inverse());
}
function
polymod_eq(a, b) {
return
(a.mod == b.mod && a.res == b.res);
}
operators_set(PolyMod.prototype,
{
"+"
: polymod_add,
"-"
: polymod_sub,
"*"
: polymod_mul,
"/"
: polymod_div,
"**"
: generic_pow,
"=="
: polymod_eq,
"pos"
(a) {
return
a;
},
"neg"
(a) {
return
PolyMod(-a.res, a.mod);
},
},
{
left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
"+"
: polymod_add,
"-"
: polymod_sub,
"*"
: polymod_mul,
"/"
: polymod_div,
"**"
: generic_pow,
});
add_props(PolyMod.prototype, {
inverse() {
var
a =
this
, m = a.mod;
if
(m
instanceof
Polynomial) {
return
PolyMod(Polynomial.invmod(a.res, m), m);
}
else
{
throw
TypeError(
"unsupported type"
);
}
},
toString() {
return
"PolyMod("
+
this
.res +
","
+
this
.mod +
")"
;
},
});
RationalFunction =
function
RationalFunction(a, b)
{
var
t, r, d, obj;
if
(
new
.target)
throw
TypeError(
"not a constructor"
);
if
(!(a
instanceof
Polynomial) ||
!(b
instanceof
Polynomial))
throw
TypeError(
"polynomial expected"
);
t = Polynomial.divrem(a, b);
r = t[1];
if
(r.deg() < 0)
return
t[0];
d = Polynomial.gcd(b, r);
if
(d.deg() > 0) {
a = Polynomial.divrem(a, d)[0];
b = Polynomial.divrem(b, d)[0];
}
obj = Object.create(RationalFunction.prototype);
obj.num = a;
obj.den = b;
return
obj;
}
add_props(RationalFunction.prototype, {
inverse() {
return
RationalFunction(
this
.den,
this
.num);
},
conj() {
return
RationalFunction(
this
.num.conj(),
this
.den.conj());
},
toString() {
var
str;
if
(
this
.num.deg() <= 0 &&
!number_need_paren(
this
.num[0]))
str =
this
.num.toString();
else
str =
"("
+
this
.num.toString() +
")"
;
str +=
"/("
+
this
.den.toString() +
")"
return
str;
},
apply(b) {
return
this
.num.apply(b) /
this
.den.apply(b);
},
deriv() {
var
n =
this
.num, d =
this
.den;
return
RationalFunction(n.deriv() * d - n * d.deriv(), d * d);
},
});
function
ratfunc_add(a, b) {
a = RationalFunction.toRationalFunction(a);
b = RationalFunction.toRationalFunction(b);
return
RationalFunction(a.num * b.den + a.den * b.num, a.den * b.den);
}
function
ratfunc_sub(a, b) {
a = RationalFunction.toRationalFunction(a);
b = RationalFunction.toRationalFunction(b);
return
RationalFunction(a.num * b.den - a.den * b.num, a.den * b.den);
}
function
ratfunc_mul(a, b) {
a = RationalFunction.toRationalFunction(a);
b = RationalFunction.toRationalFunction(b);
return
RationalFunction(a.num * b.num, a.den * b.den);
}
function
ratfunc_div(a, b) {
a = RationalFunction.toRationalFunction(a);
b = RationalFunction.toRationalFunction(b);
return
RationalFunction(a.num * b.den, a.den * b.num);
}
function
ratfunc_eq(a, b) {
a = RationalFunction.toRationalFunction(a);
b = RationalFunction.toRationalFunction(b);
return
(a.num == b.num && a.den == b.den);
}
operators_set(RationalFunction.prototype,
{
"+"
: ratfunc_add,
"-"
: ratfunc_sub,
"*"
: ratfunc_mul,
"/"
: ratfunc_div,
"**"
: generic_pow,
"=="
: ratfunc_eq,
"pos"
(a) {
return
a;
},
"neg"
(a) {
return
RationalFunction(-
this
.num,
this
.den);
},
},
{
left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
"+"
: ratfunc_add,
"-"
: ratfunc_sub,
"*"
: ratfunc_mul,
"/"
: ratfunc_div,
"**"
: generic_pow,
});
add_props(RationalFunction, {
toRationalFunction(a) {
var
obj;
if
(a
instanceof
RationalFunction) {
return
a;
}
else
{
obj = Object.create(RationalFunction.prototype);
obj.num = Polynomial(a);
obj.den = Polynomial(1);
return
obj;
}
},
});
function
get_emin(a) {
var
i, n;
n = a.length;
for
(i = 0; i < n; i++) {
if
(a[i] != 0)
return
i;
}
return
n;
};
function
series_is_scalar_or_polynomial(a)
{
return
polynomial_is_scalar(a) ||
(a
instanceof
Polynomial);
}
Series =
function
Series(a, n) {
var
emin, r, i;
if
(a
instanceof
Series) {
return
a;
}
else
if
(series_is_scalar_or_polynomial(a)) {
if
(n <= 0) {
return
Series.zero(0, 0);
}
else
{
a = Polynomial(a);
emin = get_emin(a);
r = Series.zero(n, emin);
n = Math.min(a.length - emin, n);
for
(i = 0; i < n; i++)
r[i] = a[i + emin];
return
r;
}
}
else
if
(a
instanceof
RationalFunction) {
return
Series(a.num, n) / a.den;
}
else
{
throw
TypeError(
"invalid type"
);
}
};
function
series_add(v1, v2) {
var
tmp, d, emin, n, r, i, j, v2_emin, c1, c2;
if
(!(v1
instanceof
Series)) {
tmp = v1;
v1 = v2;
v2 = tmp;
}
d = v1.emin + v1.length;
if
(series_is_scalar_or_polynomial(v2)) {
v2 = Polynomial(v2);
if
(d <= 0)
return
v1;
v2_emin = 0;
}
else
if
(v2
instanceof
RationalFunction) {
i = get_emin(v2.num) - get_emin(v2.den);
if
(d <= i)
return
v1;
v2 = Series(v2, d - i);
v2_emin = v2.emin;
}
else
{
v2_emin = v2.emin;
d = Math.min(d, v2_emin + v2.length);
}
emin = Math.min(v1.emin, v2_emin);
n = d - emin;
r = Series.zero(n, emin);
for
(i = emin; i < d; i++) {
j = i - v1.emin;
if
(j >= 0 && j < v1.length)
c1 = v1[j];
else
c1 = 0;
j = i - v2_emin;
if
(j >= 0 && j < v2.length)
c2 = v2[j];
else
c2 = 0;
r[i - emin] = c1 + c2;
}
return
r.trim();
}
function
series_sub(a, b) {
return
series_add(a, -b);
}
function
series_mul(v1, v2) {
var
n, i, j, r, n, emin, n1, n2, k;
if
(!(v1
instanceof
Series))
v1 = Series(v1, v2.length);
else
if
(!(v2
instanceof
Series))
v2 = Series(v2, v1.length);
emin = v1.emin + v2.emin;
n = Math.min(v1.length, v2.length);
n1 = v1.length;
n2 = v2.length;
r = Series.zero(n, emin);
for
(i = 0; i < n1; i++) {
k = Math.min(n2, n - i);
for
(j = 0; j < k; j++) {
r[i + j] += v1[i] * v2[j];
}
}
return
r.trim();
}
function
series_div(v1, v2) {
if
(!(v2
instanceof
Series))
v2 = Series(v2, v1.length);
return
series_mul(v1, v2.inverse());
}
function
series_pow(a, b) {
if
(Integer.isInteger(b)) {
return
generic_pow(a, b);
}
else
{
if
(!(a
instanceof
Series))
a = Series(a, b.length);
return
exp(log(a) * b);
}
}
function
series_eq(a, b) {
var
n, i;
if
(a.emin != b.emin)
return
false
;
n = a.length;
if
(n != b.length)
return
false
;
for
(i = 0; i < n; i++) {
if
(a[i] != b[i])
return
false
;
}
return
true
;
}
operators_set(Series.prototype,
{
"+"
: series_add,
"-"
: series_sub,
"*"
: series_mul,
"/"
: series_div,
"**"
: series_pow,
"=="
: series_eq,
"pos"
(a) {
return
a;
},
"neg"
(a) {
var
obj, n, i;
n = a.length;
obj = Series.zero(a.length, a.emin);
for
(i = 0; i < n; i++) {
obj[i] = -a[i];
}
return
obj;
},
},
{
left: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
right: [Number, BigInt, Float, Fraction, Complex, Mod, Polynomial],
"+"
: series_add,
"-"
: series_sub,
"*"
: series_mul,
"/"
: series_div,
"**"
: series_pow,
});
add_props(Series.prototype, {
conj() {
var
obj, n, i;
n =
this
.length;
obj = Series.zero(
this
.length,
this
.emin);
for
(i = 0; i < n; i++) {
obj[i] =
this
[i].conj();
}
return
obj;
},
inverse() {
var
r, n, i, j, sum, v1 =
this
;
n = v1.length;
if
(n == 0)
throw
RangeError(
"division by zero"
);
r = Series.zero(n, -v1.emin);
r[0] = 1 / v1[0];
for
(i = 1; i < n; i++) {
sum = 0;
for
(j = 1; j <= i; j++) {
sum += v1[j] * r[i - j];
}
r[i] = -sum * r[0];
}
return
r;
},
trim() {
var
i, j, n, r, v1 =
this
;
n = v1.length;
i = 0;
while
(i < n && v1[i] == 0)
i++;
if
(i == 0)
return
v1;
for
(j = i; j < n; j++)
v1[j - i] = v1[j];
v1.length = n - i;
v1.__proto__.emin += i;
return
v1;
},
toString() {
var
i, j, str, str1, c, a =
this
, emin, n;
str=
""
;
emin =
this
.emin;
n =
this
.length;
for
(j = 0; j < n; j++) {
i = j + emin;
c = a[j];
if
(c != 0) {
str1 = monomial_toString(c, i);
if
(str1[0] !=
"-"
) {
if
(str !=
""
)
str +=
"+"
;
}
str += str1;
}
}
if
(str !=
""
)
str +=
"+"
;
str +=
"O("
+ monomial_toString(1, n + emin) +
")"
;
return
str;
},
apply(b) {
var
i, n, r, a =
this
;
n = a.length;
if
(n == 0)
return
0;
r = a[--n];
while
(n > 0) {
n--;
r = r * b + a[n];
}
if
(a.emin != 0)
r *= b ^ a.emin;
return
r;
},
deriv() {
var
a =
this
, n = a.length, emin = a.emin, r, i, j;
if
(n == 0 && emin == 0) {
return
Series.zero(0, 0);
}
else
{
r = Series.zero(n, emin - 1);
for
(i = 0; i < n; i++) {
j = emin + i;
if
(j == 0)
r[i] = 0;
else
r[i] = j * a[i];
}
return
r.trim();
}
},
integ() {
var
a =
this
, n = a.length, emin = a.emin, i, j, r;
r = Series.zero(n, emin + 1);
for
(i = 0; i < n; i++) {
j = emin + i;
if
(j == -1) {
if
(a[i] != 0)
throw
RangeError(
"cannot represent integ(1/X)"
);
}
else
{
r[i] = a[i] / (j + 1);
}
}
return
r.trim();
},
exp() {
var
c, i, r, n, a =
this
;
if
(a.emin < 0)
throw
RangeError(
"negative exponent in exp"
);
n = a.emin + a.length;
if
(a.emin > 0 || a[0] == 0) {
c = 1;
}
else
{
c = global.exp(a[0]);
a -= a[0];
}
r = Series.zero(n, 0);
for
(i = 0; i < n; i++) {
r[i] = c / fact(i);
}
return
r.apply(a);
},
log() {
var
a =
this
, r;
if
(a.emin != 0)
throw
RangeError(
"log argument must have a non zero constant term"
);
r = integ(deriv(a) / a);
r += global.log(a[0]);
return
r;
},
});
add_props(Series, {
zero(n, emin) {
var
r, i, obj;
r = [];
for
(i = 0; i < n; i++)
r[i] = 0;
obj = Object.create(Series.prototype);
obj.emin = emin;
Object.setPrototypeOf(r, obj);
return
r;
},
O(a) {
function
ErrorO() {
return
TypeError(
"invalid O() argument"
);
}
var
n;
if
(series_is_scalar_or_polynomial(a)) {
a = Polynomial(a);
n = a.deg();
if
(n < 0)
throw
ErrorO();
}
else
if
(a
instanceof
RationalFunction) {
if
(a.num.deg() != 0)
throw
ErrorO();
n = a.den.deg();
if
(n < 0)
throw
ErrorO();
n = -n;
}
else
throw
ErrorO();
return
Series.zero(0, n);
},
});
Matrix =
function
Matrix(h, w) {
var
i, j, r, rl;
if
(
typeof
w ===
"undefined"
)
w = h;
r = [];
for
(i = 0; i < h; i++) {
rl = [];
for
(j = 0; j < w; j++)
rl[j] = 0;
r[i] = rl;
}
return
r;
};
add_props(Matrix, {
idn(n) {
var
r, i;
r = Matrix(n, n);
for
(i = 0; i < n; i++)
r[i][i] = 1;
return
r;
},
diag(a) {
var
r, i, n;
n = a.length;
r = Matrix(n, n);
for
(i = 0; i < n; i++)
r[i][i] = a[i];
return
r;
},
hilbert(n) {
var
i, j, r;
r = Matrix(n);
for
(i = 0; i < n; i++) {
for
(j = 0; j < n; j++) {
r[i][j] = 1 / (1 + i + j);
}
}
return
r;
},
trans(a) {
var
h, w, r, i, j;
if
(!Array.isArray(a))
throw
TypeError(
"matrix expected"
);
h = a.length;
if
(!Array.isArray(a[0])) {
w = 1;
r = Matrix(w, h);
for
(i = 0; i < h; i++) {
r[0][i] = a[i];
}
}
else
{
w = a[0].length;
r = Matrix(w, h);
for
(i = 0; i < h; i++) {
for
(j = 0; j < w; j++) {
r[j][i] = a[i][j];
}
}
}
return
r;
},
check_square(a) {
var
a, n;
if
(!Array.isArray(a))
throw
TypeError(
"array expected"
);
n = a.length;
if
(!Array.isArray(a[0]) || n != a[0].length)
throw
TypeError(
"square matrix expected"
);
return
n;
},
trace(a) {
var
n, r, i;
n = Matrix.check_square(a);
r = a[0][0];
for
(i = 1; i < n; i++) {
r += a[i][i];
}
return
r;
},
charpoly(a) {
var
n, p, c, i, j, coef;
n = Matrix.check_square(a);
p = [];
for
(i = 0; i < n + 1; i++)
p[i] = 0;
p[n] = 1;
c = Matrix.idn(n);
for
(i = 0; i < n; i++) {
c = c * a;
coef = -trace(c) / (i + 1);
p[n - i - 1] = coef;
for
(j = 0; j < n; j++)
c[j][j] += coef;
}
return
Polynomial(p);
},
eigenvals(a) {
return
Polynomial.roots(Matrix.charpoly(a));
},
det(a) {
var
n, i, j, k, s, src, v, c;
n = Matrix.check_square(a);
s = 1;
src = a.dup();
for
(i=0;i<n;i++) {
for
(j = i; j < n; j++) {
if
(src[j][i] != 0)
break
;
}
if
(j == n)
return
0;
if
(j != i) {
for
(k = 0;k < n; k++) {
v = src[j][k];
src[j][k] = src[i][k];
src[i][k] = v;
}
s = -s;
}
c = src[i][i].inverse();
for
(j = i + 1; j < n; j++) {
v = c * src[j][i];
for
(k = 0;k < n; k++) {
src[j][k] -= src[i][k] * v;
}
}
}
c = s;
for
(i=0;i<n;i++)
c *= src[i][i];
return
c;
},
inverse(a) {
var
n, dst, src, i, j, k, n2, r, c, v;
n = Matrix.check_square(a);
src = a.dup();
dst = Matrix.idn(n);
for
(i=0;i<n;i++) {
for
(j = i; j < n; j++) {
if
(src[j][i] != 0)
break
;
}
if
(j == n)
throw
RangeError(
"matrix is not invertible"
);
if
(j != i) {
v = src[j];
src[j] = src[i];
src[i] = v;
v = dst[j];
dst[j] = dst[i];
dst[i] = v;
}
c = src[i][i].inverse();
for
(k = 0; k < n; k++) {
src[i][k] *= c;
dst[i][k] *= c;
}
for
(j = 0; j < n; j++) {
if
(j != i) {
c = src[j][i];
for
(k = i; k < n; k++) {
src[j][k] -= src[i][k] * c;
}
for
(k = 0; k < n; k++) {
dst[j][k] -= dst[i][k] * c;
}
}
}
}
return
dst;
},
rank(a) {
var
src, i, j, k, w, h, l, c;
if
(!Array.isArray(a) ||
!Array.isArray(a[0]))
throw
TypeError(
"matrix expected"
);
h = a.length;
w = a[0].length;
src = a.dup();
l = 0;
for
(i=0;i<w;i++) {
for
(j = l; j < h; j++) {
if
(src[j][i] != 0)
break
;
}
if
(j == h)
continue
;
if
(j != l) {
for
(k = 0; k < w; k++) {
v = src[j][k];
src[j][k] = src[l][k];
src[l][k] = v;
}
}
c = src[l][i].inverse();
for
(k = 0; k < w; k++) {
src[l][k] *= c;
}
for
(j = l + 1; j < h; j++) {
c = src[j][i];
for
(k = i; k < w; k++) {
src[j][k] -= src[l][k] * c;
}
}
l++;
}
return
l;
},
ker(a) {
var
src, i, j, k, w, h, l, m, r, im_cols, ker_dim, c;
if
(!Array.isArray(a) ||
!Array.isArray(a[0]))
throw
TypeError(
"matrix expected"
);
h = a.length;
w = a[0].length;
src = a.dup();
im_cols = [];
l = 0;
for
(i=0;i<w;i++) {
im_cols[i] =
false
;
for
(j = l; j < h; j++) {
if
(src[j][i] != 0)
break
;
}
if
(j == h)
continue
;
im_cols[i] =
true
;
if
(j != l) {
for
(k = 0; k < w; k++) {
v = src[j][k];
src[j][k] = src[l][k];
src[l][k] = v;
}
}
c = src[l][i].inverse();
for
(k = 0; k < w; k++) {
src[l][k] *= c;
}
for
(j = 0; j < h; j++) {
if
(j != l) {
c = src[j][i];
for
(k = i; k < w; k++) {
src[j][k] -= src[l][k] * c;
}
}
}
l++;
}
ker_dim = w - l;
r = Matrix(w, ker_dim);
k = 0;
for
(i = 0; i < w; i++) {
if
(!im_cols[i]) {
l = 0;
m = 0;
for
(j = 0; j < w; j++) {
if
(im_cols[j]) {
r[j][k] = -src[m][i];
m++;
}
else
{
if
(l == k) {
r[j][k] = 1;
}
else
{
r[j][k] = 0;
}
l++;
}
}
k++;
}
}
return
r;
},
dp(a, b) {
var
i, n, r;
n = a.length;
if
(n != b.length)
throw
TypeError(
"incompatible array length"
);
r = 0;
for
(i = 0; i < n; i++) {
r += a[i] * b[i];
}
return
r;
},
cp(v1, v2) {
var
r;
if
(v1.length != 3 || v2.length != 3)
throw
TypeError(
"vectors must have 3 elements"
);
r = [];
r[0] = v1[1] * v2[2] - v1[2] * v2[1];
r[1] = v1[2] * v2[0] - v1[0] * v2[2];
r[2] = v1[0] * v2[1] - v1[1] * v2[0];
return
r;
},
});
function
array_add(a, b) {
var
r, i, n;
n = a.length;
if
(n != b.length)
throw
TypeError(
"incompatible array size"
);
r = [];
for
(i = 0; i < n; i++)
r[i] = a[i] + b[i];
return
r;
}
function
array_sub(a, b) {
var
r, i, n;
n = a.length;
if
(n != b.length)
throw
TypeError(
"incompatible array size"
);
r = [];
for
(i = 0; i < n; i++)
r[i] = a[i] - b[i];
return
r;
}
function
array_scalar_mul(a, b) {
var
r, i, n;
n = a.length;
r = [];
for
(i = 0; i < n; i++)
r[i] = a[i] * b;
return
r;
}
function
array_mul(a, b) {
var
h, w, l, i, j, k, r, rl, sum, a_mat, b_mat;
h = a.length;
a_mat = Array.isArray(a[0]);
if
(a_mat) {
l = a[0].length;
}
else
{
l = 1;
}
if
(l != b.length)
throw
RangeError(
"incompatible matrix size"
);
b_mat = Array.isArray(b[0]);
if
(b_mat)
w = b[0].length;
else
w = 1;
r = [];
if
(a_mat && b_mat) {
for
(i = 0; i < h; i++) {
rl = [];
for
(j = 0; j < w; j++) {
sum = 0;
for
(k = 0; k < l; k++) {
sum += a[i][k] * b[k][j];
}
rl[j] = sum;
}
r[i] = rl;
}
}
else
if
(a_mat && !b_mat) {
for
(i = 0; i < h; i++) {
sum = 0;
for
(k = 0; k < l; k++) {
sum += a[i][k] * b[k];
}
r[i] = sum;
}
}
else
if
(!a_mat && b_mat) {
for
(i = 0; i < h; i++) {
rl = [];
for
(j = 0; j < w; j++) {
rl[j] = a[i] * b[0][j];
}
r[i] = rl;
}
}
else
{
for
(i = 0; i < h; i++) {
r[i] = a[i] * b[0];
}
}
return
r;
}
function
array_div(a, b) {
return
array_mul(a, b.inverse());
}
function
array_element_wise_inverse(a) {
var
r, i, n;
n = a.length;
r = [];
for
(i = 0; i < n; i++)
r[i] = a[i].inverse();
return
r;
}
function
array_eq(a, b) {
var
n, i;
n = a.length;
if
(n != b.length)
return
false
;
for
(i = 0; i < n; i++) {
if
(a[i] != b[i])
return
false
;
}
return
true
;
}
operators_set(Array.prototype,
{
"+"
: array_add,
"-"
: array_sub,
"*"
: array_mul,
"/"
: array_div,
"=="
: array_eq,
"pos"
(a) {
return
a;
},
"neg"
(a) {
var
i, n, r;
n = a.length;
r = [];
for
(i = 0; i < n; i++)
r[i] = -a[i];
return
r;
}
},
{
right: [Number, BigInt, Float, Fraction, Complex, Mod,
Polynomial, PolyMod, RationalFunction, Series],
"*"
: array_scalar_mul,
"/"
(a, b) {
return
a * b.inverse(); },
"**"
: generic_pow,
},
{
left: [Number, BigInt, Float, Fraction, Complex, Mod,
Polynomial, PolyMod, RationalFunction, Series],
"*"
(a, b) {
return
array_scalar_mul(b, a); },
"/"
(a, b) {
return
a * array_element_wise_inverse(b); },
});
add_props(Array.prototype, {
conj() {
var
i, n, r;
n =
this
.length;
r = [];
for
(i = 0; i < n; i++)
r[i] =
this
[i].conj();
return
r;
},
dup() {
var
r, i, n, el, a =
this
;
r = [];
n = a.length;
for
(i = 0; i < n; i++) {
el = a[i];
if
(Array.isArray(el))
el = el.dup();
r[i] = el;
}
return
r;
},
inverse() {
return
Matrix.inverse(
this
);
},
norm2: Polynomial.prototype.norm2,
});
})(
this
);
var
I = Complex(0, 1);
var
X = Polynomial([0, 1]);
var
O = Series.O;
Object.defineProperty(
this
,
"PI"
, { get:
function
() {
return
Float.PI } });
var
gcd = Integer.gcd;
var
fact = Integer.fact;
var
comb = Integer.comb;
var
pmod = Integer.pmod;
var
invmod = Integer.invmod;
var
factor = Integer.factor;
var
isprime = Integer.isPrime;
var
nextprime = Integer.nextPrime;
function
deriv(a)
{
return
a.deriv();
}
function
integ(a)
{
return
a.integ();
}
function
norm2(a)
{
return
a.norm2();
}
function
abs(a)
{
return
a.abs();
}
function
conj(a)
{
return
a.conj();
}
function
arg(a)
{
return
a.arg();
}
function
inverse(a)
{
return
a.inverse();
}
function
trunc(a)
{
if
(Integer.isInteger(a)) {
return
a;
}
else
if
(a
instanceof
Fraction) {
return
Integer.tdiv(a.num, a.den);
}
else
if
(a
instanceof
Polynomial) {
return
a;
}
else
if
(a
instanceof
RationalFunction) {
return
Polynomial.divrem(a.num, a.den)[0];
}
else
{
return
Float.ceil(a);
}
}
function
floor(a)
{
if
(Integer.isInteger(a)) {
return
a;
}
else
if
(a
instanceof
Fraction) {
return
Integer.fdiv(a.num, a.den);
}
else
{
return
Float.floor(a);
}
}
function
ceil(a)
{
if
(Integer.isInteger(a)) {
return
a;
}
else
if
(a
instanceof
Fraction) {
return
Integer.cdiv(a.num, a.den);
}
else
{
return
Float.ceil(a);
}
}
function
sqrt(a)
{
var
t, u, re, im;
if
(a
instanceof
Series) {
return
a ^ (1/2);
}
else
if
(a
instanceof
Complex) {
t = abs(a);
u = a.re;
re = sqrt((t + u) / 2);
im = sqrt((t - u) / 2);
if
(a.im < 0)
im = -im;
return
Complex.toComplex(re, im);
}
else
{
a = Float(a);
if
(a < 0) {
return
Complex(0, Float.sqrt(-a));
}
else
{
return
Float.sqrt(a);
}
}
}
function
exp(a)
{
return
a.exp();
}
function
log(a)
{
return
a.log();
}
function
log2(a)
{
return
log(a) * Float.LOG2E;
}
function
log10(a)
{
return
log(a) * Float.LOG10E;
}
function
todb(a)
{
return
log10(a) * 10;
}
function
fromdb(a)
{
return
10 ^ (a / 10);
}
function
sin(a)
{
var
t;
if
(a
instanceof
Complex || a
instanceof
Series) {
t = exp(a * I);
return
(t - 1/t) / (2 * I);
}
else
{
return
Float.sin(Float(a));
}
}
function
cos(a)
{
var
t;
if
(a
instanceof
Complex || a
instanceof
Series) {
t = exp(a * I);
return
(t + 1/t) / 2;
}
else
{
return
Float.cos(Float(a));
}
}
function
tan(a)
{
if
(a
instanceof
Complex || a
instanceof
Series) {
return
sin(a) / cos(a);
}
else
{
return
Float.tan(Float(a));
}
}
function
asin(a)
{
return
Float.asin(Float(a));
}
function
acos(a)
{
return
Float.acos(Float(a));
}
function
atan(a)
{
return
Float.atan(Float(a));
}
function
atan2(a, b)
{
return
Float.atan2(Float(a), Float(b));
}
function
sinc(a)
{
if
(a == 0) {
return
1;
}
else
{
a *= Float.PI;
return
sin(a) / a;
}
}
function
todeg(a)
{
return
a * 180 / Float.PI;
}
function
fromdeg(a)
{
return
a * Float.PI / 180;
}
function
sinh(a)
{
var
e = Float.exp(Float(a));
return
(e - 1/e) * 0.5;
}
function
cosh(a)
{
var
e = Float.exp(Float(a));
return
(e + 1/e) * 0.5;
}
function
tanh(a)
{
var
e = Float.exp(Float(a) * 2);
return
(e - 1) / (e + 1);
}
function
asinh(a)
{
var
x = Float(a);
return
log(sqrt(x * x + 1) + x);
}
function
acosh(a)
{
var
x = Float(a);
return
log(sqrt(x * x - 1) + x);
}
function
atanh(a)
{
var
x = Float(a);
return
0.5 * log((1 + x) / (1 - x));
}
function
sigmoid(x)
{
x = Float(x);
return
1 / (1 + exp(-x));
}
function
lerp(a, b, t)
{
return
a + (b - a) * t;
}
var
idn = Matrix.idn;
var
diag = Matrix.diag;
var
trans = Matrix.trans;
var
trace = Matrix.trace;
var
charpoly = Matrix.charpoly;
var
eigenvals = Matrix.eigenvals;
var
det = Matrix.det;
var
rank = Matrix.rank;
var
ker = Matrix.ker;
var
cp = Matrix.cp;
var
dp = Matrix.dp;
var
polroots = Polynomial.roots;
var
bestappr = Float.bestappr;