123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138 |
- export const epsilon = 1.1102230246251565e-16;
- export const splitter = 134217729;
- export const resulterrbound = (3 + 8 * epsilon) * epsilon;
-
- // fast_expansion_sum_zeroelim routine from oritinal code
- export function sum(elen, e, flen, f, h) {
- let Q, Qnew, hh, bvirt;
- let enow = e[0];
- let fnow = f[0];
- let eindex = 0;
- let findex = 0;
- if ((fnow > enow) === (fnow > -enow)) {
- Q = enow;
- enow = e[++eindex];
- } else {
- Q = fnow;
- fnow = f[++findex];
- }
- let hindex = 0;
- if (eindex < elen && findex < flen) {
- if ((fnow > enow) === (fnow > -enow)) {
- Qnew = enow + Q;
- hh = Q - (Qnew - enow);
- enow = e[++eindex];
- } else {
- Qnew = fnow + Q;
- hh = Q - (Qnew - fnow);
- fnow = f[++findex];
- }
- Q = Qnew;
- if (hh !== 0) {
- h[hindex++] = hh;
- }
- while (eindex < elen && findex < flen) {
- if ((fnow > enow) === (fnow > -enow)) {
- Qnew = Q + enow;
- bvirt = Qnew - Q;
- hh = Q - (Qnew - bvirt) + (enow - bvirt);
- enow = e[++eindex];
- } else {
- Qnew = Q + fnow;
- bvirt = Qnew - Q;
- hh = Q - (Qnew - bvirt) + (fnow - bvirt);
- fnow = f[++findex];
- }
- Q = Qnew;
- if (hh !== 0) {
- h[hindex++] = hh;
- }
- }
- }
- while (eindex < elen) {
- Qnew = Q + enow;
- bvirt = Qnew - Q;
- hh = Q - (Qnew - bvirt) + (enow - bvirt);
- enow = e[++eindex];
- Q = Qnew;
- if (hh !== 0) {
- h[hindex++] = hh;
- }
- }
- while (findex < flen) {
- Qnew = Q + fnow;
- bvirt = Qnew - Q;
- hh = Q - (Qnew - bvirt) + (fnow - bvirt);
- fnow = f[++findex];
- Q = Qnew;
- if (hh !== 0) {
- h[hindex++] = hh;
- }
- }
- if (Q !== 0 || hindex === 0) {
- h[hindex++] = Q;
- }
- return hindex;
- }
-
- export function sum_three(alen, a, blen, b, clen, c, tmp, out) {
- return sum(sum(alen, a, blen, b, tmp), tmp, clen, c, out);
- }
-
- // scale_expansion_zeroelim routine from oritinal code
- export function scale(elen, e, b, h) {
- let Q, sum, hh, product1, product0;
- let bvirt, c, ahi, alo, bhi, blo;
-
- c = splitter * b;
- bhi = c - (c - b);
- blo = b - bhi;
- let enow = e[0];
- Q = enow * b;
- c = splitter * enow;
- ahi = c - (c - enow);
- alo = enow - ahi;
- hh = alo * blo - (Q - ahi * bhi - alo * bhi - ahi * blo);
- let hindex = 0;
- if (hh !== 0) {
- h[hindex++] = hh;
- }
- for (let i = 1; i < elen; i++) {
- enow = e[i];
- product1 = enow * b;
- c = splitter * enow;
- ahi = c - (c - enow);
- alo = enow - ahi;
- product0 = alo * blo - (product1 - ahi * bhi - alo * bhi - ahi * blo);
- sum = Q + product0;
- bvirt = sum - Q;
- hh = Q - (sum - bvirt) + (product0 - bvirt);
- if (hh !== 0) {
- h[hindex++] = hh;
- }
- Q = product1 + sum;
- hh = sum - (Q - product1);
- if (hh !== 0) {
- h[hindex++] = hh;
- }
- }
- if (Q !== 0 || hindex === 0) {
- h[hindex++] = Q;
- }
- return hindex;
- }
-
- export function negate(elen, e) {
- for (let i = 0; i < elen; i++) e[i] = -e[i];
- return elen;
- }
-
- export function estimate(elen, e) {
- let Q = e[0];
- for (let i = 1; i < elen; i++) Q += e[i];
- return Q;
- }
-
- export function vec(n) {
- return new Float64Array(n);
- }
|