From 73ee33e4023ed62347caf249e7faae7c7b932609 Mon Sep 17 00:00:00 2001 From: SheetJS Date: Sat, 24 Sep 2016 02:15:52 -0400 Subject: [PATCH] version bump 0.3.0: cleanup - eliminated array reduce (performance improvement) - explicit treatment of divergent values - better handling of non-integral order (fixes #3, h/t @vron) - reorganized source, removed voc dependency (fixes #2, h/t @hmalphettes) - new browser demo - more detailed test suite and coverage tests - updated travis versions for test - miscellaneous adjustments to tooling --- .flowconfig | 19 ++++ .gitignore | 1 + .jscs.json | 6 ++ .jshintrc | 4 + .travis.yml | 12 ++- Makefile | 69 +++++++++++- README.md | 94 ++++++++++++---- bessel.flow.js | 252 +++++++++++++++++++++++++++++++++++++++++++ bessel.js | 141 +++++++++++++++--------- bessel.md | 10 +- bits/00_header.js | 25 +++++ bits/01_version.js | 1 + bits/02_helpers.js | 2 + bits/10_types.js | 4 + bits/20_helpers.js | 1 + bits/30_iter.js | 10 ++ bits/40_wrapper.js | 14 +++ bits/50_besselj.js | 83 ++++++++++++++ bits/60_bessely.js | 45 ++++++++ bits/70_besseli.js | 43 ++++++++ bits/80_besselk.js | 19 ++++ bits/90_exports.js | 4 + bits/99_footer.js | 1 + index.html | 84 +++++++++++++++ misc/help.sh | 42 ++++++++ misc/spin.sh | 14 +++ package.json | 29 +++-- test.html | 4 - test.js | 142 ++++++++++++++++++++---- test_files/excel.tsv | 60 +++++++++++ test_files/mma.tsv | 60 +++++++++++ 31 files changed, 1177 insertions(+), 118 deletions(-) create mode 100644 .flowconfig create mode 100644 .jscs.json create mode 100644 .jshintrc create mode 100644 bessel.flow.js create mode 100644 bits/00_header.js create mode 100644 bits/01_version.js create mode 100644 bits/02_helpers.js create mode 100644 bits/10_types.js create mode 100644 bits/20_helpers.js create mode 100644 bits/30_iter.js create mode 100644 bits/40_wrapper.js create mode 100644 bits/50_besselj.js create mode 100644 bits/60_bessely.js create mode 100644 bits/70_besseli.js create mode 100644 bits/80_besselk.js create mode 100644 bits/90_exports.js create mode 100644 bits/99_footer.js create mode 100644 index.html create mode 100755 misc/help.sh create mode 100755 misc/spin.sh delete mode 100644 test.html create mode 100644 test_files/excel.tsv create mode 100644 test_files/mma.tsv diff --git a/.flowconfig b/.flowconfig new file mode 100644 index 0000000..39d82f6 --- /dev/null +++ b/.flowconfig @@ -0,0 +1,19 @@ +[ignore] +.*/node_modules/.* +.*/dist/.* +.*/test.js + +.*/bits/.* +.*/ctest/.* +.*/misc/.* +.*/perf/.* + +.*/bessel.js + +[include] +bessel.flow.js + +[libs] +bits/10_types.js + +[options] diff --git a/.gitignore b/.gitignore index 3c3629e..865fd19 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ node_modules +misc/coverage.html diff --git a/.jscs.json b/.jscs.json new file mode 100644 index 0000000..5965004 --- /dev/null +++ b/.jscs.json @@ -0,0 +1,6 @@ +{ + "requireCommaBeforeLineBreak": true, + "disallowTrailingWhitespace": true, + "disallowTrailingComma": true +} + diff --git a/.jshintrc b/.jshintrc new file mode 100644 index 0000000..5cd205b --- /dev/null +++ b/.jshintrc @@ -0,0 +1,4 @@ +{ + "bitwise": false, + "curly": false +} diff --git a/.travis.yml b/.travis.yml index dba96e7..8beb852 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,7 +1,15 @@ language: node_js node_js: - - "0.11" + - "6" + - "5" + - "4.2" + - "0.12" - "0.10" - "0.8" before_install: - - "npm install -g mocha" + - "npm install -g npm@next" + - "npm install -g mocha voc" + - "npm install blanket" + - "npm install coveralls mocha-lcov-reporter" +after_success: + - "make coveralls-spin" diff --git a/Makefile b/Makefile index bc8c57c..756e598 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,66 @@ -LIBRARY=bessel +LIB=bessel +REQS= +ADDONS= +AUXTARGETS= +HTMLLINT=index.html -$(LIBRARY).js: bessel.md - node_modules/.bin/voc $^ > $@ +ULIB=$(shell echo $(LIB) | tr a-z A-Z) +DEPS=$(sort $(wildcard bits/*.js)) +TARGET=$(LIB).js +FLOWTARGET=$(LIB).flow.js -test mocha: - mocha -R spec +## Main Targets + +.PHONY: all +all: $(TARGET) $(AUXTARGETS) ## Build library and auxiliary scripts + +$(TARGET) $(AUXTARGETS): %.js : %.flow.js + node -e 'process.stdout.write(require("fs").readFileSync("$<","utf8").replace(/^[ \t]*\/\*[:#][^*]*\*\/\s*(\n)?/gm,"").replace(/\/\*[:#][^*]*\*\//gm,""))' > $@ + +$(FLOWTARGET): $(DEPS) + cat $^ | tr -d '\15\32' > $@ + +bits/01_version.js: package.json + echo "BESSEL.version = '"`grep version package.json | awk '{gsub(/[^0-9a-z\.-]/,"",$$2); print $$2}'`"';" > $@ + +.PHONY: clean +clean: ## Remove targets and build artifacts + rm -f $(TARGET) $(FLOWTARGET) + +## Testing + +.PHONY: test mocha +test mocha: test.js $(TARGET) ## Run test suite + mocha -R spec -t 20000 + +## Code Checking + +.PHONY: lint +lint: $(TARGET) $(AUXTARGETS) ## Run jshint and jscs checks + @jshint --show-non-errors $(TARGET) $(AUXTARGETS) + @jshint --show-non-errors package.json + @jshint --show-non-errors --extract=always $(HTMLLINT) + @jscs $(TARGET) $(AUXTARGETS) + +.PHONY: flow +flow: lint ## Run flow checker + @flow check --all --show-all-errors + +.PHONY: cov +cov: misc/coverage.html ## Run coverage test + +misc/coverage.html: $(TARGET) test.js + mocha --require blanket -R html-cov -t 20000 > $@ + +.PHONY: coveralls +coveralls: ## Coverage Test + Send to coveralls.io + mocha --require blanket --reporter mocha-lcov-reporter -t 20000 | node ./node_modules/coveralls/bin/coveralls.js + + +.PHONY: help +help: + @grep -hE '(^[a-zA-Z_-][ a-zA-Z_-]*:.*?|^#[#*])' $(MAKEFILE_LIST) | bash misc/help.sh + +#* To show a spinner, append "-spin" to any target e.g. cov-spin +%-spin: + @make $* & bash misc/spin.sh $$! diff --git a/README.md b/README.md index 3ce6439..3e8133d 100644 --- a/README.md +++ b/README.md @@ -1,36 +1,92 @@ -# Bessel.JS +# bessel -Pure-JS implementation of the Bessel functions (J,Y,I,K), for node and browser +Pure-JS implementation of Bessel functions J,Y,I,K (for the browser and nodejs). +Emphasis on correctness and performance for integer order. The standard notation is used here: - - J is the Bessel function of the first kind - - Y is the Bessel function of the second kind - - I is the modified Bessel function of the first kind - - K is the modified Bessel function of the first kind + - `J` is the Bessel function of the first kind + - `Y` is the Bessel function of the second kind + - `I` is the modified Bessel function of the first kind + - `K` is the modified Bessel function of the second kind -# Usage +## Installation -The functions `besselj`, `bessely`, `besseli`, `besselk` are exposed when you include -the script `bessel.js`: +With [npm](https://www.npmjs.org/package/bessel): -```html> +```bash +$ npm install bessel +``` + +In the browser: + +```html - ``` -See `test.html` for an example +The browser exposes a variable `BESSEL` -In node, those four functions are exported: +The script will manipulate `module.exports` if available (e.g. in a CommonJS +`require` context). This is not always desirable. To prevent the behavior, +define `DO_NOT_EXPORT_BESSEL` -```js> -var besselj01 = require('bessel').besselj(0,1); +## Usage + +In all cases, the relevant function takes two arguments (`value`, `order`). + +The return value is a JS number. `NaN` signals an error in calculation. + +- `BESSEL.besselj(x, n)` computes `J_{n}(x)` + +- `BESSEL.bessely(x, n)` computes `Y_{n}(x)` + +- `BESSEL.besseli(x, n)` computes `I_{n}(x)` + +- `BESSEL.besselk(x, n)` computes `K_{n}(x)` + +For example: + +```js +> // var BESSEL = require('bessel'); // uncomment this line if in node +> BESSEL.besselj(1.5,0) // 0.5118276712499389 +> BESSEL.bessely(1.5,0) // 0.38244892476502895 +> BESSEL.besseli(1.5,0) // 1.6467232021476754 +> BESSEL.besselk(1.5,0) // 0.2138055693236539 + +> var Y = BESSEL.bessely +> Y(Math.PI, 5) + Y(Math.PI, 3) - (2 * 4 / Math.PI) * Y(Math.PI, 4) // 0 ``` -Each function follows Excel semantics `(value, function-order)`. For example, +## Testing -```js> -bessel.besselj(1.5, 1) +`make test` will run the nodejs-based test. + +To generate the `excel.tsv` test cases, make a 6-column Excel Sheet: + +| x | n | `BESSELI` | `BESSELJ` | `BESSELK` | `BESSELY` | +|---|---|:--------------:|:--------------:|:--------------:|:--------------:| +| x | n |`BESSELI(A1,B1)`|`BESSELJ(A1,B1)`|`BESSELK(A1,B1)`|`BESSELY(A1,B1)`| + +To generate the `mma.tsv` test cases, use the Mathematica Bessel Functions: + +```mathematica +(* Bessel_[value, order] *) +F[x_,n_]:={x/2,n,BesselI[n,x/2], BesselJ[n,x/2], BesselK[n,x/2], BesselY[n,x/2]} +Do[Print[ExportString[N[F[x,n],10],"csv"]], {n,1,3}, {x,1,20} ] ``` -is the value of the bessel function J1 at the point x=1.5 +Note: Each function follows Excel semantics `(value, order)`. Other platforms +like Mathematica reverse the argument order. + +## License + +Please consult the attached LICENSE file for details. All rights not explicitly +granted by the MIT License are reserved by the Original Author. + +## Badges + +[![Build Status](https://travis-ci.org/SheetJS/bessel.svg?branch=master)](https://travis-ci.org/SheetJS/bessel) + +[![Coverage Status](http://img.shields.io/coveralls/SheetJS/bessel/master.svg)](https://coveralls.io/r/SheetJS/bessel?branch=master) + +[![Analytics](https://ga-beacon.appspot.com/UA-36810333-1/SheetJS/bessel?pixel)](https://github.com/SheetJS/bessel) diff --git a/bessel.flow.js b/bessel.flow.js new file mode 100644 index 0000000..2336c91 --- /dev/null +++ b/bessel.flow.js @@ -0,0 +1,252 @@ +/* bessel.js (C) 2013-present SheetJS -- http://sheetjs.com */ +/* vim: set ts=2: */ +/*exported BESSEL */ +var BESSEL; +/*:: declare var DO_NOT_EXPORT_BESSEL: any; */ +/*:: declare var define: any; */ +(function (factory) { + /*jshint ignore:start */ + if(typeof DO_NOT_EXPORT_BESSEL === 'undefined') { + if('object' === typeof exports) { + factory(exports); + } else if ('function' === typeof define && define.amd) { + define(function () { + var module = {}; + factory(module); + return module; + }); + } else { + factory(BESSEL = {}); + } + } else { + factory(BESSEL = {}); + } + /*jshint ignore:end */ +}(function(BESSEL) { +BESSEL.version = '0.3.0'; +var M = Math; + +/*:: +type BesselN = (x:number) => number; +type BesselF = (x:number, n:number) => number; +*/ +function _horner(arr/*:Array*/, v/*:number*/)/*:number*/ { for(var i = 0, z = 0; i < arr.length; ++i) z = v * z + arr[i]; return z; } +function _bessel_iter(x/*:number*/, n/*:number*/, f0/*:number*/, f1/*:number*/, sign/*:number*/)/*:number*/ { + if(n === 0) return f0; + if(n === 1) return f1; + var tdx = 2 / x, f2 = f1; + for(var o = 1; o < n; ++o) { + f2 = f1 * o * tdx + sign * f0; + f0 = f1; f1 = f2; + } + return f2; +} +function _bessel_wrap(bessel0/*:BesselN*/, bessel1/*:BesselN*/, name/*:string*/, nonzero/*:number*/, sign/*:number*/)/*:BesselF*/ { + return function bessel(x/*:number*/,n/*:number*/) { + if(nonzero) { + if(x === 0) return (nonzero == 1 ? -Infinity : Infinity); + else if(x < 0) return NaN; + } + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(n < 0) return NaN; + n|=0; + var b0/*:number*/ = bessel0(x), b1/*:number*/ = bessel1(x); + return _bessel_iter(x, n, b0, b1, sign); + }; +} +var besselj/*:BesselF*/ = (function() { + var W = 0.636619772; // 2 / Math.PI + + var b0_a1a = [57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456].reverse(); + var b0_a2a = [57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse(); + + function bessel0(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x * x; + if(x < 8) { + a1 = _horner(b0_a1a, y); + a2 = _horner(b0_a2a, y); + a = a1 / a2; + } else { + var xx = x - 0.785398164; + y = 64 / y; + a1 = _horner(b0_a1b, y); + a2 = _horner(b0_a2b, y); + a = M.sqrt(W/x)*(M.cos(xx)*a1-M.sin(xx)*a2*8/x); + } + return a; + } + + var b1_a1a = [72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606].reverse(); + var b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + + function bessel1(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x*x, xx = M.abs(x) - 2.356194491; + if(Math.abs(x)< 8) { + a1 = x*_horner(b1_a1a, y); + a2 = _horner(b1_a2a, y); + a = a1 / a2; + } else { + y = 64 / 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; + } + return a; + } + + return function besselj(x/*:number*/, n/*:number*/)/*:number*/ { + n = Math.round(n); + if(!isFinite(x)) return isNaN(x) ? x : 0; + if(n < 0) return ((n%2)?-1:1)*besselj(x, -n); + if(x < 0) return ((n%2)?-1:1)*besselj(-x, n); + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(x === 0) return 0; + + var ret=0.0; + if(x > n) { + ret = _bessel_iter(x, n, bessel0(x), bessel1(x),-1); + } else { + var m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2); + var jsum=false; + var bjp=0.0, sum=0.0; + var bj=1.0, bjm = 0.0; + var tox = 2 / x; + for (var j=m;j>0;j--) { + bjm=j*tox*bj-bjp; + bjp=bj; + bj=bjm; + if (M.abs(bj) > 1E10) { + bj *= 1E-10; + bjp *= 1E-10; + ret *= 1E-10; + sum *= 1E-10; + } + if (jsum) sum += bj; + jsum=!jsum; + if (j == n) ret=bjp; + } + sum=2.0*sum-bj; + ret /= sum; + } + return ret; + }; +})(); +var bessely/*:BesselF*/ = (function() { + var W = 0.636619772; + + var b0_a1a = [-2957821389.0, 7062834065.0, -512359803.6, 10879881.29, -86327.92757, 228.4622733].reverse(); + var b0_a2a = [40076544269.0, 745249964.8, 7189466.438, 47447.26470, 226.1030244, 1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934945152e-7].reverse(); + + function bessel0(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x * x, xx = x - 0.785398164; + if(x < 8) { + 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 = _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; + } + + var b1_a1a = [-0.4900604943e13, 0.1275274390e13, -0.5153438139e11, 0.7349264551e9, -0.4237922726e7, 0.8511937935e4].reverse(); + var b1_a2a = [0.2499580570e14, 0.4244419664e12, 0.3733650367e10, 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + + function bessel1(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x*x, xx = x - 2.356194491; + if(x < 8) { + 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=_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; + } + + return _bessel_wrap(bessel0, bessel1, 'BESSELY', 1, -1); +})(); +var besseli/*:BesselF*/ = (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/*:number*/)/*:number*/ { + 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/*:number*/)/*:number*/ { + 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/*:number*/, n/*:number*/)/*:number*/ { + n = Math.round(n); + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(n < 0) return NaN; + if(M.abs(x) === 0) return 0; + if(x == Infinity) return Infinity; + + var ret = 0.0, j, tox = 2 / M.abs(x), bip = 0.0, bi=1.0, bim=0.0; + var m=2*M.round((n+M.round(M.sqrt(40*n)))/2); + for (j=m;j>0;j--) { + bim=j*tox*bi + bip; + bip=bi; bi=bim; + if (M.abs(bi) > 1E10) { + bi *= 1E-10; + bip *= 1E-10; + ret *= 1E-10; + } + if(j == n) ret = bip; + } + ret *= besseli(x, 0) / bi; + return x < 0 && (n%2) ? -ret : ret; + }; + +})(); + +var besselk/*:BesselF*/ = (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/*:number*/)/*:number*/ { + 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/*:number*/)/*:number*/ { + 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); +})(); +BESSEL.besselj = besselj; +BESSEL.bessely = bessely; +BESSEL.besseli = besseli; +BESSEL.besselk = besselk; +})); diff --git a/bessel.js b/bessel.js index c260b52..c69f73e 100644 --- a/bessel.js +++ b/bessel.js @@ -1,55 +1,85 @@ +/* bessel.js (C) 2013-present SheetJS -- http://sheetjs.com */ +/* vim: set ts=2: */ +/*exported BESSEL */ +var BESSEL; +(function (factory) { + /*jshint ignore:start */ + if(typeof DO_NOT_EXPORT_BESSEL === 'undefined') { + if('object' === typeof exports) { + factory(exports); + } else if ('function' === typeof define && define.amd) { + define(function () { + var module = {}; + factory(module); + return module; + }); + } else { + factory(BESSEL = {}); + } + } else { + factory(BESSEL = {}); + } + /*jshint ignore:end */ +}(function(BESSEL) { +BESSEL.version = '0.3.0'; var M = Math; -function _horner(arr, v) { return arr.reduce(function(z,w){return v * z + w;},0); }; + +function _horner(arr, v) { for(var i = 0, z = 0; i < arr.length; ++i) z = v * z + arr[i]; return z; } function _bessel_iter(x, n, f0, f1, sign) { - if(!sign) sign = -1; - var tdx = 2 / x, f2; if(n === 0) return f0; if(n === 1) return f1; - for(var o = 1; o != n; ++o) { + var tdx = 2 / x, f2 = f1; + for(var o = 1; o < n; ++o) { f2 = f1 * o * tdx + sign * f0; f0 = f1; f1 = f2; } - return f1; + return f2; } function _bessel_wrap(bessel0, bessel1, name, nonzero, sign) { return function bessel(x,n) { + if(nonzero) { + if(x === 0) return (nonzero == 1 ? -Infinity : Infinity); + else if(x < 0) return NaN; + } if(n === 0) return bessel0(x); if(n === 1) return bessel1(x); - if(n < 0) throw name + ': Order (' + n + ') must be nonnegative'; - if(nonzero == 1 && x === 0) throw name + ': Undefined when x == 0'; - if(nonzero == 2 && x <= 0) throw name + ': Undefined when x <= 0'; + if(n < 0) return NaN; + n|=0; var b0 = bessel0(x), b1 = bessel1(x); return _bessel_iter(x, n, b0, b1, sign); }; } var besselj = (function() { - var b0_a1a = [57568490574.0,-13362590354.0,651619640.7,-11214424.18,77392.33017,-184.9052456].reverse(); - var b0_a2a = [57568490411.0,1029532985.0,9494680.718,59272.64853,267.8532712,1.0].reverse(); - var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); - var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse(); var W = 0.636619772; // 2 / Math.PI + var b0_a1a = [57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456].reverse(); + var b0_a2a = [57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse(); + function bessel0(x) { - var a, a1, a2, y = x * x, xx = M.abs(x) - 0.785398164; - if(M.abs(x) < 8) { + var a=0, a1=0, a2=0, y = x * x; + if(x < 8) { a1 = _horner(b0_a1a, y); a2 = _horner(b0_a2a, y); - a = a1/a2; - } - else { + a = a1 / a2; + } else { + var xx = x - 0.785398164; y = 64 / 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)); + a = M.sqrt(W/x)*(M.cos(xx)*a1-M.sin(xx)*a2*8/x); } return a; } - var b1_a1a = [72362614232.0,-7895059235.0,242396853.1,-2972611.439, 15704.48260, -30.16036606].reverse(); + + var b1_a1a = [72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606].reverse(); var b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse(); var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + function bessel1(x) { - var a, a1, a2, y = x*x, xx = M.abs(x) - 2.356194491; + var a=0, a1=0, a2=0, y = x*x, xx = M.abs(x) - 2.356194491; if(Math.abs(x)< 8) { a1 = x*_horner(b1_a1a, y); a2 = _horner(b1_a2a, y); @@ -63,22 +93,26 @@ var besselj = (function() { } return a; } + return function besselj(x, n) { n = Math.round(n); - if(n === 0) return bessel0(M.abs(x)); - if(n === 1) return bessel1(M.abs(x)); - if(n < 0) throw 'BESSELJ: Order (' + n + ') must be nonnegative'; - if(M.abs(x) === 0) return 0; + if(!isFinite(x)) return isNaN(x) ? x : 0; + if(n < 0) return ((n%2)?-1:1)*besselj(x, -n); + if(x < 0) return ((n%2)?-1:1)*besselj(-x, n); + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(x === 0) return 0; - var ret, j, tox = 2 / M.abs(x), m, jsum, sum, bjp, bj, bjm; - if(M.abs(x) > n) { - ret = _bessel_iter(x, n, bessel0(M.abs(x)), bessel1(M.abs(x)),-1); + var ret=0.0; + if(x > n) { + ret = _bessel_iter(x, n, bessel0(x), bessel1(x),-1); } else { - m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2); - jsum=0; - bjp=ret=sum=0.0; - bj=1.0; - for (j=m;j>0;j--) { + var m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2); + var jsum=false; + var bjp=0.0, sum=0.0; + var bj=1.0, bjm = 0.0; + var tox = 2 / x; + for (var j=m;j>0;j--) { bjm=j*tox*bj-bjp; bjp=bj; bj=bjm; @@ -95,18 +129,19 @@ var besselj = (function() { sum=2.0*sum-bj; ret /= sum; } - return x < 0 && (n%2) ? -ret : ret; + return ret; }; })(); var bessely = (function() { + var W = 0.636619772; + var b0_a1a = [-2957821389.0, 7062834065.0, -512359803.6, 10879881.29, -86327.92757, 228.4622733].reverse(); var b0_a2a = [40076544269.0, 745249964.8, 7189466.438, 47447.26470, 226.1030244, 1.0].reverse(); var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934945152e-7].reverse(); - var W = 0.636619772; function bessel0(x) { - var a, a1, a2, y = x * x, xx = x - 0.785398164; + var a=0, a1=0, a2=0, y = x * x, xx = x - 0.785398164; if(x < 8) { a1 = _horner(b0_a1a, y); a2 = _horner(b0_a2a, y); @@ -124,8 +159,9 @@ var bessely = (function() { var b1_a2a = [0.2499580570e14, 0.4244419664e12, 0.3733650367e10, 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1].reverse(); var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + function bessel1(x) { - var a, a1, a2, y = x*x, xx = x - 2.356194491; + var a=0, a1=0, a2=0, y = x*x, xx = x - 2.356194491; if(x < 8) { a1 = x*_horner(b1_a1a, y); a2 = _horner(b1_a2a, y); @@ -144,6 +180,7 @@ var bessely = (function() { 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 _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)); @@ -151,6 +188,7 @@ var besseli = (function() { 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 * _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)); @@ -159,14 +197,13 @@ var besseli = (function() { return function besseli(x, n) { n = Math.round(n); if(n === 0) return bessel0(x); - if(n == 1) return bessel1(x); - if(n < 0) throw 'BESSELI Order (' + n + ') must be nonnegative'; + if(n === 1) return bessel1(x); + if(n < 0) return NaN; if(M.abs(x) === 0) return 0; + if(x == Infinity) return Infinity; - var ret, j, tox = 2 / M.abs(x), m, bip, bi, bim; - m=2*M.round((n+M.round(M.sqrt(40*n)))/2); - bip=ret=0.0; - bi=1.0; + var ret = 0.0, j, tox = 2 / M.abs(x), bip = 0.0, bi=1.0, bim=0.0; + var m=2*M.round((n+M.round(M.sqrt(40*n)))/2); for (j=m;j>0;j--) { bim=j*tox*bi + bip; bip=bi; bi=bim; @@ -186,24 +223,24 @@ var besseli = (function() { 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) + _horner(b0_a, x*x/4); - return M.exp(-x)/M.sqrt(x)*_horner(b0_b, 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)*_horner(b1_a, x*x/4); + 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); })(); -if(typeof exports !== "undefined") { - exports.besselj = besselj; - exports.bessely = bessely; - exports.besseli = besseli; - exports.besselk = besselk; -} - +BESSEL.besselj = besselj; +BESSEL.bessely = bessely; +BESSEL.besseli = besseli; +BESSEL.besselk = besselk; +})); diff --git a/bessel.md b/bessel.md index b532abc..c31f947 100644 --- a/bessel.md +++ b/bessel.md @@ -14,7 +14,7 @@ var M = Math; The methods use an approximating polynomial and evaluate using Horner's method: ``` -function _horner(arr, v) { return arr.reduce(function(z,w){return v * z + w;},0); }; +function _horner(arr/*:Array*/, v/*:number*/) { return arr.reduce(function(z,w){return v * z + w;},0); } ``` ## Recurrence @@ -29,7 +29,7 @@ So rather than go back and try to find solution for each order, we will build solutions for `n=0` and `n=1` and then apply the recurrence. The helper: ```js -function _bessel_iter(x, n, f0, f1, sign) { +function _bessel_iter(x/*:number*/, n/*:number*/, f0/*:number*/, f1/*:number*/, sign/*:?number*/) { if(!sign) sign = -1; var tdx = 2 / x, f2; if(n === 0) return f0; @@ -48,7 +48,7 @@ checks since `Y_n` is undefined at 0 and `K_n` is real only when `x>0` ``` function _bessel_wrap(bessel0, bessel1, name, nonzero, sign) { - return function bessel(x,n) { + return function bessel(x/*:number*/,n/*:number*/) { if(n === 0) return bessel0(x); if(n === 1) return bessel1(x); if(n < 0) throw name + ': Order (' + n + ') must be nonnegative'; @@ -143,7 +143,7 @@ For large values of x, the aforementioned iteration is fine, but for small values the expressions quickly blow up. Hence a more careful iteration is used: ``` - return function besselj(x, n) { + return function besselj(x/*:number*/, n/*:number*/) { n = Math.round(n); if(n === 0) return bessel0(M.abs(x)); if(n === 1) return bessel1(M.abs(x)); @@ -246,7 +246,7 @@ var besseli = (function() { 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) { + return function besseli(x/*:number*/, n/*:number*/) { n = Math.round(n); if(n === 0) return bessel0(x); if(n == 1) return bessel1(x); diff --git a/bits/00_header.js b/bits/00_header.js new file mode 100644 index 0000000..1b6f408 --- /dev/null +++ b/bits/00_header.js @@ -0,0 +1,25 @@ +/* bessel.js (C) 2013-present SheetJS -- http://sheetjs.com */ +/* vim: set ts=2: */ +/*exported BESSEL */ +var BESSEL; +/*:: declare var DO_NOT_EXPORT_BESSEL: any; */ +/*:: declare var define: any; */ +(function (factory) { + /*jshint ignore:start */ + if(typeof DO_NOT_EXPORT_BESSEL === 'undefined') { + if('object' === typeof exports) { + factory(exports); + } else if ('function' === typeof define && define.amd) { + define(function () { + var module = {}; + factory(module); + return module; + }); + } else { + factory(BESSEL = {}); + } + } else { + factory(BESSEL = {}); + } + /*jshint ignore:end */ +}(function(BESSEL) { diff --git a/bits/01_version.js b/bits/01_version.js new file mode 100644 index 0000000..025a37a --- /dev/null +++ b/bits/01_version.js @@ -0,0 +1 @@ +BESSEL.version = '0.3.0'; diff --git a/bits/02_helpers.js b/bits/02_helpers.js new file mode 100644 index 0000000..726d059 --- /dev/null +++ b/bits/02_helpers.js @@ -0,0 +1,2 @@ +var M = Math; + diff --git a/bits/10_types.js b/bits/10_types.js new file mode 100644 index 0000000..3b845db --- /dev/null +++ b/bits/10_types.js @@ -0,0 +1,4 @@ +/*:: +type BesselN = (x:number) => number; +type BesselF = (x:number, n:number) => number; +*/ diff --git a/bits/20_helpers.js b/bits/20_helpers.js new file mode 100644 index 0000000..441adec --- /dev/null +++ b/bits/20_helpers.js @@ -0,0 +1 @@ +function _horner(arr/*:Array*/, v/*:number*/)/*:number*/ { for(var i = 0, z = 0; i < arr.length; ++i) z = v * z + arr[i]; return z; } diff --git a/bits/30_iter.js b/bits/30_iter.js new file mode 100644 index 0000000..308feba --- /dev/null +++ b/bits/30_iter.js @@ -0,0 +1,10 @@ +function _bessel_iter(x/*:number*/, n/*:number*/, f0/*:number*/, f1/*:number*/, sign/*:number*/)/*:number*/ { + if(n === 0) return f0; + if(n === 1) return f1; + var tdx = 2 / x, f2 = f1; + for(var o = 1; o < n; ++o) { + f2 = f1 * o * tdx + sign * f0; + f0 = f1; f1 = f2; + } + return f2; +} diff --git a/bits/40_wrapper.js b/bits/40_wrapper.js new file mode 100644 index 0000000..c12ab9d --- /dev/null +++ b/bits/40_wrapper.js @@ -0,0 +1,14 @@ +function _bessel_wrap(bessel0/*:BesselN*/, bessel1/*:BesselN*/, name/*:string*/, nonzero/*:number*/, sign/*:number*/)/*:BesselF*/ { + return function bessel(x/*:number*/,n/*:number*/) { + if(nonzero) { + if(x === 0) return (nonzero == 1 ? -Infinity : Infinity); + else if(x < 0) return NaN; + } + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(n < 0) return NaN; + n|=0; + var b0/*:number*/ = bessel0(x), b1/*:number*/ = bessel1(x); + return _bessel_iter(x, n, b0, b1, sign); + }; +} diff --git a/bits/50_besselj.js b/bits/50_besselj.js new file mode 100644 index 0000000..0d18dc3 --- /dev/null +++ b/bits/50_besselj.js @@ -0,0 +1,83 @@ +var besselj/*:BesselF*/ = (function() { + var W = 0.636619772; // 2 / Math.PI + + var b0_a1a = [57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456].reverse(); + var b0_a2a = [57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse(); + + function bessel0(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x * x; + if(x < 8) { + a1 = _horner(b0_a1a, y); + a2 = _horner(b0_a2a, y); + a = a1 / a2; + } else { + var xx = x - 0.785398164; + y = 64 / y; + a1 = _horner(b0_a1b, y); + a2 = _horner(b0_a2b, y); + a = M.sqrt(W/x)*(M.cos(xx)*a1-M.sin(xx)*a2*8/x); + } + return a; + } + + var b1_a1a = [72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606].reverse(); + var b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + + function bessel1(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x*x, xx = M.abs(x) - 2.356194491; + if(Math.abs(x)< 8) { + a1 = x*_horner(b1_a1a, y); + a2 = _horner(b1_a2a, y); + a = a1 / a2; + } else { + y = 64 / 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; + } + return a; + } + + return function besselj(x/*:number*/, n/*:number*/)/*:number*/ { + n = Math.round(n); + if(!isFinite(x)) return isNaN(x) ? x : 0; + if(n < 0) return ((n%2)?-1:1)*besselj(x, -n); + if(x < 0) return ((n%2)?-1:1)*besselj(-x, n); + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(x === 0) return 0; + + var ret=0.0; + if(x > n) { + ret = _bessel_iter(x, n, bessel0(x), bessel1(x),-1); + } else { + var m=2*M.floor((n+M.floor(M.sqrt(40*n)))/2); + var jsum=false; + var bjp=0.0, sum=0.0; + var bj=1.0, bjm = 0.0; + var tox = 2 / x; + for (var j=m;j>0;j--) { + bjm=j*tox*bj-bjp; + bjp=bj; + bj=bjm; + if (M.abs(bj) > 1E10) { + bj *= 1E-10; + bjp *= 1E-10; + ret *= 1E-10; + sum *= 1E-10; + } + if (jsum) sum += bj; + jsum=!jsum; + if (j == n) ret=bjp; + } + sum=2.0*sum-bj; + ret /= sum; + } + return ret; + }; +})(); diff --git a/bits/60_bessely.js b/bits/60_bessely.js new file mode 100644 index 0000000..c8a6013 --- /dev/null +++ b/bits/60_bessely.js @@ -0,0 +1,45 @@ +var bessely/*:BesselF*/ = (function() { + var W = 0.636619772; + + var b0_a1a = [-2957821389.0, 7062834065.0, -512359803.6, 10879881.29, -86327.92757, 228.4622733].reverse(); + var b0_a2a = [40076544269.0, 745249964.8, 7189466.438, 47447.26470, 226.1030244, 1.0].reverse(); + var b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse(); + var b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934945152e-7].reverse(); + + function bessel0(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x * x, xx = x - 0.785398164; + if(x < 8) { + 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 = _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; + } + + var b1_a1a = [-0.4900604943e13, 0.1275274390e13, -0.5153438139e11, 0.7349264551e9, -0.4237922726e7, 0.8511937935e4].reverse(); + var b1_a2a = [0.2499580570e14, 0.4244419664e12, 0.3733650367e10, 0.2245904002e8, 0.1020426050e6, 0.3549632885e3, 1].reverse(); + var b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse(); + var b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse(); + + function bessel1(x/*:number*/)/*:number*/ { + var a=0, a1=0, a2=0, y = x*x, xx = x - 2.356194491; + if(x < 8) { + 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=_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; + } + + return _bessel_wrap(bessel0, bessel1, 'BESSELY', 1, -1); +})(); diff --git a/bits/70_besseli.js b/bits/70_besseli.js new file mode 100644 index 0000000..c67f1a9 --- /dev/null +++ b/bits/70_besseli.js @@ -0,0 +1,43 @@ +var besseli/*:BesselF*/ = (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/*:number*/)/*:number*/ { + 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/*:number*/)/*:number*/ { + 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/*:number*/, n/*:number*/)/*:number*/ { + n = Math.round(n); + if(n === 0) return bessel0(x); + if(n === 1) return bessel1(x); + if(n < 0) return NaN; + if(M.abs(x) === 0) return 0; + if(x == Infinity) return Infinity; + + var ret = 0.0, j, tox = 2 / M.abs(x), bip = 0.0, bi=1.0, bim=0.0; + var m=2*M.round((n+M.round(M.sqrt(40*n)))/2); + for (j=m;j>0;j--) { + bim=j*tox*bi + bip; + bip=bi; bi=bim; + if (M.abs(bi) > 1E10) { + bi *= 1E-10; + bip *= 1E-10; + ret *= 1E-10; + } + if(j == n) ret = bip; + } + ret *= besseli(x, 0) / bi; + return x < 0 && (n%2) ? -ret : ret; + }; + +})(); + diff --git a/bits/80_besselk.js b/bits/80_besselk.js new file mode 100644 index 0000000..ec35b61 --- /dev/null +++ b/bits/80_besselk.js @@ -0,0 +1,19 @@ +var besselk/*:BesselF*/ = (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/*:number*/)/*:number*/ { + 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/*:number*/)/*:number*/ { + 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/bits/90_exports.js b/bits/90_exports.js new file mode 100644 index 0000000..4940334 --- /dev/null +++ b/bits/90_exports.js @@ -0,0 +1,4 @@ +BESSEL.besselj = besselj; +BESSEL.bessely = bessely; +BESSEL.besseli = besseli; +BESSEL.besselk = besselk; diff --git a/bits/99_footer.js b/bits/99_footer.js new file mode 100644 index 0000000..be4600a --- /dev/null +++ b/bits/99_footer.js @@ -0,0 +1 @@ +})); diff --git a/index.html b/index.html new file mode 100644 index 0000000..ce125bf --- /dev/null +++ b/index.html @@ -0,0 +1,84 @@ + + + + + + Bessel Live Demo + + + + + Bessel Live Demo
+ Source Code
+ Issues? Something look weird? Click here and report an issue
+
+ + + + + + + + +
Number:
Order:
 
J
Y
I
K
+ + + + + diff --git a/misc/help.sh b/misc/help.sh new file mode 100755 index 0000000..7a1c4c4 --- /dev/null +++ b/misc/help.sh @@ -0,0 +1,42 @@ +#!/bin/bash +# make_help.sh -- process listing of targets and special items in Makefile +# Copyright (C) 2016-present SheetJS +# +# usage in makefile: pipe the output of the following command: +# @grep -hE '(^[a-zA-Z_-][ a-zA-Z_-]*:.*?|^#[#*])' $(MAKEFILE_LIST) +# +# lines starting with "## " are treated as subtitles +# lines starting with "#* " are treated as plaintext comments +# multiple targets with "## " after the ":" are rendered as separate targets +# if the presumed default target is labeled, it will be assigned a unique color + +awk ' +BEGIN{recipes=0;} + !/#[#*] .*$/ {next;} + {multi=0; isrecipe=0;} + /^[^#]*:/ {isrecipe=1; ++recipes;} + /^[^ :]* .*:/ {multi=1} + multi==0 && isrecipe>0 { if(recipes > 1) print; else print $0, "[default]"; next} + isrecipe == 0 {print; next} + multi>0 { + k=split($0, msg, "##"); m=split($0, a, ":"); n=split(a[1], b, " "); + for(i=1; i<=n; ++i) print b[i] ":", "##" msg[2], (recipes==1 && i==1 ? "[default]" : "") + } +END {} +' | if [[ -t 1 ]]; then +awk ' +BEGIN {FS = ":.*?## "} + {color=36} + /\[default\]/ {color=35} + NF==1 && /^##/ {color=34} + NF==1 && /^#\*/ {color=20; $1 = substr($1, 4)} + {printf "\033[" color "m%-20s\033[0m %s\n", $1, $2;} +END{}' - +else +awk ' +BEGIN {FS = ":.*?## "} + /^#\* / {$1 = substr($1, 4)} + {printf "%-20s %s\n", $1, $2;} +END{}' - +fi + diff --git a/misc/spin.sh b/misc/spin.sh new file mode 100755 index 0000000..471dfee --- /dev/null +++ b/misc/spin.sh @@ -0,0 +1,14 @@ +#!/bin/bash +# spin.sh -- show a spinner (for coverage test) +# Copyright (C) 2014-present SheetJS + +wpid=$1 +delay=1 +str="|/-\\" +while [ $(ps -a|awk '$1=='$wpid' {print $1}') ]; do + t=${str#?} + printf " [%c]" "$str" + str=$t${str%"$t"} + sleep $delay + printf "\b\b\b\b" +done diff --git a/package.json b/package.json index ec15861..6176285 100644 --- a/package.json +++ b/package.json @@ -1,19 +1,26 @@ { "name": "bessel", - "version": "0.2.0", + "version": "0.3.0", "author": "SheetJS", - "description": "Bessel Functions in pure JS", - "keywords": ["bessel", "math", "specfun"], - "main": "bessel.js", - "repository": { "type": "git", "url": "https://github.com/SheetJS/bessel" }, - "license": "MIT", + "description": "Pure-JS Bessel Functions", + "keywords": [ "bessel", "math", "specfun" ], + "main": "./bessel", + "dependencies": { + }, + "devDependencies": { + "mocha":"", + "voc":"" + }, + "repository": { "type":"git", "url":"git://github.com/SheetJS/bessel.git" }, "scripts": { "test": "make test" }, - "dependencies": { - "voc":"" + "config": { + "blanket": { + "pattern": "bessel.js" + } }, - "devDependencies": { - "mocha":"" - } + "bugs": { "url": "https://github.com/SheetJS/bessel/issues" }, + "license": "MIT", + "engines": { "node": ">=0.8" } } diff --git a/test.html b/test.html deleted file mode 100644 index 53c67ef..0000000 --- a/test.html +++ /dev/null @@ -1,4 +0,0 @@ - - diff --git a/test.js b/test.js index bb45acb..7dd4150 100644 --- a/test.js +++ b/test.js @@ -1,22 +1,124 @@ -var assert = require('assert'); -var bessel; -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); +/* vim: set ts=2: */ +var X; +if(typeof require !== 'undefined') { + assert = require('assert'); + describe('source',function(){ + it('should load',function(){X=require('./');}); + it('all bessel functions', function() { + assert.equal(typeof X.besseli, 'function'); + assert.equal(typeof X.besselj, 'function'); + assert.equal(typeof X.besselk, 'function'); + assert.equal(typeof X.bessely, 'function'); + }); }); +} else { X = BESSEL; } + +function approx(a, b, t) { return Math.abs((a - b) / (Math.abs(a) + Math.abs(b))) <= t; } + +function chk(func, x, n, v, F, t) { + var V = func(x,n); + if(!approx(V, v, t) && v != V && !(isNaN(v) && isNaN(V))) throw new Error(F + "(" + x + "," + n + ") = " + V + " !~ " + v); +} + +var funcs = ['besseli', 'besselj', 'besselk', 'bessely']; + +var tests = { +'besseli': [ + [1.5, 1, 0.981666], + [1.9, 2, 0.603272435], + [2.5, 1, 2.516716242] +], +'besselj': [ + [1.5, 1, 0.557936508], + [1.9, 2, 0.329925829], + [2.5, 1, 0.497094103] +], +'besselk': [ + [1.5, 1, 0.277388], + [1.9, 2, 0.296909301], + [2.5, 1, 0.073890816] +], +'bessely': [ + [1.5, 1, -0.412308627], + [1.9, 2, -0.669878674], + [2.5, 1, 0.14591814] +] }; + +describe('correctness', function() { + funcs.forEach(function(F) { it(F, function() { + tests[F].forEach(function(t) { chk(X[F], t[0], t[1], t[2], F, 1e-6); }); + });}); + + if(typeof require != 'undefined') { + var fs = require('fs'); + it('excel', function() { + var xl = fs.readFileSync("test_files/excel.tsv", 'ascii').split("\n").map(function(l) { return l.split("\t");}); + xl.forEach(function(l) { if(l.length < 6) return; + var x=Number(l[0]), n=Number(l[1]); + var i=Number(l[2]), j=Number(l[3]), k=Number(l[4]), y=Number(l[5]); + chk(X.besseli, x, n, i, 'besseli', 1e-4); + chk(X.besselj, x, n, j, 'besselj', 1e-4); + chk(X.besselk, x, n, k, 'besselk', 1e-4); + chk(X.bessely, x, n, y, 'bessely', 1e-4); + }); + }); + it('mma', function() { + var mma = fs.readFileSync("test_files/mma.tsv", 'ascii').split("\n").map(function(l) { return l.split("\t");}); + mma.forEach(function(l) { if(l.length < 6) return; + var x=Number(l[0]), n=Number(l[1]); + var i=Number(l[2]), j=Number(l[3]), k=Number(l[4]), y=Number(l[5]); + chk(X.besseli, x, n, i, 'besseli', 1e-5); + chk(X.besselj, x, n, j, 'besselj', 1e-5); + chk(X.besselk, x, n, k, 'besselk', 1e-5); + chk(X.bessely, x, n, y, 'bessely', 1e-5); + }); + }); + } +}); + +var rand = [ + 0.3946086812775882, 3.946086812775882, 39.46086812775882, + 0.03119419917289612, 0.3119419917289612, 3.119419917289612, + 0.18311903629136062, 1.8311903629136062, 18.311903629136062, + 0.5317681126599452, 5.317681126599452, 53.17681126599452, + 0.37724541457032235, 3.7724541457032235, 37.724541457032235, + 0.8287802926845655, 8.287802926845655, 82.87802926845655, + 0.7197724866379658, 7.197724866379658, 71.97724866379658 +]; + +describe('properties', function() { + var t = 1e-4; + it('limiting behavior', function() { + chk(X.besselj, 0, 0, 1, 'besselj', 1e-6); + chk(X.besseli, 0, 0, 1, 'besseli', 1e-6); + chk(X.besselk, 0, 0, Infinity, 'besselk', 1e-6); + chk(X.bessely, 0, 0, -Infinity, 'bessely', 1e-6); + for(var i = 1; i < 20; ++i) { + chk(X.besselj, 0, i, 0, 'besselj', 1e-6); + chk(X.besseli, 0, i, 0, 'besseli', 1e-6); + } + }); + it('besselj', function() { + var F = "BESSELJ", f = X.besselj; + rand.forEach(function(r) { + for(var n = 0; n < 20; ++n) { + var pp = f( r, n), pn = f( r, -n), np = f(-r, n), nn = f(-r, -n); + if((n%2)) { + /* besselj odd if n is odd */ + if(!approx(np, -pp, t)) throw new Error(F + " np[" + np + "] != -pp[" + (-pp) + "] (" + r + "," + n + ")"); + if(!approx(nn, -pn, t)) throw new Error(F + " nn[" + nn + "] != -pn[" + (-pn) + "] (" + r + "," + n + ")"); + /* asymmetric in n */ + if(!approx(pn, -pp, t)) throw new Error(F + " pn[" + pn + "] != -pp[" + (-pp) + "] (" + r + "," + n + ")"); + if(!approx(nn, -np, t)) throw new Error(F + " nn[" + nn + "] != -pn[" + (-pn) + "] (" + r + "," + n + ")"); + } else { + /* besselj even if n is even */ + if(!approx(np, pp, t)) throw new Error(F + " np[" + np + "] != pp[" + (pp) + "] (" + r + "," + n + ")"); + if(!approx(nn, pn, t)) throw new Error(F + " nn[" + nn + "] != pn[" + (pn) + "] (" + r + "," + n + ")"); + /* symmetric in n */ + if(!approx(pn, pp, t)) throw new Error(F + " pn[" + pn + "] != pp[" + (pp) + "] (" + r + "," + n + ")"); + if(!approx(nn, np, t)) throw new Error(F + " nn[" + nn + "] != pn[" + (pn) + "] (" + r + "," + n + ")"); + } + } + }); + }); }); diff --git a/test_files/excel.tsv b/test_files/excel.tsv new file mode 100644 index 0000000..0f03b09 --- /dev/null +++ b/test_files/excel.tsv @@ -0,0 +1,60 @@ +0.5 1 0.257894303 0.242268458 1.656441128 -1.471472392 +1 1 0.565159098 0.440050586 0.601907232 -0.781212821 +1.5 1 0.981666428 0.557936508 0.277387804 -0.412308627 +2 1 1.590636857 0.576724808 0.13986588 -0.107032432 +2.5 1 2.516716242 0.497094103 0.073890816 0.145918138 +3 1 3.953370217 0.339058958 0.040156431 0.324674424 +3.5 1 6.205834932 0.137377527 0.022239393 0.410188417 +4 1 9.759465157 -0.066043328 0.012483499 0.397925709 +4.5 1 15.38922293 -0.231060432 0.007078095 0.300997321 +5 1 24.33564185 -0.327579139 0.004044613 0.14786314 +5.5 1 38.58816378 -0.341438216 0.002325569 -0.023758243 +6 1 61.34193694 -0.276683859 0.00134392 -0.175010349 +6.5 1 97.73501304 -0.153841304 0.000779894 -0.274091281 +7 1 156.0390965 -0.004682826 0.000454182 -0.302667244 +7.5 1 249.5843679 0.135248424 0.000265297 -0.25912852 +8 1 399.8731348 0.234628387 0.000155369 -0.158060462 +8.5 1 641.6198938 0.273115291 9.11973E-05 -0.02616868 +9 1 1030.914709 0.245307412 5.3637E-05 0.104314575 +9.5 1 1658.453068 0.161262274 3.1602E-05 0.203179899 +10 1 2670.988321 0.04347225 1.86488E-05 0.249015424 +0.5 2 0.031906148 0.030604023 7.550183547 -5.441370834 +1 2 0.135747667 0.114903485 1.624838884 -1.650682613 +1.5 2 0.337834621 0.232087679 0.583655974 -0.932193761 +2 2 0.688948449 0.352834208 0.25375976 -0.617408098 +2.5 2 1.276466159 0.446059058 0.121460207 -0.381335848 +3 2 2.245212431 0.486091263 0.061510459 -0.160400399 +3.5 2 3.832012072 0.458629185 0.032307122 0.045371436 +4 2 6.422189499 0.364128143 0.017401426 0.215903599 +4.5 2 10.64151731 0.217848984 0.009545677 0.328481596 +5 2 17.50561501 0.046565119 0.005308944 0.367662879 +5.5 2 28.66258571 -0.117315483 0.00298437 0.330841237 +6 2 46.7870949 -0.242873213 0.001691968 0.229857908 +6.5 2 76.22054669 -0.307430389 0.000965899 0.088906659 +7 2 124.01131 -0.301417224 0.000554562 -0.0605266 +7.5 2 201.6054825 -0.230273411 0.000319924 -0.186414216 +8 2 327.5958389 -0.112993711 0.000185313 -0.263036605 +8.5 2 532.1925545 0.022323169 0.000107716 -0.276362442 +9 2 864.4962206 0.144846369 6.28006E-05 -0.226755682 +9.5 2 1404.333008 0.2278787 3.67109E-05 -0.128435911 +10 2 2281.518997 0.254630214 2.15098E-05 -0.005868083 +0.5 3 0.002645112 0.00256373 62.0579095 -42.05949428 +1 3 0.022168424 0.019563354 7.101262769 -5.821517632 +1.5 3 0.080774114 0.060963951 1.833803735 -2.073541401 +2 3 0.21273996 0.12894325 0.6473854 -1.127783765 +2.5 3 0.474370413 0.216600403 0.268227146 -0.756055495 +3 3 0.959753625 0.309062865 0.122170376 -0.538541623 +3.5 3 1.826392593 0.386770113 0.059161818 -0.358335346 +4 3 3.337275843 0.430171471 0.029884924 -0.18202211 +4.5 3 5.930096274 0.424703973 0.015563141 -0.009013679 +5 3 10.3311502 0.364831234 0.008291768 0.146267163 +5.5 3 17.7426483 0.256117865 0.00449602 0.264370052 +6 3 30.15054042 0.114768384 0.002471898 0.328248955 +6.5 3 50.83005827 -0.035346628 0.001374294 0.328803071 +7 3 85.17548647 -0.167555588 0.000771075 0.268080616 +7.5 3 142.0614436 -0.25806091 0.000435923 0.159707605 +8 3 236.0752263 -0.291125242 0.000248026 0.026542159 +8.5 3 391.1763672 -0.26261027 0.000141887 -0.103884234 +9 3 646.6942118 -0.180931248 8.15484E-05 -0.205094878 +9.5 3 1067.155009 -0.065313348 4.70593E-05 -0.257258177 +10 3 1758.380739 0.058379836 2.72527E-05 -0.251362657 diff --git a/test_files/mma.tsv b/test_files/mma.tsv new file mode 100644 index 0000000..5bd4c8c --- /dev/null +++ b/test_files/mma.tsv @@ -0,0 +1,60 @@ +0.5 1 0.2578943054 0.2422684577 1.65644112 -1.471472393 +1 1 0.565159104 0.4400505857 0.6019072302 -0.7812128213 +1.5 1 0.9816664286 0.5579365079 0.2773878005 -0.412308627 +2 1 1.590636855 0.5767248078 0.1398658818 -0.1070324315 +2.5 1 2.516716245 0.4970941025 0.07389081635 0.145918138 +3 1 3.953370217 0.3390589585 0.04015643113 0.3246744248 +3.5 1 6.205834922 0.1373775274 0.02223939293 0.4101884179 +4 1 9.759465154 -0.06604332802 0.01248349889 0.3979257106 +4.5 1 15.38922275 -0.2310604319 0.007078094909 0.3009973231 +5 1 24.33564214 -0.3275791376 0.004044613445 0.1478631434 +5.5 1 38.58816462 -0.3414382154 0.002325569009 -0.02375823896 +6 1 61.34193678 -0.2766838581 0.001343919718 -0.1750103443 +6.5 1 97.73501077 -0.1538413014 0.0007798943982 -0.274091274 +7 1 156.0390929 -0.004682823482 0.0004541824869 -0.302667237 +7.5 1 249.5843654 0.1352484276 0.0002652973901 -0.2591285105 +8 1 399.8731368 0.2346363469 0.0001553692118 -0.1580604617 +8.5 1 641.6199025 0.2731219637 0.00009119724775 -0.0261686794 +9 1 1030.914723 0.2453117866 0.00005363701638 0.1043145752 +9.5 1 1658.453078 0.1612644308 0.00003160203411 0.2031798994 +10 1 2670.988304 0.04347274617 0.00001864877345 0.2490154242 +0.5 2 0.03190614918 0.03060402346 7.550183551 -5.441370837 +1 2 0.1357476698 0.1149034849 1.624838899 -1.650682607 +1.5 2 0.3378346183 0.2320876721 0.5836559633 -0.9321937598 +2 2 0.6889484477 0.3528340286 0.2537597546 -0.6174081042 +2.5 2 1.276466148 0.4460590584 0.1214602063 -0.3813358492 +3 2 2.245212441 0.4860912606 0.06151045847 -0.1604003935 +3.5 2 3.832012048 0.4586291842 0.0323071217 0.04537143773 +4 2 6.422189375 0.3641281459 0.01740142553 0.2159035946 +4.5 2 10.6415173 0.2178489837 0.009545677203 0.3284815967 +5 2 17.50561497 0.04656511628 0.005308943712 0.3676628826 +5.5 2 28.66258529 -0.1173154816 0.002984370024 0.3308412333 +6 2 46.78709472 -0.24287321 0.001691967567 0.2298579025 +6.5 2 76.22054724 -0.3074303906 0.0009658992748 0.08890665832 +7 2 124.0113105 -0.3014172201 0.0005545621667 -0.06052660947 +7.5 2 201.6054807 -0.2302734105 0.0003199235871 -0.1864142223 +8 2 327.5958315 -0.1129917204 0.0001853130082 -0.2630366048 +8.5 2 532.1925382 0.02232473961 0.0001077157423 -0.2763624417 +9 2 864.496194 0.1448473415 0.00006280064993 -0.2267556816 +9.5 2 1404.332974 0.2278791542 0.00003671094477 -0.1284359105 +10 2 2281.518968 0.2546303137 0.00002150981701 -0.005868082442 +0.5 3 0.002645111969 0.002563729995 62.05790953 -42.0594943 +1 3 0.02216842492 0.01956335398 7.101262825 -5.821517606 +1.5 3 0.08077411302 0.06096395114 1.833803702 -2.073541399 +2 3 0.2127399592 0.1289432495 0.6473853909 -1.127783777 +2.5 3 0.4743704088 0.216600391 0.2682271464 -0.7560554968 +3 3 0.9597536295 0.3090627223 0.1221703758 -0.5385416161 +3.5 3 1.826392582 0.3867701117 0.05916181773 -0.3583353462 +4 3 3.337275778 0.4301714739 0.02988492442 -0.182022116 +4.5 3 5.930096266 0.424703973 0.01556314131 -0.009013681594 +5 3 10.33115017 0.3648312306 0.008291768415 0.1462671627 +5.5 3 17.74264804 0.2561178651 0.004496019935 0.264370045 +6 3 30.1505403 0.1147683848 0.002471898096 0.328248946 +6.5 3 50.83005863 -0.03534663129 0.001374293952 0.3288030637 +7 3 85.17548684 -0.167555588 0.0007710751536 0.268080603 +7.5 3 142.0614424 -0.2580609132 0.0004359233032 0.1597075919 +8 3 236.075221 -0.2911322071 0.0002480257159 0.02654215932 +8.5 3 391.1763552 -0.2626162039 0.0001418870088 -0.1038842343 +9 3 646.6941919 -0.1809351903 0.00008154841635 -0.2050948781 +9.5 3 1067.154983 -0.06531531322 0.00004705927401 -0.2572581775 +10 3 1758.380717 0.05837937931 0.00002725270026 -0.2513626572