Code: Select all
function biased(n){return Math.random()<(1/n)};
function unbiased(n){var a;while((a=biased(n))==biased(n));return a};
Code: Select all
function KoggeStoneAdd (A, B)
{
// Addition without "+" operator
var G, P;
G = A & B;
P = A ^ B;
G |= P & (G << 1);
P = P & (P << 1);
G |= P & (G << 2);
P = P & (P << 2);
G |= P & (G << 4);
P = P & (P << 4);
G |= P & (G << 8);
P = P & (P << 8);
G |= P & (G << 16);
return (A ^ B) ^ (G << 1);
}
function KoggeStoneSub (A, B)
{
// Subtraction without "-" operator
var G, P, T;
T = ~A;
G = T & B;
P = T ^ B;
G |= P & (G << 1);
P = P & (P << 1);
G |= P & (G << 2);
P = P & (P << 2);
G |= P & (G << 4);
P = P & (P << 4);
G |= P & (G << 8);
P = P & (P << 8);
G |= P & (G << 16);
return (A ^ B) ^ (G << 1);
}
function KoggeStoneAdd3 (A, B, C)
{
// Addition without "+" operator
var T = A & B;
B ^= A;
A = B & C;
return KoggeStoneAdd((T | A) << 1, B ^ C);
}
function RevIntKoggeStoneAdd (A, B)
{
// Reversed integer addition
var G, P;
G = A & B;
P = A ^ B;
G |= P & (G >>> 1);
P = P & (P >>> 1);
G |= P & (G >>> 2);
P = P & (P >>> 2);
G |= P & (G >>> 4);
P = P & (P >>> 4);
G |= P & (G >>> 8);
P = P & (P >>> 8);
G |= P & (G >>> 16);
return (A ^ B) ^ (G >>> 1);
}
function RevIntKoggeStoneSub (A, B)
{
// Reversed integer subtraction
var G, P, T;
T = ~A;
G = T & B;
P = T ^ B;
G |= P & (G >>> 1);
P = P & (P >>> 1);
G |= P & (G >>> 2);
P = P & (P >>> 2);
G |= P & (G >>> 4);
P = P & (P >>> 4);
G |= P & (G >>> 8);
P = P & (P >>> 8);
G |= P & (G >>> 16);
return (A ^ B) ^ (G >>> 1);
}
function RevIntKoggeStoneAdd3 (A, B, C)
{
// Reversed integer addition
var T = A & B;
B ^= A;
A = B & C;
return RevIntKoggeStoneAdd((T | A) >>> 1, B ^ C);
}
function BigIntKoggeStoneAdd (A, B)
{
// Addition without "+" operator
// A and B must be same sign
var G, P, i, n1 = BigInt(1);
A = BigInt(A);
B = BigInt(B);
i = n1;
G = A & B;
P = A ^ B;
while (P)
{
G |= P & (G << i);
P = P & (P << i);
i <<= n1;
}
return (A ^ B) ^ (G << n1);
}
Code: Select all
function pcg32_constructor (state, inc)
{
// PCG-XSH-RR 64/32 (LCG)
let mult = BigInt("6364136223846793005");
let n1 = BigInt(1), n18 = BigInt(18), n27 = BigInt(27), n59 = BigInt(59), mask = BigInt("18446744073709551615");
inc = (BigInt(inc) & mask) | n1;
state = BigInt(state) & mask;
return function ()
{
let oldstate = state;
state = (oldstate * mult + inc) & mask;
let xorshifted = Number(((oldstate >> n18) ^ oldstate) >> n27) | 0;
let rot = Number(oldstate >> n59);
return (xorshifted >>> rot) | (xorshifted << (-rot & 31));
};
};
// Usage:
pcg32 = pcg32_constructor("5573589319906701683", "1442695040888963407");
pcg32 (); // 676697322
pcg32 (); // 420258633
pcg32 (); // -876335118
Code: Select all
// http://www.netlib.org/cephes/
var ellipj = (function() {
let sqrt = Math.sqrt;
let fabs = Math.abs;
let sin = Math.sin;
let cos = Math.cos;
let asin = Math.asin;
let tanh = Math.tanh;
let sinh = Math.sinh;
let cosh = Math.cosh;
let atan = Math.atan;
let exp = Math.exp;
let PIO2 = Math.PI / 2;
let MACHEP = Number.EPSILON / 2;
return function(u, m) {
let sn, cn, dn, ph;
let ai, b, phi, t, twon;
let a = []
, c = [];
let i;
/* Check for special cases */
if (m < 0.0 || m > 1.0) {
return ({
sn: NaN,
cn: NaN,
dn: NaN,
ph: NaN
});
}
if (m < 1.0e-9) {
t = sin(u);
b = cos(u);
ai = 0.25 * m * (u - t * b);
return ({
sn: t - ai * b,
cn: b + ai * t,
ph: u - ai,
dn: 1.0 - 0.5 * m * t * t
});
}
if (m >= 0.9999999999) {
ai = 0.25 * (1.0 - m);
b = cosh(u);
t = tanh(u);
phi = 1.0 / b;
twon = b * sinh(u);
sn = t + ai * (twon - u) / (b * b);
ph = 2.0 * atan(exp(u)) - PIO2 + ai * (twon - u) / b;
ai *= t * phi;
cn = phi - ai * (twon - u);
dn = phi + ai * (twon + u);
return ({
sn: sn,
cn: cn,
dn: dn,
ph: ph
});
}
/* A. G. M. scale */
a[0] = 1.0;
b = sqrt(1.0 - m);
c[0] = sqrt(m);
twon = 1.0;
i = 0;
while (fabs(c[i] / a[i]) > MACHEP) {
if (i > 7) {
console.warn("Overflow in ellipj");
break;
}
ai = a[i];
++i;
c[i] = (ai - b) / 2.0;
t = sqrt(ai * b);
a[i] = (ai + b) / 2.0;
b = t;
twon *= 2.0;
}
/* backward recurrence */
phi = twon * a[i] * u;
do {
t = c[i] * sin(phi) / a[i];
b = phi;
phi = (asin(t) + phi) / 2.0;
} while (--i);
t = sin(phi);
sn = t;
cn = cos(phi);
/* Thanks to Hartmut Henkel for reporting a bug here: */
dn = sqrt(1.0 - m * t * t);
ph = phi;
return ({
sn: sn,
cn: cn,
dn: dn,
ph: ph
});
}
})();
Code: Select all
var ellik = (function(phi, m) {
if (m > 1) {
let k = Math.sqrt(m);
return ellik(Math.asin(k * Math.sin(phi)), 1 / m) / k
}
let a, b, r, t, a1, b1, eps, i;
phi /= 2;
t = Math.cos(phi) / Math.sin(phi);
a = t * t + (1 - 2 * m);
b = Math.sign(a + 1 / a) * Math.sqrt(a * a + 4 * m * (1 - m));
r = Math.sqrt(1 + 2 * m / (a + b));
t = t / r;
a1 = 1;
b1 = Math.sqrt(1 - m);
eps = Math.max(a1, b1) * Number.EPSILON;
i = 0;
do {
a = a1;
b = b1;
a1 = (a + b) / 2;
b1 = Math.sqrt(a * b);
r = Math.sqrt(a1 * (r + 1) / (b * r + a));
t *= r;
} while (i++ < 100 && (Math.abs(a - b) > eps || Math.abs(r - 1) > Number.EPSILON));
t = Math.atan(a / t);
b = phi - t;
r = Math.sign(b) * Math.PI;
t += r * Math.round(b / r);
return 2 * t / a;
});
Code: Select all
var nearest = (function() {
"use strict";
function norm(w) {
return w[0] * w[0] + w[1] * w[1];
}
function proj(w1, w2) {
return (w1[0] * w2[0] + w1[1] * w2[1]) / norm(w2);
}
function nint(x) {
return Math.sign(x) * Math.ceil(Math.abs(x) - 0.5);
}
function reduce(w1, w2) {
var i, t;
w1 = [+w1[0], +w1[1]];
w2 = [+w2[0], +w2[1]];
if (norm(w1) > norm(w2)) t = w1, w1 = w2, w2 = t;
for (i = 0; i < 1000; i++) {
t = nint(proj(w2, w1));
if (!t) break;
w2[0] -= t * w1[0];
w2[1] -= t * w1[1];
t = w1, w1 = w2, w2 = t;
}
if (norm(w1) > norm(w2)) t = w1, w1 = w2, w2 = t;
return [w1, w2];
}
function nearesti(p, w) {
var n = nint(proj(p, w));
p[0] -= n * w[0];
p[1] -= n * w[1];
return [n, norm(p)];
}
return function (p, w1, w2) {
var w = reduce(w1, w2);
var a = proj(w[1], w[0]);
var b = proj(p, [w[1][0] - a * w[0][0], w[1][1] - a * w[0][1]]);
var c = nint(b);
var d = c + Math.sign(b - c);
var e = nearesti([p[0] - c * w[1][0], p[1] - c * w[1][1]], w[0]);
if (c !== d) {
a = nearesti([p[0] - d * w[1][0], p[1] - d * w[1][1]], w[0]);
if (e[1] > a[1]) c = d, e = a;
}
d = e[0];
return [d * w[0][0] + c * w[0][1], d * w[1][0] + c * w[1][1]];
}
})();Code: Select all
function twosum(a, b) {
let s = a + b;
let v = s - a;
let e = (a - (s - v)) + (b - v);
return [s, e]
}
function split(a) {
let t = 134217729 * a;
let hi = t - (t - a);
let lo = a - hi;
return [hi, lo];
}
function twoprod(a, b) {
let p = a * b;
let [ahi, alo] = split(a);
let [bhi, blo] = split(b);
let e = ((ahi * bhi - p) + ahi * blo + alo * bhi) + alo * blo;
return [p, e];
}
function fma(a, b, c) {
[a, b] = twoprod(a, b);
[b, c] = twosum(b, c);
return a + b;
}Code: Select all
function GiniIndex (list) {
let i, n, sum, eps, max;
list = Array.from(list).sort((x, y) => y - x);
max = list[0];
for (i = sum = eps = 0, n = list.length; i < n; i++) {
let y = (max - list[i]) - eps;
let t = sum + y;
eps = (t - sum) - y;
sum = t;
}
return sum / ((n - 1) * max);
}Code: Select all
function integrator1 (a, b, h, f, steps) {
a = +a; b = +b;
for (let i = 0; i < steps; i++) {
let c = f(a, b);
a += h*(b+c*h/2);
b += c*h;
let d = f(a, b);
b += (d-c)*h/2;
}
return [a, b];
}
function integrator2 (a, b, h, f, steps) {
a = +a; b = +b;
for (let i = 0; i < steps; i++) {
let c = f(a, b);
a += h*(b+c*h/2);
b += c*h;
let d = f(a, b);
b += (d-c)*h/2;
let e = f(a, b) - d;
if (e !== 0) b += e/(1-e/(d-c))*h/2;
}
return [a, b];
}
