正規分布の近似(14)

結局、分数クラスを作った。
思っていたほどは苦労しなかったが、
やはりデバッグがいやな作業だった。

p10,0 = 1/362880x^9
p10,1 = 1/36288-1/4032x+1/1008x^2-1/432x^3+1/288x^4-1/288x^5+1/432x^6-1/1008x^7+1/4032x^8-1/40320x^9
p10,2 = -329/5184+1151/4032x-575/1008x^2+287/432x^3-143/288x^4+71/288x^5-35/432x^6+17/1008x^7-1/504x^8+1/10080x^9
p10,3 = 233893/36288-11083/576x+3667/144x^2-8461/432x^3+2773/288x^4-901/288x^5+289/432x^6-13/144x^7+1/144x^8-1/4320x^9
p10,4 = -5271131/36288+185525/576x-45485/144x^2+77555/432x^3-18731/288x^4+4475/288x^5-1055/432x^6+35/144x^7-1/72x^8+1/2880x^9
p10,5 = 43947619/36288-1220725/576x+235765/144x^2-316195/432x^3+60019/288x^4-11275/288x^5+2095/432x^6-55/144x^7+5/288x^8-1/2880x^9
p10,6 = -167683997/36288+3818123/576x-604043/144x^2+663581/432x^3-103277/288x^4+15941/288x^5-2441/432x^6+53/144x^7-1/72x^8+1/4320x^9
p10,7 = 316559287/36288-6064393/576x+807745/144x^2-748207/432x^3+98407/288x^4-12871/288x^5+1675/432x^6-31/144x^7+1/144x^8-1/10080x^9
p10,8 = -287420489/36288+33046721/4032x-3782969/1008x^2+431441/432x^3-49049/288x^4+5561/288x^5-629/432x^6+71/1008x^7-1/504x^8+1/40320x^9
p10,9 = 1562500/567-156250/63x+62500/63x^2-6250/27x^3+625/18x^4-125/36x^5+25/108x^6-5/504x^7+1/4032x^8-1/362880x^9

定数項がかなり違う。
まあ、循環にもなっていない小数を分数にするのは無理がある。
20秒もかかったが、速くする方法はないのか?


ソースは隠しておく。



function poly() {
var args = poly.arguments;
if(args.length == 0)
this.a = [ new fraction() ];
else {
this.a = [ ];
for(var i = 0; i < args.length; i++) {
var b = args[i];
if(typeof b == "number")
this.a.push(new fraction(b));
else if(b instanceof fraction)
this.a.push(b.copy());
else if(b instanceof Array) {
for(var j = 0; j < b.length; j++)
this.a.push(b[j]);
}
}
}

// コピー
this.copy = function() {
return new poly(this.a);
}

// 加算
this.add = function(p) {
var r = this.copy();
if(typeof p == "number" || p instanceof fraction)
r.a[0] = r.a[0].add(p);
else if(p instanceof poly) {
for(var i = 0; i < p.a.length; i++) {
if(i < r.a.length)
r.a[i] = r.a[i].add(p.a[i]);
else
r.a[i] = p.a[i].copy();
}
}
else
throw("unsupported argument in poly::add");

r.normalize();
return r;
}

// 減算
this.subtract = function(p) {
var r = this.copy();
if(typeof p == "number" || p instanceof fraction)
r.a[0] = r.a[0].subtract(p);
else if(p instanceof poly) {
for(var i = 0; i < p.a.length; i++) {
if(i < r.a.length)
r.a[i] = r.a[i].subtract(p.a[i]);
else {
r.a[i] = p.a[i].copy();
r.a[i].a = -r.a[i].a;
}
}
}
else
throw("unsupported argument in poly::subtract");

r.normalize();
return r;
}

// 乗算
this.multiply = function(p) {
var r = new poly();
if(typeof p == "number" || p instanceof fraction) {
for(var i = 0; i < this.a.length; i++)
r.a[i] = r.a[i].multiply(p);
}
else if(p instanceof poly) {
for(var i = 0; i < this.a.length; i++) {
for(var j = 0; j < p.a.length; j++) {
if(r.a[i+j] == undefined)
r.a[i+j] = this.a[i].multiply(p.a[j]);
else
r.a[i+j] = r.a[i+j].add(this.a[i].multiply(p.a[j]));
}
}
}
else
throw("unsupported argument in poly::multiply");

r.normalize();
return r;
}

// 関数の値
this.value = function(x) {
if(!(typeof x == "number" || x instanceof fraction))
throw("unsupported argument in poly::value");

var r = new fraction();
for(var i = this.a.length - 1; i >= 0; i--)
r = r.multiply(x).add(this.a[i]);

return r;
}

// 合成
this.synthesize = function(p) {
if(!(p instanceof poly))
throw("unsupported argument in poly::synthesize");

var r = new poly();
for(var i = this.a.length - 1; i >= 0; i--)
r = r.multiply(p).add(this.a[i]);

return r;
}

// 積分
this.integral = function() {
var args = this.integral.arguments;
if(args.length == 0) { // 不定積分
var r = new poly();
for(var i = 0; i < this.a.length; i++)
r.a[i+1] = this.a[i].divide(i + 1);
return r;
}
else if(args.length == 2) { // 定積分
var r1 = [ ];
var q = this.integral(); // とりあえず不定積分する

for(var i = 0; i < 2; i++) {
var a = args[i];
if(a instanceof poly) // 多項式なら合成
r1[i] = q.synthesize(a);
else
r1[i] = q.value(a);
}

// 種類によって引き算を場合分け
if(r1[1] instanceof poly) {
return r1[1].subtract(r1[0]);
}
else {
if(r1[0] instanceof poly)
return (new poly(r1[1]).subtract(r1[0]));
else
return r1[1].subtract(r1[0]);
}
}
else
throw("wrong number of arguments in poly::integral");
}

// 上位の係数0の項は削除する(内部関数)
this.normalize = function() {
// 真の配列の長さを求める
var len = this.a.length;
for( ; len > 1; len--) {
if(this.a[len-1].a != 0)
break;
}
this.a.length = len;
}
this.toString = function() {
if(this.a.length == 1) {
return "" + this.a[0];
}
var str = "";
var first = true;
for(var i = 0; i < this.a.length; i++) {
var c = this.a[i];
if(c.a == 0)
continue;
else if(c.a > 0) {
if(!first)
str += "+";
if(c != 1 || i == 0)
str += c;
}
else {
if(c.a == -c.b && i > 0)
str += "-";
else
str += c;
}
if(i == 1)
str += "x";
else if(i > 1)
str += "x^" + i;

first = false;
}
return str;
}
}

function fraction() {
this.a; // 分子
this.b; // 分母

// 引数の個数によって変える
var args = fraction.arguments;
if(args.length == 0) {
this.a = 0;
this.b = 1;
}
else if(args.length == 1) {
this.a = args[0];
this.b = 1;
}
else {
this.a = args[0]; // 分子
this.b = args[1]; // 分母
}

this.copy = function() {
return new fraction(this.a, this.b);
}
this.add = function(f) {
if(typeof f == "number")
return new fraction(this.a + this.b * f, this.b);
else if(f instanceof fraction)
return new fraction(this.a * f.b + this.b * f.a, this.b * f.b);
else
throw("wrong argument in fraction::add");
}
this.subtract = function(f) {
if(typeof f == "number")
return new fraction(this.a - this.b * f, this.b);
else if(f instanceof fraction)
return new fraction(this.a * f.b - this.b * f.a, this.b * f.b);
else
throw("wrong argument in fraction::add");
}
this.multiply = function(f) {
if(typeof f == "number")
return new fraction(this.a * f, this.b);
else if(f instanceof fraction)
return new fraction(this.a * f.a, this.b * f.b);
else
throw("wrong argument in fraction::add");
}
this.divide = function(f) {
if(typeof f == "number")
return new fraction(this.a, this.b * f);
else if(f instanceof fraction)
return new fraction(this.a * f.b, this.b * f.a);
else
throw("wrong argument in fraction::add");
}

this.normalize = function() {
// ユークリッドの互除法
var a = this.a;
var b = this.b;
var a2, f;
while(true) {
if((a2 = a % b) == 0) {
f = b;
break;
}
if((b %= a) == 0) {
f = a;
break;
}
a = a2;
}
this.a /= f;
this.b /= f;
if(this.b < 0) {
this.a = -this.a;
this.b = -this.b;
}
}
this.toString = function() {
if(this.b == 1)
return this.a;
else
return this.a + "/" + this.b;
}

this.normalize();
}