From ce446fc526d75d8215d585e279812ec553f32132 Mon Sep 17 00:00:00 2001 From: Hugues Malphettes Date: Mon, 4 Aug 2014 09:32:18 +0800 Subject: [PATCH] Avoid modifying the Array.prototype; add a test for each bessel function --- bessel.js | 50 +++++++++++++++++++++++++------------------------- test.js | 22 +++++++++++++++++++++- 2 files changed, 46 insertions(+), 26 deletions(-) diff --git a/bessel.js b/bessel.js index 88c0eef..c260b52 100644 --- a/bessel.js +++ b/bessel.js @@ -1,5 +1,5 @@ var M = Math; -Array.prototype.horner = function(v) { return this.reduce(function(z,w){return v * z + w;},0); }; +function _horner(arr, v) { return arr.reduce(function(z,w){return v * z + w;},0); }; function _bessel_iter(x, n, f0, f1, sign) { if(!sign) sign = -1; var tdx = 2 / x, f2; @@ -32,14 +32,14 @@ var besselj = (function() { function bessel0(x) { var a, a1, a2, y = x * x, xx = M.abs(x) - 0.785398164; if(M.abs(x) < 8) { - a1 = b0_a1a.horner(y); - a2 = b0_a2a.horner(y); + a1 = _horner(b0_a1a, y); + a2 = _horner(b0_a2a, y); a = a1/a2; } else { y = 64 / y; - a1 = b0_a1b.horner(y); - a2 = b0_a2b.horner(y); + a1 = _horner(b0_a1b, y); + a2 = _horner(b0_a2b, y); a = M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x)); } return a; @@ -51,13 +51,13 @@ var besselj = (function() { function bessel1(x) { var a, a1, a2, y = x*x, xx = M.abs(x) - 2.356194491; if(Math.abs(x)< 8) { - a1 = x*b1_a1a.horner(y); - a2 = b1_a2a.horner(y); + a1 = x*_horner(b1_a1a, y); + a2 = _horner(b1_a2a, y); a = a1 / a2; } else { y = 64 / y; - a1=b1_a1b.horner(y); - a2=b1_a2b.horner(y); + a1=_horner(b1_a1b, y); + a2=_horner(b1_a2b, y); a=M.sqrt(W/M.abs(x))*(M.cos(xx)*a1-M.sin(xx)*a2*8/M.abs(x)); if(x < 0) a = -a; } @@ -108,13 +108,13 @@ var bessely = (function() { function bessel0(x) { var a, a1, a2, y = x * x, xx = x - 0.785398164; if(x < 8) { - a1 = b0_a1a.horner(y); - a2 = b0_a2a.horner(y); + a1 = _horner(b0_a1a, y); + a2 = _horner(b0_a2a, y); a = a1/a2 + W * besselj(x,0) * M.log(x); } else { y = 64 / y; - a1 = b0_a1b.horner(y); - a2 = b0_a2b.horner(y); + a1 = _horner(b0_a1b, y); + a2 = _horner(b0_a2b, y); a = M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x); } return a; @@ -127,13 +127,13 @@ var bessely = (function() { function bessel1(x) { var a, a1, a2, y = x*x, xx = x - 2.356194491; if(x < 8) { - a1 = x*b1_a1a.horner(y); - a2 = b1_a2a.horner(y); + a1 = x*_horner(b1_a1a, y); + a2 = _horner(b1_a2a, y); a = a1/a2 + W * (besselj(x,1) * M.log(x) - 1 / x); } else { y = 64 / y; - a1=b1_a1b.horner(y); - a2=b1_a2b.horner(y); + a1=_horner(b1_a1b, y); + a2=_horner(b1_a2b, y); a=M.sqrt(W/x)*(M.sin(xx)*a1+M.cos(xx)*a2*8/x); } return a; @@ -145,15 +145,15 @@ var besseli = (function() { var b0_a = [1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732, 0.360768e-1, 0.45813e-2].reverse(); var b0_b = [0.39894228, 0.1328592e-1, 0.225319e-2, -0.157565e-2, 0.916281e-2, -0.2057706e-1, 0.2635537e-1, -0.1647633e-1, 0.392377e-2].reverse(); function bessel0(x) { - if(x <= 3.75) return b0_a.horner(x*x/(3.75*3.75)); - return M.exp(M.abs(x))/M.sqrt(M.abs(x))*b0_b.horner(3.75/M.abs(x)); + if(x <= 3.75) return _horner(b0_a, x*x/(3.75*3.75)); + return M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b0_b, 3.75/M.abs(x)); } var b1_a = [0.5, 0.87890594, 0.51498869, 0.15084934, 0.2658733e-1, 0.301532e-2, 0.32411e-3].reverse(); var b1_b = [0.39894228, -0.3988024e-1, -0.362018e-2, 0.163801e-2, -0.1031555e-1, 0.2282967e-1, -0.2895312e-1, 0.1787654e-1, -0.420059e-2].reverse(); function bessel1(x) { - if(x < 3.75) return x * b1_a.horner(x*x/(3.75*3.75)); - return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*b1_b.horner(3.75/M.abs(x)); + if(x < 3.75) return x * _horner(b1_a, x*x/(3.75*3.75)); + return (x < 0 ? -1 : 1) * M.exp(M.abs(x))/M.sqrt(M.abs(x))*_horner(b1_b, 3.75/M.abs(x)); } return function besseli(x, n) { @@ -187,15 +187,15 @@ var besselk = (function() { var b0_a = [-0.57721566, 0.42278420, 0.23069756, 0.3488590e-1, 0.262698e-2, 0.10750e-3, 0.74e-5].reverse(); var b0_b = [1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1, 0.587872e-2, -0.251540e-2, 0.53208e-3].reverse(); function bessel0(x) { - if(x <= 2) return -M.log(x/2)*besseli(x,0) + b0_a.horner(x*x/4); - return M.exp(-x)/M.sqrt(x)*b0_b.horner(2/x); + if(x <= 2) return -M.log(x/2)*besseli(x,0) + _horner(b0_a, x*x/4); + return M.exp(-x)/M.sqrt(x)*_horner(b0_b, 2/x); } var b1_a = [1.0, 0.15443144, -0.67278579, -0.18156897, -0.1919402e-1, -0.110404e-2, -0.4686e-4].reverse(); var b1_b = [1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1, -0.780353e-2, 0.325614e-2, -0.68245e-3].reverse(); function bessel1(x) { - if(x <= 2) return M.log(x/2)*besseli(x,1) + (1/x)*b1_a.horner(x*x/4); - return M.exp(-x)/M.sqrt(x)*b1_b.horner(2/x); + if(x <= 2) return M.log(x/2)*besseli(x,1) + (1/x)*_horner(b1_a, x*x/4); + return M.exp(-x)/M.sqrt(x)*_horner(b1_b, 2/x); } return _bessel_wrap(bessel0, bessel1, 'BESSELK', 2, 1); diff --git a/test.js b/test.js index 836aec3..bb45acb 100644 --- a/test.js +++ b/test.js @@ -1,2 +1,22 @@ +var assert = require('assert'); var bessel; -describe('source', function() { it('should load', function() { bessel = require('./'); }); }); +describe('When using bessel functions', function() { + before(function() { + bessel = require('./bessel'); + }); + it('It must load besseli', function() { + assert.equal(typeof bessel.besseli, 'function'); + }); + it('It must compute the Modified Bessel function at 1.5 with an order of 1 (0.981666)', function() { + assert(Math.floor((bessel.besseli(1.5, 1) - 0.981666) * 10e6) < 10); + }); + it('It must compute the Bessel function at 1.9 with an order of 2 (0.329926)', function() { + assert(Math.floor((bessel.besselj(1.9, 2) - 0.329926) * 10e6) < 10); + }); + it('It must compute the Modified Bessel function at 1.5 with an order of 1 (0.277388)', function() { + assert(Math.floor((bessel.besselk(1.5, 1) - 0.277388) * 10e6) < 10); + }); + it('It must compute Weber\'s Bessel function at 2.5 with an order of 1 (0.145918)', function() { + assert(Math.floor((bessel.besselk(2.5, 1) - 0.981666) * 10e6) < 10); + }); +}); -- 2.34.1