Brunobkr/llama.cpp_AlgMor24_github

1import {epsilon, splitter, resulterrbound, estimate, vec, sum, sum_three, scale, negate} from './util.js';2 3const isperrboundA = (16 + 224 * epsilon) * epsilon;4const isperrboundB = (5 + 72 * epsilon) * epsilon;5const isperrboundC = (71 + 1408 * epsilon) * epsilon * epsilon;6 7const ab = vec(4);8const bc = vec(4);9const cd = vec(4);10const de = vec(4);11const ea = vec(4);12const ac = vec(4);13const bd = vec(4);14const ce = vec(4);15const da = vec(4);16const eb = vec(4);17 18const abc = vec(24);19const bcd = vec(24);20const cde = vec(24);21const dea = vec(24);22const eab = vec(24);23const abd = vec(24);24const bce = vec(24);25const cda = vec(24);26const deb = vec(24);27const eac = vec(24);28 29const adet = vec(1152);30const bdet = vec(1152);31const cdet = vec(1152);32const ddet = vec(1152);33const edet = vec(1152);34const abdet = vec(2304);35const cddet = vec(2304);36const cdedet = vec(3456);37const deter = vec(5760);38 39const _8 = vec(8);40const _8b = vec(8);41const _8c = vec(8);42const _16 = vec(16);43const _24 = vec(24);44const _48 = vec(48);45const _48b = vec(48);46const _96 = vec(96);47const _192 = vec(192);48const _384x = vec(384);49const _384y = vec(384);50const _384z = vec(384);51const _768 = vec(768);52 53function sum_three_scale(a, b, c, az, bz, cz, out) {54    return sum_three(55        scale(4, a, az, _8), _8,56        scale(4, b, bz, _8b), _8b,57        scale(4, c, cz, _8c), _8c, _16, out);58}59 60function liftexact(alen, a, blen, b, clen, c, dlen, d, x, y, z, out) {61    const len = sum(62        sum(alen, a, blen, b, _48), _48,63        negate(sum(clen, c, dlen, d, _48b), _48b), _48b, _96);64 65    return sum_three(66        scale(scale(len, _96, x, _192), _192, x, _384x), _384x,67        scale(scale(len, _96, y, _192), _192, y, _384y), _384y,68        scale(scale(len, _96, z, _192), _192, z, _384z), _384z, _768, out);69}70 71function insphereexact(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez) {72    let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0, u3;73 74    s1 = ax * by;75    c = splitter * ax;76    ahi = c - (c - ax);77    alo = ax - ahi;78    c = splitter * by;79    bhi = c - (c - by);80    blo = by - bhi;81    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);82    t1 = bx * ay;83    c = splitter * bx;84    ahi = c - (c - bx);85    alo = bx - ahi;86    c = splitter * ay;87    bhi = c - (c - ay);88    blo = ay - bhi;89    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);90    _i = s0 - t0;91    bvirt = s0 - _i;92    ab[0] = s0 - (_i + bvirt) + (bvirt - t0);93    _j = s1 + _i;94    bvirt = _j - s1;95    _0 = s1 - (_j - bvirt) + (_i - bvirt);96    _i = _0 - t1;97    bvirt = _0 - _i;98    ab[1] = _0 - (_i + bvirt) + (bvirt - t1);99    u3 = _j + _i;100    bvirt = u3 - _j;101    ab[2] = _j - (u3 - bvirt) + (_i - bvirt);102    ab[3] = u3;103    s1 = bx * cy;104    c = splitter * bx;105    ahi = c - (c - bx);106    alo = bx - ahi;107    c = splitter * cy;108    bhi = c - (c - cy);109    blo = cy - bhi;110    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);111    t1 = cx * by;112    c = splitter * cx;113    ahi = c - (c - cx);114    alo = cx - ahi;115    c = splitter * by;116    bhi = c - (c - by);117    blo = by - bhi;118    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);119    _i = s0 - t0;120    bvirt = s0 - _i;121    bc[0] = s0 - (_i + bvirt) + (bvirt - t0);122    _j = s1 + _i;123    bvirt = _j - s1;124    _0 = s1 - (_j - bvirt) + (_i - bvirt);125    _i = _0 - t1;126    bvirt = _0 - _i;127    bc[1] = _0 - (_i + bvirt) + (bvirt - t1);128    u3 = _j + _i;129    bvirt = u3 - _j;130    bc[2] = _j - (u3 - bvirt) + (_i - bvirt);131    bc[3] = u3;132    s1 = cx * dy;133    c = splitter * cx;134    ahi = c - (c - cx);135    alo = cx - ahi;136    c = splitter * dy;137    bhi = c - (c - dy);138    blo = dy - bhi;139    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);140    t1 = dx * cy;141    c = splitter * dx;142    ahi = c - (c - dx);143    alo = dx - ahi;144    c = splitter * cy;145    bhi = c - (c - cy);146    blo = cy - bhi;147    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);148    _i = s0 - t0;149    bvirt = s0 - _i;150    cd[0] = s0 - (_i + bvirt) + (bvirt - t0);151    _j = s1 + _i;152    bvirt = _j - s1;153    _0 = s1 - (_j - bvirt) + (_i - bvirt);154    _i = _0 - t1;155    bvirt = _0 - _i;156    cd[1] = _0 - (_i + bvirt) + (bvirt - t1);157    u3 = _j + _i;158    bvirt = u3 - _j;159    cd[2] = _j - (u3 - bvirt) + (_i - bvirt);160    cd[3] = u3;161    s1 = dx * ey;162    c = splitter * dx;163    ahi = c - (c - dx);164    alo = dx - ahi;165    c = splitter * ey;166    bhi = c - (c - ey);167    blo = ey - bhi;168    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);169    t1 = ex * dy;170    c = splitter * ex;171    ahi = c - (c - ex);172    alo = ex - ahi;173    c = splitter * dy;174    bhi = c - (c - dy);175    blo = dy - bhi;176    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);177    _i = s0 - t0;178    bvirt = s0 - _i;179    de[0] = s0 - (_i + bvirt) + (bvirt - t0);180    _j = s1 + _i;181    bvirt = _j - s1;182    _0 = s1 - (_j - bvirt) + (_i - bvirt);183    _i = _0 - t1;184    bvirt = _0 - _i;185    de[1] = _0 - (_i + bvirt) + (bvirt - t1);186    u3 = _j + _i;187    bvirt = u3 - _j;188    de[2] = _j - (u3 - bvirt) + (_i - bvirt);189    de[3] = u3;190    s1 = ex * ay;191    c = splitter * ex;192    ahi = c - (c - ex);193    alo = ex - ahi;194    c = splitter * ay;195    bhi = c - (c - ay);196    blo = ay - bhi;197    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);198    t1 = ax * ey;199    c = splitter * ax;200    ahi = c - (c - ax);201    alo = ax - ahi;202    c = splitter * ey;203    bhi = c - (c - ey);204    blo = ey - bhi;205    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);206    _i = s0 - t0;207    bvirt = s0 - _i;208    ea[0] = s0 - (_i + bvirt) + (bvirt - t0);209    _j = s1 + _i;210    bvirt = _j - s1;211    _0 = s1 - (_j - bvirt) + (_i - bvirt);212    _i = _0 - t1;213    bvirt = _0 - _i;214    ea[1] = _0 - (_i + bvirt) + (bvirt - t1);215    u3 = _j + _i;216    bvirt = u3 - _j;217    ea[2] = _j - (u3 - bvirt) + (_i - bvirt);218    ea[3] = u3;219    s1 = ax * cy;220    c = splitter * ax;221    ahi = c - (c - ax);222    alo = ax - ahi;223    c = splitter * cy;224    bhi = c - (c - cy);225    blo = cy - bhi;226    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);227    t1 = cx * ay;228    c = splitter * cx;229    ahi = c - (c - cx);230    alo = cx - ahi;231    c = splitter * ay;232    bhi = c - (c - ay);233    blo = ay - bhi;234    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);235    _i = s0 - t0;236    bvirt = s0 - _i;237    ac[0] = s0 - (_i + bvirt) + (bvirt - t0);238    _j = s1 + _i;239    bvirt = _j - s1;240    _0 = s1 - (_j - bvirt) + (_i - bvirt);241    _i = _0 - t1;242    bvirt = _0 - _i;243    ac[1] = _0 - (_i + bvirt) + (bvirt - t1);244    u3 = _j + _i;245    bvirt = u3 - _j;246    ac[2] = _j - (u3 - bvirt) + (_i - bvirt);247    ac[3] = u3;248    s1 = bx * dy;249    c = splitter * bx;250    ahi = c - (c - bx);251    alo = bx - ahi;252    c = splitter * dy;253    bhi = c - (c - dy);254    blo = dy - bhi;255    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);256    t1 = dx * by;257    c = splitter * dx;258    ahi = c - (c - dx);259    alo = dx - ahi;260    c = splitter * by;261    bhi = c - (c - by);262    blo = by - bhi;263    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);264    _i = s0 - t0;265    bvirt = s0 - _i;266    bd[0] = s0 - (_i + bvirt) + (bvirt - t0);267    _j = s1 + _i;268    bvirt = _j - s1;269    _0 = s1 - (_j - bvirt) + (_i - bvirt);270    _i = _0 - t1;271    bvirt = _0 - _i;272    bd[1] = _0 - (_i + bvirt) + (bvirt - t1);273    u3 = _j + _i;274    bvirt = u3 - _j;275    bd[2] = _j - (u3 - bvirt) + (_i - bvirt);276    bd[3] = u3;277    s1 = cx * ey;278    c = splitter * cx;279    ahi = c - (c - cx);280    alo = cx - ahi;281    c = splitter * ey;282    bhi = c - (c - ey);283    blo = ey - bhi;284    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);285    t1 = ex * cy;286    c = splitter * ex;287    ahi = c - (c - ex);288    alo = ex - ahi;289    c = splitter * cy;290    bhi = c - (c - cy);291    blo = cy - bhi;292    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);293    _i = s0 - t0;294    bvirt = s0 - _i;295    ce[0] = s0 - (_i + bvirt) + (bvirt - t0);296    _j = s1 + _i;297    bvirt = _j - s1;298    _0 = s1 - (_j - bvirt) + (_i - bvirt);299    _i = _0 - t1;300    bvirt = _0 - _i;301    ce[1] = _0 - (_i + bvirt) + (bvirt - t1);302    u3 = _j + _i;303    bvirt = u3 - _j;304    ce[2] = _j - (u3 - bvirt) + (_i - bvirt);305    ce[3] = u3;306    s1 = dx * ay;307    c = splitter * dx;308    ahi = c - (c - dx);309    alo = dx - ahi;310    c = splitter * ay;311    bhi = c - (c - ay);312    blo = ay - bhi;313    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);314    t1 = ax * dy;315    c = splitter * ax;316    ahi = c - (c - ax);317    alo = ax - ahi;318    c = splitter * dy;319    bhi = c - (c - dy);320    blo = dy - bhi;321    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);322    _i = s0 - t0;323    bvirt = s0 - _i;324    da[0] = s0 - (_i + bvirt) + (bvirt - t0);325    _j = s1 + _i;326    bvirt = _j - s1;327    _0 = s1 - (_j - bvirt) + (_i - bvirt);328    _i = _0 - t1;329    bvirt = _0 - _i;330    da[1] = _0 - (_i + bvirt) + (bvirt - t1);331    u3 = _j + _i;332    bvirt = u3 - _j;333    da[2] = _j - (u3 - bvirt) + (_i - bvirt);334    da[3] = u3;335    s1 = ex * by;336    c = splitter * ex;337    ahi = c - (c - ex);338    alo = ex - ahi;339    c = splitter * by;340    bhi = c - (c - by);341    blo = by - bhi;342    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);343    t1 = bx * ey;344    c = splitter * bx;345    ahi = c - (c - bx);346    alo = bx - ahi;347    c = splitter * ey;348    bhi = c - (c - ey);349    blo = ey - bhi;350    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);351    _i = s0 - t0;352    bvirt = s0 - _i;353    eb[0] = s0 - (_i + bvirt) + (bvirt - t0);354    _j = s1 + _i;355    bvirt = _j - s1;356    _0 = s1 - (_j - bvirt) + (_i - bvirt);357    _i = _0 - t1;358    bvirt = _0 - _i;359    eb[1] = _0 - (_i + bvirt) + (bvirt - t1);360    u3 = _j + _i;361    bvirt = u3 - _j;362    eb[2] = _j - (u3 - bvirt) + (_i - bvirt);363    eb[3] = u3;364 365    const abclen = sum_three_scale(ab, bc, ac, cz, az, -bz, abc);366    const bcdlen = sum_three_scale(bc, cd, bd, dz, bz, -cz, bcd);367    const cdelen = sum_three_scale(cd, de, ce, ez, cz, -dz, cde);368    const dealen = sum_three_scale(de, ea, da, az, dz, -ez, dea);369    const eablen = sum_three_scale(ea, ab, eb, bz, ez, -az, eab);370    const abdlen = sum_three_scale(ab, bd, da, dz, az, bz, abd);371    const bcelen = sum_three_scale(bc, ce, eb, ez, bz, cz, bce);372    const cdalen = sum_three_scale(cd, da, ac, az, cz, dz, cda);373    const deblen = sum_three_scale(de, eb, bd, bz, dz, ez, deb);374    const eaclen = sum_three_scale(ea, ac, ce, cz, ez, az, eac);375 376    const deterlen = sum_three(377        liftexact(cdelen, cde, bcelen, bce, deblen, deb, bcdlen, bcd, ax, ay, az, adet), adet,378        liftexact(dealen, dea, cdalen, cda, eaclen, eac, cdelen, cde, bx, by, bz, bdet), bdet,379        sum_three(380            liftexact(eablen, eab, deblen, deb, abdlen, abd, dealen, dea, cx, cy, cz, cdet), cdet,381            liftexact(abclen, abc, eaclen, eac, bcelen, bce, eablen, eab, dx, dy, dz, ddet), ddet,382            liftexact(bcdlen, bcd, abdlen, abd, cdalen, cda, abclen, abc, ex, ey, ez, edet), edet, cddet, cdedet), cdedet, abdet, deter);383 384    return deter[deterlen - 1];385}386 387const xdet = vec(96);388const ydet = vec(96);389const zdet = vec(96);390const fin = vec(1152);391 392function liftadapt(a, b, c, az, bz, cz, x, y, z, out) {393    const len = sum_three_scale(a, b, c, az, bz, cz, _24);394    return sum_three(395        scale(scale(len, _24, x, _48), _48, x, xdet), xdet,396        scale(scale(len, _24, y, _48), _48, y, ydet), ydet,397        scale(scale(len, _24, z, _48), _48, z, zdet), zdet, _192, out);398}399 400function insphereadapt(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez, permanent) {401    let ab3, bc3, cd3, da3, ac3, bd3;402 403    let aextail, bextail, cextail, dextail;404    let aeytail, beytail, ceytail, deytail;405    let aeztail, beztail, ceztail, deztail;406 407    let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0;408 409    const aex = ax - ex;410    const bex = bx - ex;411    const cex = cx - ex;412    const dex = dx - ex;413    const aey = ay - ey;414    const bey = by - ey;415    const cey = cy - ey;416    const dey = dy - ey;417    const aez = az - ez;418    const bez = bz - ez;419    const cez = cz - ez;420    const dez = dz - ez;421 422    s1 = aex * bey;423    c = splitter * aex;424    ahi = c - (c - aex);425    alo = aex - ahi;426    c = splitter * bey;427    bhi = c - (c - bey);428    blo = bey - bhi;429    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);430    t1 = bex * aey;431    c = splitter * bex;432    ahi = c - (c - bex);433    alo = bex - ahi;434    c = splitter * aey;435    bhi = c - (c - aey);436    blo = aey - bhi;437    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);438    _i = s0 - t0;439    bvirt = s0 - _i;440    ab[0] = s0 - (_i + bvirt) + (bvirt - t0);441    _j = s1 + _i;442    bvirt = _j - s1;443    _0 = s1 - (_j - bvirt) + (_i - bvirt);444    _i = _0 - t1;445    bvirt = _0 - _i;446    ab[1] = _0 - (_i + bvirt) + (bvirt - t1);447    ab3 = _j + _i;448    bvirt = ab3 - _j;449    ab[2] = _j - (ab3 - bvirt) + (_i - bvirt);450    ab[3] = ab3;451    s1 = bex * cey;452    c = splitter * bex;453    ahi = c - (c - bex);454    alo = bex - ahi;455    c = splitter * cey;456    bhi = c - (c - cey);457    blo = cey - bhi;458    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);459    t1 = cex * bey;460    c = splitter * cex;461    ahi = c - (c - cex);462    alo = cex - ahi;463    c = splitter * bey;464    bhi = c - (c - bey);465    blo = bey - bhi;466    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);467    _i = s0 - t0;468    bvirt = s0 - _i;469    bc[0] = s0 - (_i + bvirt) + (bvirt - t0);470    _j = s1 + _i;471    bvirt = _j - s1;472    _0 = s1 - (_j - bvirt) + (_i - bvirt);473    _i = _0 - t1;474    bvirt = _0 - _i;475    bc[1] = _0 - (_i + bvirt) + (bvirt - t1);476    bc3 = _j + _i;477    bvirt = bc3 - _j;478    bc[2] = _j - (bc3 - bvirt) + (_i - bvirt);479    bc[3] = bc3;480    s1 = cex * dey;481    c = splitter * cex;482    ahi = c - (c - cex);483    alo = cex - ahi;484    c = splitter * dey;485    bhi = c - (c - dey);486    blo = dey - bhi;487    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);488    t1 = dex * cey;489    c = splitter * dex;490    ahi = c - (c - dex);491    alo = dex - ahi;492    c = splitter * cey;493    bhi = c - (c - cey);494    blo = cey - bhi;495    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);496    _i = s0 - t0;497    bvirt = s0 - _i;498    cd[0] = s0 - (_i + bvirt) + (bvirt - t0);499    _j = s1 + _i;500    bvirt = _j - s1;501    _0 = s1 - (_j - bvirt) + (_i - bvirt);502    _i = _0 - t1;503    bvirt = _0 - _i;504    cd[1] = _0 - (_i + bvirt) + (bvirt - t1);505    cd3 = _j + _i;506    bvirt = cd3 - _j;507    cd[2] = _j - (cd3 - bvirt) + (_i - bvirt);508    cd[3] = cd3;509    s1 = dex * aey;510    c = splitter * dex;511    ahi = c - (c - dex);512    alo = dex - ahi;513    c = splitter * aey;514    bhi = c - (c - aey);515    blo = aey - bhi;516    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);517    t1 = aex * dey;518    c = splitter * aex;519    ahi = c - (c - aex);520    alo = aex - ahi;521    c = splitter * dey;522    bhi = c - (c - dey);523    blo = dey - bhi;524    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);525    _i = s0 - t0;526    bvirt = s0 - _i;527    da[0] = s0 - (_i + bvirt) + (bvirt - t0);528    _j = s1 + _i;529    bvirt = _j - s1;530    _0 = s1 - (_j - bvirt) + (_i - bvirt);531    _i = _0 - t1;532    bvirt = _0 - _i;533    da[1] = _0 - (_i + bvirt) + (bvirt - t1);534    da3 = _j + _i;535    bvirt = da3 - _j;536    da[2] = _j - (da3 - bvirt) + (_i - bvirt);537    da[3] = da3;538    s1 = aex * cey;539    c = splitter * aex;540    ahi = c - (c - aex);541    alo = aex - ahi;542    c = splitter * cey;543    bhi = c - (c - cey);544    blo = cey - bhi;545    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);546    t1 = cex * aey;547    c = splitter * cex;548    ahi = c - (c - cex);549    alo = cex - ahi;550    c = splitter * aey;551    bhi = c - (c - aey);552    blo = aey - bhi;553    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);554    _i = s0 - t0;555    bvirt = s0 - _i;556    ac[0] = s0 - (_i + bvirt) + (bvirt - t0);557    _j = s1 + _i;558    bvirt = _j - s1;559    _0 = s1 - (_j - bvirt) + (_i - bvirt);560    _i = _0 - t1;561    bvirt = _0 - _i;562    ac[1] = _0 - (_i + bvirt) + (bvirt - t1);563    ac3 = _j + _i;564    bvirt = ac3 - _j;565    ac[2] = _j - (ac3 - bvirt) + (_i - bvirt);566    ac[3] = ac3;567    s1 = bex * dey;568    c = splitter * bex;569    ahi = c - (c - bex);570    alo = bex - ahi;571    c = splitter * dey;572    bhi = c - (c - dey);573    blo = dey - bhi;574    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);575    t1 = dex * bey;576    c = splitter * dex;577    ahi = c - (c - dex);578    alo = dex - ahi;579    c = splitter * bey;580    bhi = c - (c - bey);581    blo = bey - bhi;582    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);583    _i = s0 - t0;584    bvirt = s0 - _i;585    bd[0] = s0 - (_i + bvirt) + (bvirt - t0);586    _j = s1 + _i;587    bvirt = _j - s1;588    _0 = s1 - (_j - bvirt) + (_i - bvirt);589    _i = _0 - t1;590    bvirt = _0 - _i;591    bd[1] = _0 - (_i + bvirt) + (bvirt - t1);592    bd3 = _j + _i;593    bvirt = bd3 - _j;594    bd[2] = _j - (bd3 - bvirt) + (_i - bvirt);595    bd[3] = bd3;596 597    const finlen = sum(598        sum(599            negate(liftadapt(bc, cd, bd, dez, bez, -cez, aex, aey, aez, adet), adet), adet,600            liftadapt(cd, da, ac, aez, cez, dez, bex, bey, bez, bdet), bdet, abdet), abdet,601        sum(602            negate(liftadapt(da, ab, bd, bez, dez, aez, cex, cey, cez, cdet), cdet), cdet,603            liftadapt(ab, bc, ac, cez, aez, -bez, dex, dey, dez, ddet), ddet, cddet), cddet, fin);604 605    let det = estimate(finlen, fin);606    let errbound = isperrboundB * permanent;607    if (det >= errbound || -det >= errbound) {608        return det;609    }610 611    bvirt = ax - aex;612    aextail = ax - (aex + bvirt) + (bvirt - ex);613    bvirt = ay - aey;614    aeytail = ay - (aey + bvirt) + (bvirt - ey);615    bvirt = az - aez;616    aeztail = az - (aez + bvirt) + (bvirt - ez);617    bvirt = bx - bex;618    bextail = bx - (bex + bvirt) + (bvirt - ex);619    bvirt = by - bey;620    beytail = by - (bey + bvirt) + (bvirt - ey);621    bvirt = bz - bez;622    beztail = bz - (bez + bvirt) + (bvirt - ez);623    bvirt = cx - cex;624    cextail = cx - (cex + bvirt) + (bvirt - ex);625    bvirt = cy - cey;626    ceytail = cy - (cey + bvirt) + (bvirt - ey);627    bvirt = cz - cez;628    ceztail = cz - (cez + bvirt) + (bvirt - ez);629    bvirt = dx - dex;630    dextail = dx - (dex + bvirt) + (bvirt - ex);631    bvirt = dy - dey;632    deytail = dy - (dey + bvirt) + (bvirt - ey);633    bvirt = dz - dez;634    deztail = dz - (dez + bvirt) + (bvirt - ez);635    if (aextail === 0 && aeytail === 0 && aeztail === 0 &&636        bextail === 0 && beytail === 0 && beztail === 0 &&637        cextail === 0 && ceytail === 0 && ceztail === 0 &&638        dextail === 0 && deytail === 0 && deztail === 0) {639        return det;640    }641 642    errbound = isperrboundC * permanent + resulterrbound * Math.abs(det);643 644    const abeps = (aex * beytail + bey * aextail) - (aey * bextail + bex * aeytail);645    const bceps = (bex * ceytail + cey * bextail) - (bey * cextail + cex * beytail);646    const cdeps = (cex * deytail + dey * cextail) - (cey * dextail + dex * ceytail);647    const daeps = (dex * aeytail + aey * dextail) - (dey * aextail + aex * deytail);648    const aceps = (aex * ceytail + cey * aextail) - (aey * cextail + cex * aeytail);649    const bdeps = (bex * deytail + dey * bextail) - (bey * dextail + dex * beytail);650    det +=651        (((bex * bex + bey * bey + bez * bez) * ((cez * daeps + dez * aceps + aez * cdeps) +652        (ceztail * da3 + deztail * ac3 + aeztail * cd3)) + (dex * dex + dey * dey + dez * dez) *653        ((aez * bceps - bez * aceps + cez * abeps) + (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) -654        ((aex * aex + aey * aey + aez * aez) * ((bez * cdeps - cez * bdeps + dez * bceps) +655        (beztail * cd3 - ceztail * bd3 + deztail * bc3)) + (cex * cex + cey * cey + cez * cez) *656        ((dez * abeps + aez * bdeps + bez * daeps) + (deztail * ab3 + aeztail * bd3 + beztail * da3)))) +657        2 * (((bex * bextail + bey * beytail + bez * beztail) * (cez * da3 + dez * ac3 + aez * cd3) +658        (dex * dextail + dey * deytail + dez * deztail) * (aez * bc3 - bez * ac3 + cez * ab3)) -659        ((aex * aextail + aey * aeytail + aez * aeztail) * (bez * cd3 - cez * bd3 + dez * bc3) +660        (cex * cextail + cey * ceytail + cez * ceztail) * (dez * ab3 + aez * bd3 + bez * da3)));661 662    if (det >= errbound || -det >= errbound) {663        return det;664    }665 666    return insphereexact(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez);667}668 669export function insphere(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez) {670    const aex = ax - ex;671    const bex = bx - ex;672    const cex = cx - ex;673    const dex = dx - ex;674    const aey = ay - ey;675    const bey = by - ey;676    const cey = cy - ey;677    const dey = dy - ey;678    const aez = az - ez;679    const bez = bz - ez;680    const cez = cz - ez;681    const dez = dz - ez;682 683    const aexbey = aex * bey;684    const bexaey = bex * aey;685    const ab = aexbey - bexaey;686    const bexcey = bex * cey;687    const cexbey = cex * bey;688    const bc = bexcey - cexbey;689    const cexdey = cex * dey;690    const dexcey = dex * cey;691    const cd = cexdey - dexcey;692    const dexaey = dex * aey;693    const aexdey = aex * dey;694    const da = dexaey - aexdey;695    const aexcey = aex * cey;696    const cexaey = cex * aey;697    const ac = aexcey - cexaey;698    const bexdey = bex * dey;699    const dexbey = dex * bey;700    const bd = bexdey - dexbey;701 702    const alift = aex * aex + aey * aey + aez * aez;703    const blift = bex * bex + bey * bey + bez * bez;704    const clift = cex * cex + cey * cey + cez * cez;705    const dlift = dex * dex + dey * dey + dez * dez;706 707    const det =708        (clift * (dez * ab + aez * bd + bez * da) - dlift * (aez * bc - bez * ac + cez * ab)) +709        (alift * (bez * cd - cez * bd + dez * bc) - blift * (cez * da + dez * ac + aez * cd));710 711    const aezplus = Math.abs(aez);712    const bezplus = Math.abs(bez);713    const cezplus = Math.abs(cez);714    const dezplus = Math.abs(dez);715    const aexbeyplus = Math.abs(aexbey) + Math.abs(bexaey);716    const bexceyplus = Math.abs(bexcey) + Math.abs(cexbey);717    const cexdeyplus = Math.abs(cexdey) + Math.abs(dexcey);718    const dexaeyplus = Math.abs(dexaey) + Math.abs(aexdey);719    const aexceyplus = Math.abs(aexcey) + Math.abs(cexaey);720    const bexdeyplus = Math.abs(bexdey) + Math.abs(dexbey);721    const permanent =722        (cexdeyplus * bezplus + bexdeyplus * cezplus + bexceyplus * dezplus) * alift +723        (dexaeyplus * cezplus + aexceyplus * dezplus + cexdeyplus * aezplus) * blift +724        (aexbeyplus * dezplus + bexdeyplus * aezplus + dexaeyplus * bezplus) * clift +725        (bexceyplus * aezplus + aexceyplus * bezplus + aexbeyplus * cezplus) * dlift;726 727    const errbound = isperrboundA * permanent;728    if (det > errbound || -det > errbound) {729        return det;730    }731    return -insphereadapt(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez, permanent);732}733 734export function inspherefast(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez) {735    const aex = ax - ex;736    const bex = bx - ex;737    const cex = cx - ex;738    const dex = dx - ex;739    const aey = ay - ey;740    const bey = by - ey;741    const cey = cy - ey;742    const dey = dy - ey;743    const aez = az - ez;744    const bez = bz - ez;745    const cez = cz - ez;746    const dez = dz - ez;747 748    const ab = aex * bey - bex * aey;749    const bc = bex * cey - cex * bey;750    const cd = cex * dey - dex * cey;751    const da = dex * aey - aex * dey;752    const ac = aex * cey - cex * aey;753    const bd = bex * dey - dex * bey;754 755    const abc = aez * bc - bez * ac + cez * ab;756    const bcd = bez * cd - cez * bd + dez * bc;757    const cda = cez * da + dez * ac + aez * cd;758    const dab = dez * ab + aez * bd + bez * da;759 760    const alift = aex * aex + aey * aey + aez * aez;761    const blift = bex * bex + bey * bey + bez * bez;762    const clift = cex * cex + cey * cey + cez * cez;763    const dlift = dex * dex + dey * dey + dez * dez;764 765    return (clift * dab - dlift * abc) + (alift * bcd - blift * cda);766}767