Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
calcurate Bernoulli number with rational computation
(function(n){
var ORD0 = '0'.charCodeAt(0);
var addu = function(u, v) {
var i, j, a, c, w;
w = ''
j = u.length;
i = v.length;
c = 0;
while (i > 0) {
a = u.charCodeAt(--j) - ORD0 + v.charCodeAt(--i) - ORD0 + c;
c = a >= 10 ? 1 : 0;
if (c == 1)
a -= 10;
w = String.fromCharCode(a + ORD0) + w;
}
while (j > 0) {
a = u.charCodeAt(--j) - ORD0 + c;
c = a >= 10 ? 1 : 0;
if (c == 1)
a -= 10;
w = String.fromCharCode(a + ORD0) + w;
}
if (c != 0)
w = String.fromCharCode(c + ORD0) + w;
return w;
};
var subu = function(u, v) {
var i, j, a, c, w;
w = ''
j = u.length;
i = v.length;
c = 1;
while (i > 0) {
a = u.charCodeAt(--j) - v.charCodeAt(--i) + c + 9;
c = a >= 10 ? 1 : 0;
if (c == 1)
a -= 10;
w = String.fromCharCode(a + ORD0) + w;
}
while (j > 0) {
a = u.charCodeAt(--j) - ORD0 + c + 9;
c = a >= 10 ? 1 : 0;
if (c == 1)
a -= 10;
w = String.fromCharCode(a + ORD0) + w;
}
return w;
};
var mulu = function(uu, vv) {
var m, i, j, a, c, x, u=[], w=[], ww;
m = uu.length;
j = vv.length;
for (i = 0; i < m; ++i)
u[i] = uu.charCodeAt(i) - ORD0;
for (i = 0; i < m + j; ++i)
w[i] = 0;
while (--j >= 0) {
x = vv.charCodeAt(j) - ORD0;
if (x > 0) {
i = m;
c = 0;
while (--i >= 0) {
a = u[i] * x + w[i + j + 1] + c;
c = Math.floor(a / 10);
w[i + j + 1] = a - c * 10;
}
w[j] = c;
}
}
ww = '';
for (i = 0, m = w.length; i < m; ++i)
ww = ww + String.fromCharCode(w[i] + ORD0);
return ww;
};
var divmodus = function(u, v) {
var q, r, i, n, a;
q = '';
r = 0;
for (i = 0, n = u.length; i < n; ++i) {
r = r * 10 + u.charCodeAt(i) - ORD0;
a = Math.floor(r / v);
r -= a * v;
q = q + String.fromCharCode(a + ORD0);
}
return [q, '' + r];
};
var divmodu = function(uu, vv) {
var d, i, j, c, c1, c2, m, n, qhat, u01, u2, a, a1, a2, c, u=[0], v=[0], q='';
m = uu.length - vv.length;
n = vv.length;
for (i = 0; i < n + m; ++i)
u[i + 1] = uu.charCodeAt(i) - ORD0;
for (i = 0; i < n; ++i)
v[i] = vv.charCodeAt(i) - ORD0;
d = Math.floor(10 / (vv.charCodeAt(0) - ORD0 + 1));
if (d > 1) {
i = m + n + 1;
c = 0;
while (--i >= 0) {
a = u[i] * d + c;
c = Math.floor(a / 10);
u[i] = a - c * 10;
}
i = n;
c = 0;
while (--i >= 0) {
a = v[i] * d + c;
c = Math.floor(a / 10);
v[i] = a - c * 10;
}
}
for (j = 0; j <= m; ++j) {
u01 = u[j] * 10 + u[j + 1];
u2 = j + 2 >= u.length ? 0 : u[j + 2];
qhat = Math.floor(u01 / v[0]);
if (qhat > 9)
qhat = 9;
while (qhat * v[1] > (u01 - qhat * v[0]) * 10 + u2)
--qhat;
if (qhat > 0) {
c1 = 0;
c2 = 1;
i = n;
while (--i >= 0) {
a1 = qhat * v[i] + c1;
c1 = Math.floor(a1 / 10);
a1 -= c1 * 10;
a2 = u[j + i + 1] - a1 + c2 + 9;
c2 = a2 >= 10 ? 1 : 0;
u[j + i + 1] = c2 == 1 ? (a2 - 10) : a2;
}
if (u[j] - c1 + c2 < 1) {
--qhat;
c = 0;
i = n;
while (--i >= 0) {
a = u[j + i + 1] + v[i] + c;
c = a >= 10 ? 1 : 0;
u[j + i + 1] = c > 0 ? (a - 10) : a;
}
}
}
u[j] = 0;
q = q + String.fromCharCode(qhat + ORD0);
}
uu = '';
if (d <= 1) {
for (i = m + 1; i <= m + n; ++i)
uu = uu + String.fromCharCode(u[i] + ORD0);
}
else {
c = 0;
for (i = 0; i < n; ++i) {
a = c * 10 + u[m + i + 1];
a1 = Math.floor(a / d);
c = a - a1 * d;
uu = uu + String.fromCharCode(a1 + ORD0);
}
}
return [q, uu];
};
var add = function(u, v) {
var su='', sv='', un, vn;
u = '' + u; v = '' + v;
if (u.charAt(0) == '-') {
su = '-';
u = u.slice(1);
}
if (v.charAt(0) == '-') {
sv = '-';
v = v.slice(1);
}
un = u.length;
vn = v.length;
if (((su == '-') && (sv == '')) || ((su == '') && (sv == '-'))) {
if (un > vn || ((un == vn) && (('' + u) >= v))) {
u = subu(u, v);
}
else {
su = sv;
u = subu(v, u);
}
}
else {
if (un > vn || ((un == vn) && (('' + u) >= v)))
u = addu(u, v);
else
u = addu(v, u);
}
u = u.replace(/^0+/, '');
if (u == '')
u = '0';
else if (su == '-')
u = su + u;
return u;
};
var bigabs = function(u) {
if (u.charAt(0) == '-')
u = u.slice(1);
return u;
};
var neg = function(u) {
if (u.charAt(0) == '-')
u = u.slice(1);
else if (u != '0')
u = '-' + u;
return u;
};
var sub = function(u, v) { return add(u, neg(v)); };
var mul = function(u, v) {
var su='', sv='';
u = '' + u; v = '' + v;
if ((u == '0') || (v == '0'))
return '0';
if (u == '1')
return v;
if (u == '-1')
return neg(v);
if (v == '1')
return u;
if (v == '-1')
return neg(u);
if (u.charAt(0) == '-') {
su = '-';
u = u.slice(1);
}
if (v.charAt(0) == '-') {
sv = '-';
v = v.slice(1);
}
u = mulu(u, v);
u = u.replace(/^0+/, '');
if (u == '')
u = '0';
else if ((su == '-') ^ (sv == '-'))
u = '-' + u;
return u;
};
var cmp = function(u, v) {
var su='', sv='', nu, nv;
u = '' + u; v = '' + v;
if (u.charAt(0) == '-') {
su = '-';
}
if (v.charAt(0) == '-') {
sv = '-';
}
if ((su == '-') && (sv == ''))
return -1;
if ((su == '') && (sv == '-'))
return +1;
if ((su == '') && (sv == '')) {
nu = u.length; nv = v.length;
return nu < nv ? -1 : nu > nv ? +1 : u < v ? -1 : u == v ? 0 : +1;
}
else {
u = u.slice(1); v = v.slice(1);
nu = u.length; nv = v.length;
return nu < nv ? +1 : nu > nv ? -1 : u < v ? +1 : u == v ? 0 : -1;
}
};
var divmod = function(u, v) {
var su='', sv='', sq = '', sr = '', qr, cp;
if (v == '0')
throw 'DIV/0';
if (u == '0')
return ['0', '0'];
if (v == '1')
return [u, '0'];
if (v == '-1')
return [neg(u), '0'];
u = '' + u; v = '' + v;
if (u.charAt(0) == '-') {
su = '-';
u = u.slice(1);
}
if (v.charAt(0) == '-') {
sv = '-';
v = v.slice(1);
}
if (v.length == 1) {
qr = divmodus(u, v);
}
else {
cp = cmp(u, v);
qr = cp < 0 ? ['0', u] : cp == 0 ? ['1', '0'] : divmodu(u, v);
}
if ((su == '-') && (sv == '-')) {
sr = '-';
}
else if ((su == '') && (sv == '-')) {
sq = '-'; sr = '-';
qr[0] = addu(qr[0], '1');
qr[1] = subu(v, qr[1]);
}
else if ((su == '-') && (sv == '')) {
sq = '-';
qr[0] = addu(qr[0], '1');
qr[1] = subu(v, qr[1]);
}
qr[0] = qr[0].replace(/^0+/, '');
if (qr[0] == '')
qr[0] = '0';
else if (sq == '-')
qr[0] = '-' + qr[0];
qr[1] = qr[1].replace(/^0+/, '');
if (qr[1] == '')
qr[1] = '0';
else if (sr == '-')
qr[1] = '-' + qr[1];
return qr;
};
var div = function(u, v) { return divmod(u, v)[0]; };
var gcd = function(u, v) {
var t;
u = bigabs(u); v = bigabs(v);
if (u == '0')
return v;
if (v == '0')
return u;
if (cmp(u, v) < 0) {
t = u; u = v; v = t;
}
while (v != '0') {
t = divmod(u, v)[1];
u = v;
v = t;
}
return u;
};
var fraction_reduction = function(n, d) {
var sn='', sd='';
if (n.charAt(0) == '-') {
sn = '-';
n = n.slice(1);
}
if (d.charAt(0) == '-') {
sd = '-';
d = d.slice(1);
}
var sign = (sn == '-') ^ (sd == '-') ? '-' : '';
var f = gcd(n, d);
return sign + div(n, f) + '/' + div(d, f);
};
var bernoulli_number = function(n) {
var bn='0', bd='1', f='1', i, j, m, r, sn, sd, p;
for (j = 1; j <= n; ++j) {
sn = '0'; sd = '1';
for (m = 1; m <= n - j + 1; ++m) {
r = '1';
for (i = m; i <= m + j - 1; ++i)
r = mul(r, i);
sn = add(mul(sn, m + j), mul(sd, r)); sd = mul(sd, m + j);
}
f = mul(f, -j);
p = '1';
for (i = 0; i < n; ++i)
p = mul(p, j);
sn = mul(p, sn); sd = mul(f, sd);
bn = add(mul(bn, sd), mul(bd, sn)); bd = mul(bd, sd);
}
return fraction_reduction(bn, bd);
};
return bernoulli_number(n);
})(22);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.