source: node_modules/robust-predicates/umd/predicates.js@ e4c61dd

Last change on this file since e4c61dd was e4c61dd, checked in by istevanoska <ilinastevanoska@…>, 6 months ago

Prototype 1.1

  • Property mode set to 100644
File size: 72.6 KB
Line 
1(function (global, factory) {
2typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports) :
3typeof define === 'function' && define.amd ? define(['exports'], factory) :
4(global = typeof globalThis !== 'undefined' ? globalThis : global || self, factory(global.predicates = {}));
5})(this, (function (exports) { 'use strict';
6
7const epsilon = 1.1102230246251565e-16;
8const splitter = 134217729;
9const resulterrbound = (3 + 8 * epsilon) * epsilon;
10
11// fast_expansion_sum_zeroelim routine from oritinal code
12function sum(elen, e, flen, f, h) {
13 let Q, Qnew, hh, bvirt;
14 let enow = e[0];
15 let fnow = f[0];
16 let eindex = 0;
17 let findex = 0;
18 if ((fnow > enow) === (fnow > -enow)) {
19 Q = enow;
20 enow = e[++eindex];
21 } else {
22 Q = fnow;
23 fnow = f[++findex];
24 }
25 let hindex = 0;
26 if (eindex < elen && findex < flen) {
27 if ((fnow > enow) === (fnow > -enow)) {
28 Qnew = enow + Q;
29 hh = Q - (Qnew - enow);
30 enow = e[++eindex];
31 } else {
32 Qnew = fnow + Q;
33 hh = Q - (Qnew - fnow);
34 fnow = f[++findex];
35 }
36 Q = Qnew;
37 if (hh !== 0) {
38 h[hindex++] = hh;
39 }
40 while (eindex < elen && findex < flen) {
41 if ((fnow > enow) === (fnow > -enow)) {
42 Qnew = Q + enow;
43 bvirt = Qnew - Q;
44 hh = Q - (Qnew - bvirt) + (enow - bvirt);
45 enow = e[++eindex];
46 } else {
47 Qnew = Q + fnow;
48 bvirt = Qnew - Q;
49 hh = Q - (Qnew - bvirt) + (fnow - bvirt);
50 fnow = f[++findex];
51 }
52 Q = Qnew;
53 if (hh !== 0) {
54 h[hindex++] = hh;
55 }
56 }
57 }
58 while (eindex < elen) {
59 Qnew = Q + enow;
60 bvirt = Qnew - Q;
61 hh = Q - (Qnew - bvirt) + (enow - bvirt);
62 enow = e[++eindex];
63 Q = Qnew;
64 if (hh !== 0) {
65 h[hindex++] = hh;
66 }
67 }
68 while (findex < flen) {
69 Qnew = Q + fnow;
70 bvirt = Qnew - Q;
71 hh = Q - (Qnew - bvirt) + (fnow - bvirt);
72 fnow = f[++findex];
73 Q = Qnew;
74 if (hh !== 0) {
75 h[hindex++] = hh;
76 }
77 }
78 if (Q !== 0 || hindex === 0) {
79 h[hindex++] = Q;
80 }
81 return hindex;
82}
83
84function sum_three(alen, a, blen, b, clen, c, tmp, out) {
85 return sum(sum(alen, a, blen, b, tmp), tmp, clen, c, out);
86}
87
88// scale_expansion_zeroelim routine from oritinal code
89function scale(elen, e, b, h) {
90 let Q, sum, hh, product1, product0;
91 let bvirt, c, ahi, alo, bhi, blo;
92
93 c = splitter * b;
94 bhi = c - (c - b);
95 blo = b - bhi;
96 let enow = e[0];
97 Q = enow * b;
98 c = splitter * enow;
99 ahi = c - (c - enow);
100 alo = enow - ahi;
101 hh = alo * blo - (Q - ahi * bhi - alo * bhi - ahi * blo);
102 let hindex = 0;
103 if (hh !== 0) {
104 h[hindex++] = hh;
105 }
106 for (let i = 1; i < elen; i++) {
107 enow = e[i];
108 product1 = enow * b;
109 c = splitter * enow;
110 ahi = c - (c - enow);
111 alo = enow - ahi;
112 product0 = alo * blo - (product1 - ahi * bhi - alo * bhi - ahi * blo);
113 sum = Q + product0;
114 bvirt = sum - Q;
115 hh = Q - (sum - bvirt) + (product0 - bvirt);
116 if (hh !== 0) {
117 h[hindex++] = hh;
118 }
119 Q = product1 + sum;
120 hh = sum - (Q - product1);
121 if (hh !== 0) {
122 h[hindex++] = hh;
123 }
124 }
125 if (Q !== 0 || hindex === 0) {
126 h[hindex++] = Q;
127 }
128 return hindex;
129}
130
131function negate(elen, e) {
132 for (let i = 0; i < elen; i++) e[i] = -e[i];
133 return elen;
134}
135
136function estimate(elen, e) {
137 let Q = e[0];
138 for (let i = 1; i < elen; i++) Q += e[i];
139 return Q;
140}
141
142function vec(n) {
143 return new Float64Array(n);
144}
145
146const ccwerrboundA = (3 + 16 * epsilon) * epsilon;
147const ccwerrboundB = (2 + 12 * epsilon) * epsilon;
148const ccwerrboundC = (9 + 64 * epsilon) * epsilon * epsilon;
149
150const B = vec(4);
151const C1 = vec(8);
152const C2 = vec(12);
153const D = vec(16);
154const u$2 = vec(4);
155
156function orient2dadapt(ax, ay, bx, by, cx, cy, detsum) {
157 let acxtail, acytail, bcxtail, bcytail;
158 let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0, u3;
159
160 const acx = ax - cx;
161 const bcx = bx - cx;
162 const acy = ay - cy;
163 const bcy = by - cy;
164
165 s1 = acx * bcy;
166 c = splitter * acx;
167 ahi = c - (c - acx);
168 alo = acx - ahi;
169 c = splitter * bcy;
170 bhi = c - (c - bcy);
171 blo = bcy - bhi;
172 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
173 t1 = acy * bcx;
174 c = splitter * acy;
175 ahi = c - (c - acy);
176 alo = acy - ahi;
177 c = splitter * bcx;
178 bhi = c - (c - bcx);
179 blo = bcx - bhi;
180 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
181 _i = s0 - t0;
182 bvirt = s0 - _i;
183 B[0] = s0 - (_i + bvirt) + (bvirt - t0);
184 _j = s1 + _i;
185 bvirt = _j - s1;
186 _0 = s1 - (_j - bvirt) + (_i - bvirt);
187 _i = _0 - t1;
188 bvirt = _0 - _i;
189 B[1] = _0 - (_i + bvirt) + (bvirt - t1);
190 u3 = _j + _i;
191 bvirt = u3 - _j;
192 B[2] = _j - (u3 - bvirt) + (_i - bvirt);
193 B[3] = u3;
194
195 let det = estimate(4, B);
196 let errbound = ccwerrboundB * detsum;
197 if (det >= errbound || -det >= errbound) {
198 return det;
199 }
200
201 bvirt = ax - acx;
202 acxtail = ax - (acx + bvirt) + (bvirt - cx);
203 bvirt = bx - bcx;
204 bcxtail = bx - (bcx + bvirt) + (bvirt - cx);
205 bvirt = ay - acy;
206 acytail = ay - (acy + bvirt) + (bvirt - cy);
207 bvirt = by - bcy;
208 bcytail = by - (bcy + bvirt) + (bvirt - cy);
209
210 if (acxtail === 0 && acytail === 0 && bcxtail === 0 && bcytail === 0) {
211 return det;
212 }
213
214 errbound = ccwerrboundC * detsum + resulterrbound * Math.abs(det);
215 det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail);
216 if (det >= errbound || -det >= errbound) return det;
217
218 s1 = acxtail * bcy;
219 c = splitter * acxtail;
220 ahi = c - (c - acxtail);
221 alo = acxtail - ahi;
222 c = splitter * bcy;
223 bhi = c - (c - bcy);
224 blo = bcy - bhi;
225 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
226 t1 = acytail * bcx;
227 c = splitter * acytail;
228 ahi = c - (c - acytail);
229 alo = acytail - ahi;
230 c = splitter * bcx;
231 bhi = c - (c - bcx);
232 blo = bcx - bhi;
233 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
234 _i = s0 - t0;
235 bvirt = s0 - _i;
236 u$2[0] = s0 - (_i + bvirt) + (bvirt - t0);
237 _j = s1 + _i;
238 bvirt = _j - s1;
239 _0 = s1 - (_j - bvirt) + (_i - bvirt);
240 _i = _0 - t1;
241 bvirt = _0 - _i;
242 u$2[1] = _0 - (_i + bvirt) + (bvirt - t1);
243 u3 = _j + _i;
244 bvirt = u3 - _j;
245 u$2[2] = _j - (u3 - bvirt) + (_i - bvirt);
246 u$2[3] = u3;
247 const C1len = sum(4, B, 4, u$2, C1);
248
249 s1 = acx * bcytail;
250 c = splitter * acx;
251 ahi = c - (c - acx);
252 alo = acx - ahi;
253 c = splitter * bcytail;
254 bhi = c - (c - bcytail);
255 blo = bcytail - bhi;
256 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
257 t1 = acy * bcxtail;
258 c = splitter * acy;
259 ahi = c - (c - acy);
260 alo = acy - ahi;
261 c = splitter * bcxtail;
262 bhi = c - (c - bcxtail);
263 blo = bcxtail - bhi;
264 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
265 _i = s0 - t0;
266 bvirt = s0 - _i;
267 u$2[0] = s0 - (_i + bvirt) + (bvirt - t0);
268 _j = s1 + _i;
269 bvirt = _j - s1;
270 _0 = s1 - (_j - bvirt) + (_i - bvirt);
271 _i = _0 - t1;
272 bvirt = _0 - _i;
273 u$2[1] = _0 - (_i + bvirt) + (bvirt - t1);
274 u3 = _j + _i;
275 bvirt = u3 - _j;
276 u$2[2] = _j - (u3 - bvirt) + (_i - bvirt);
277 u$2[3] = u3;
278 const C2len = sum(C1len, C1, 4, u$2, C2);
279
280 s1 = acxtail * bcytail;
281 c = splitter * acxtail;
282 ahi = c - (c - acxtail);
283 alo = acxtail - ahi;
284 c = splitter * bcytail;
285 bhi = c - (c - bcytail);
286 blo = bcytail - bhi;
287 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
288 t1 = acytail * bcxtail;
289 c = splitter * acytail;
290 ahi = c - (c - acytail);
291 alo = acytail - ahi;
292 c = splitter * bcxtail;
293 bhi = c - (c - bcxtail);
294 blo = bcxtail - bhi;
295 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
296 _i = s0 - t0;
297 bvirt = s0 - _i;
298 u$2[0] = s0 - (_i + bvirt) + (bvirt - t0);
299 _j = s1 + _i;
300 bvirt = _j - s1;
301 _0 = s1 - (_j - bvirt) + (_i - bvirt);
302 _i = _0 - t1;
303 bvirt = _0 - _i;
304 u$2[1] = _0 - (_i + bvirt) + (bvirt - t1);
305 u3 = _j + _i;
306 bvirt = u3 - _j;
307 u$2[2] = _j - (u3 - bvirt) + (_i - bvirt);
308 u$2[3] = u3;
309 const Dlen = sum(C2len, C2, 4, u$2, D);
310
311 return D[Dlen - 1];
312}
313
314function orient2d(ax, ay, bx, by, cx, cy) {
315 const detleft = (ay - cy) * (bx - cx);
316 const detright = (ax - cx) * (by - cy);
317 const det = detleft - detright;
318
319 const detsum = Math.abs(detleft + detright);
320 if (Math.abs(det) >= ccwerrboundA * detsum) return det;
321
322 return -orient2dadapt(ax, ay, bx, by, cx, cy, detsum);
323}
324
325function orient2dfast(ax, ay, bx, by, cx, cy) {
326 return (ay - cy) * (bx - cx) - (ax - cx) * (by - cy);
327}
328
329const o3derrboundA = (7 + 56 * epsilon) * epsilon;
330const o3derrboundB = (3 + 28 * epsilon) * epsilon;
331const o3derrboundC = (26 + 288 * epsilon) * epsilon * epsilon;
332
333const bc$2 = vec(4);
334const ca$1 = vec(4);
335const ab$2 = vec(4);
336const at_b = vec(4);
337const at_c = vec(4);
338const bt_c = vec(4);
339const bt_a = vec(4);
340const ct_a = vec(4);
341const ct_b = vec(4);
342const bct$1 = vec(8);
343const cat$1 = vec(8);
344const abt$1 = vec(8);
345const u$1 = vec(4);
346
347const _8$2 = vec(8);
348const _8b$1 = vec(8);
349const _16$2 = vec(8);
350const _12 = vec(12);
351
352let fin$2 = vec(192);
353let fin2$1 = vec(192);
354
355function finadd$1(finlen, alen, a) {
356 finlen = sum(finlen, fin$2, alen, a, fin2$1);
357 const tmp = fin$2; fin$2 = fin2$1; fin2$1 = tmp;
358 return finlen;
359}
360
361function tailinit(xtail, ytail, ax, ay, bx, by, a, b) {
362 let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0, u3, negate;
363 if (xtail === 0) {
364 if (ytail === 0) {
365 a[0] = 0;
366 b[0] = 0;
367 return 1;
368 } else {
369 negate = -ytail;
370 s1 = negate * ax;
371 c = splitter * negate;
372 ahi = c - (c - negate);
373 alo = negate - ahi;
374 c = splitter * ax;
375 bhi = c - (c - ax);
376 blo = ax - bhi;
377 a[0] = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
378 a[1] = s1;
379 s1 = ytail * bx;
380 c = splitter * ytail;
381 ahi = c - (c - ytail);
382 alo = ytail - ahi;
383 c = splitter * bx;
384 bhi = c - (c - bx);
385 blo = bx - bhi;
386 b[0] = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
387 b[1] = s1;
388 return 2;
389 }
390 } else {
391 if (ytail === 0) {
392 s1 = xtail * ay;
393 c = splitter * xtail;
394 ahi = c - (c - xtail);
395 alo = xtail - ahi;
396 c = splitter * ay;
397 bhi = c - (c - ay);
398 blo = ay - bhi;
399 a[0] = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
400 a[1] = s1;
401 negate = -xtail;
402 s1 = negate * by;
403 c = splitter * negate;
404 ahi = c - (c - negate);
405 alo = negate - ahi;
406 c = splitter * by;
407 bhi = c - (c - by);
408 blo = by - bhi;
409 b[0] = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
410 b[1] = s1;
411 return 2;
412 } else {
413 s1 = xtail * ay;
414 c = splitter * xtail;
415 ahi = c - (c - xtail);
416 alo = xtail - ahi;
417 c = splitter * ay;
418 bhi = c - (c - ay);
419 blo = ay - bhi;
420 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
421 t1 = ytail * ax;
422 c = splitter * ytail;
423 ahi = c - (c - ytail);
424 alo = ytail - ahi;
425 c = splitter * ax;
426 bhi = c - (c - ax);
427 blo = ax - bhi;
428 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
429 _i = s0 - t0;
430 bvirt = s0 - _i;
431 a[0] = s0 - (_i + bvirt) + (bvirt - t0);
432 _j = s1 + _i;
433 bvirt = _j - s1;
434 _0 = s1 - (_j - bvirt) + (_i - bvirt);
435 _i = _0 - t1;
436 bvirt = _0 - _i;
437 a[1] = _0 - (_i + bvirt) + (bvirt - t1);
438 u3 = _j + _i;
439 bvirt = u3 - _j;
440 a[2] = _j - (u3 - bvirt) + (_i - bvirt);
441 a[3] = u3;
442 s1 = ytail * bx;
443 c = splitter * ytail;
444 ahi = c - (c - ytail);
445 alo = ytail - ahi;
446 c = splitter * bx;
447 bhi = c - (c - bx);
448 blo = bx - bhi;
449 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
450 t1 = xtail * by;
451 c = splitter * xtail;
452 ahi = c - (c - xtail);
453 alo = xtail - ahi;
454 c = splitter * by;
455 bhi = c - (c - by);
456 blo = by - bhi;
457 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
458 _i = s0 - t0;
459 bvirt = s0 - _i;
460 b[0] = s0 - (_i + bvirt) + (bvirt - t0);
461 _j = s1 + _i;
462 bvirt = _j - s1;
463 _0 = s1 - (_j - bvirt) + (_i - bvirt);
464 _i = _0 - t1;
465 bvirt = _0 - _i;
466 b[1] = _0 - (_i + bvirt) + (bvirt - t1);
467 u3 = _j + _i;
468 bvirt = u3 - _j;
469 b[2] = _j - (u3 - bvirt) + (_i - bvirt);
470 b[3] = u3;
471 return 4;
472 }
473 }
474}
475
476function tailadd(finlen, a, b, k, z) {
477 let bvirt, c, ahi, alo, bhi, blo, _i, _j, _k, _0, s1, s0, u3;
478 s1 = a * b;
479 c = splitter * a;
480 ahi = c - (c - a);
481 alo = a - ahi;
482 c = splitter * b;
483 bhi = c - (c - b);
484 blo = b - bhi;
485 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
486 c = splitter * k;
487 bhi = c - (c - k);
488 blo = k - bhi;
489 _i = s0 * k;
490 c = splitter * s0;
491 ahi = c - (c - s0);
492 alo = s0 - ahi;
493 u$1[0] = alo * blo - (_i - ahi * bhi - alo * bhi - ahi * blo);
494 _j = s1 * k;
495 c = splitter * s1;
496 ahi = c - (c - s1);
497 alo = s1 - ahi;
498 _0 = alo * blo - (_j - ahi * bhi - alo * bhi - ahi * blo);
499 _k = _i + _0;
500 bvirt = _k - _i;
501 u$1[1] = _i - (_k - bvirt) + (_0 - bvirt);
502 u3 = _j + _k;
503 u$1[2] = _k - (u3 - _j);
504 u$1[3] = u3;
505 finlen = finadd$1(finlen, 4, u$1);
506 if (z !== 0) {
507 c = splitter * z;
508 bhi = c - (c - z);
509 blo = z - bhi;
510 _i = s0 * z;
511 c = splitter * s0;
512 ahi = c - (c - s0);
513 alo = s0 - ahi;
514 u$1[0] = alo * blo - (_i - ahi * bhi - alo * bhi - ahi * blo);
515 _j = s1 * z;
516 c = splitter * s1;
517 ahi = c - (c - s1);
518 alo = s1 - ahi;
519 _0 = alo * blo - (_j - ahi * bhi - alo * bhi - ahi * blo);
520 _k = _i + _0;
521 bvirt = _k - _i;
522 u$1[1] = _i - (_k - bvirt) + (_0 - bvirt);
523 u3 = _j + _k;
524 u$1[2] = _k - (u3 - _j);
525 u$1[3] = u3;
526 finlen = finadd$1(finlen, 4, u$1);
527 }
528 return finlen;
529}
530
531function orient3dadapt(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, permanent) {
532 let finlen;
533 let adxtail, bdxtail, cdxtail;
534 let adytail, bdytail, cdytail;
535 let adztail, bdztail, cdztail;
536 let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0, u3;
537
538 const adx = ax - dx;
539 const bdx = bx - dx;
540 const cdx = cx - dx;
541 const ady = ay - dy;
542 const bdy = by - dy;
543 const cdy = cy - dy;
544 const adz = az - dz;
545 const bdz = bz - dz;
546 const cdz = cz - dz;
547
548 s1 = bdx * cdy;
549 c = splitter * bdx;
550 ahi = c - (c - bdx);
551 alo = bdx - ahi;
552 c = splitter * cdy;
553 bhi = c - (c - cdy);
554 blo = cdy - bhi;
555 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
556 t1 = cdx * bdy;
557 c = splitter * cdx;
558 ahi = c - (c - cdx);
559 alo = cdx - ahi;
560 c = splitter * bdy;
561 bhi = c - (c - bdy);
562 blo = bdy - bhi;
563 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
564 _i = s0 - t0;
565 bvirt = s0 - _i;
566 bc$2[0] = s0 - (_i + bvirt) + (bvirt - t0);
567 _j = s1 + _i;
568 bvirt = _j - s1;
569 _0 = s1 - (_j - bvirt) + (_i - bvirt);
570 _i = _0 - t1;
571 bvirt = _0 - _i;
572 bc$2[1] = _0 - (_i + bvirt) + (bvirt - t1);
573 u3 = _j + _i;
574 bvirt = u3 - _j;
575 bc$2[2] = _j - (u3 - bvirt) + (_i - bvirt);
576 bc$2[3] = u3;
577 s1 = cdx * ady;
578 c = splitter * cdx;
579 ahi = c - (c - cdx);
580 alo = cdx - ahi;
581 c = splitter * ady;
582 bhi = c - (c - ady);
583 blo = ady - bhi;
584 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
585 t1 = adx * cdy;
586 c = splitter * adx;
587 ahi = c - (c - adx);
588 alo = adx - ahi;
589 c = splitter * cdy;
590 bhi = c - (c - cdy);
591 blo = cdy - bhi;
592 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
593 _i = s0 - t0;
594 bvirt = s0 - _i;
595 ca$1[0] = s0 - (_i + bvirt) + (bvirt - t0);
596 _j = s1 + _i;
597 bvirt = _j - s1;
598 _0 = s1 - (_j - bvirt) + (_i - bvirt);
599 _i = _0 - t1;
600 bvirt = _0 - _i;
601 ca$1[1] = _0 - (_i + bvirt) + (bvirt - t1);
602 u3 = _j + _i;
603 bvirt = u3 - _j;
604 ca$1[2] = _j - (u3 - bvirt) + (_i - bvirt);
605 ca$1[3] = u3;
606 s1 = adx * bdy;
607 c = splitter * adx;
608 ahi = c - (c - adx);
609 alo = adx - ahi;
610 c = splitter * bdy;
611 bhi = c - (c - bdy);
612 blo = bdy - bhi;
613 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
614 t1 = bdx * ady;
615 c = splitter * bdx;
616 ahi = c - (c - bdx);
617 alo = bdx - ahi;
618 c = splitter * ady;
619 bhi = c - (c - ady);
620 blo = ady - bhi;
621 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
622 _i = s0 - t0;
623 bvirt = s0 - _i;
624 ab$2[0] = s0 - (_i + bvirt) + (bvirt - t0);
625 _j = s1 + _i;
626 bvirt = _j - s1;
627 _0 = s1 - (_j - bvirt) + (_i - bvirt);
628 _i = _0 - t1;
629 bvirt = _0 - _i;
630 ab$2[1] = _0 - (_i + bvirt) + (bvirt - t1);
631 u3 = _j + _i;
632 bvirt = u3 - _j;
633 ab$2[2] = _j - (u3 - bvirt) + (_i - bvirt);
634 ab$2[3] = u3;
635
636 finlen = sum(
637 sum(
638 scale(4, bc$2, adz, _8$2), _8$2,
639 scale(4, ca$1, bdz, _8b$1), _8b$1, _16$2), _16$2,
640 scale(4, ab$2, cdz, _8$2), _8$2, fin$2);
641
642 let det = estimate(finlen, fin$2);
643 let errbound = o3derrboundB * permanent;
644 if (det >= errbound || -det >= errbound) {
645 return det;
646 }
647
648 bvirt = ax - adx;
649 adxtail = ax - (adx + bvirt) + (bvirt - dx);
650 bvirt = bx - bdx;
651 bdxtail = bx - (bdx + bvirt) + (bvirt - dx);
652 bvirt = cx - cdx;
653 cdxtail = cx - (cdx + bvirt) + (bvirt - dx);
654 bvirt = ay - ady;
655 adytail = ay - (ady + bvirt) + (bvirt - dy);
656 bvirt = by - bdy;
657 bdytail = by - (bdy + bvirt) + (bvirt - dy);
658 bvirt = cy - cdy;
659 cdytail = cy - (cdy + bvirt) + (bvirt - dy);
660 bvirt = az - adz;
661 adztail = az - (adz + bvirt) + (bvirt - dz);
662 bvirt = bz - bdz;
663 bdztail = bz - (bdz + bvirt) + (bvirt - dz);
664 bvirt = cz - cdz;
665 cdztail = cz - (cdz + bvirt) + (bvirt - dz);
666
667 if (adxtail === 0 && bdxtail === 0 && cdxtail === 0 &&
668 adytail === 0 && bdytail === 0 && cdytail === 0 &&
669 adztail === 0 && bdztail === 0 && cdztail === 0) {
670 return det;
671 }
672
673 errbound = o3derrboundC * permanent + resulterrbound * Math.abs(det);
674 det +=
675 adz * (bdx * cdytail + cdy * bdxtail - (bdy * cdxtail + cdx * bdytail)) + adztail * (bdx * cdy - bdy * cdx) +
676 bdz * (cdx * adytail + ady * cdxtail - (cdy * adxtail + adx * cdytail)) + bdztail * (cdx * ady - cdy * adx) +
677 cdz * (adx * bdytail + bdy * adxtail - (ady * bdxtail + bdx * adytail)) + cdztail * (adx * bdy - ady * bdx);
678 if (det >= errbound || -det >= errbound) {
679 return det;
680 }
681
682 const at_len = tailinit(adxtail, adytail, bdx, bdy, cdx, cdy, at_b, at_c);
683 const bt_len = tailinit(bdxtail, bdytail, cdx, cdy, adx, ady, bt_c, bt_a);
684 const ct_len = tailinit(cdxtail, cdytail, adx, ady, bdx, bdy, ct_a, ct_b);
685
686 const bctlen = sum(bt_len, bt_c, ct_len, ct_b, bct$1);
687 finlen = finadd$1(finlen, scale(bctlen, bct$1, adz, _16$2), _16$2);
688
689 const catlen = sum(ct_len, ct_a, at_len, at_c, cat$1);
690 finlen = finadd$1(finlen, scale(catlen, cat$1, bdz, _16$2), _16$2);
691
692 const abtlen = sum(at_len, at_b, bt_len, bt_a, abt$1);
693 finlen = finadd$1(finlen, scale(abtlen, abt$1, cdz, _16$2), _16$2);
694
695 if (adztail !== 0) {
696 finlen = finadd$1(finlen, scale(4, bc$2, adztail, _12), _12);
697 finlen = finadd$1(finlen, scale(bctlen, bct$1, adztail, _16$2), _16$2);
698 }
699 if (bdztail !== 0) {
700 finlen = finadd$1(finlen, scale(4, ca$1, bdztail, _12), _12);
701 finlen = finadd$1(finlen, scale(catlen, cat$1, bdztail, _16$2), _16$2);
702 }
703 if (cdztail !== 0) {
704 finlen = finadd$1(finlen, scale(4, ab$2, cdztail, _12), _12);
705 finlen = finadd$1(finlen, scale(abtlen, abt$1, cdztail, _16$2), _16$2);
706 }
707
708 if (adxtail !== 0) {
709 if (bdytail !== 0) {
710 finlen = tailadd(finlen, adxtail, bdytail, cdz, cdztail);
711 }
712 if (cdytail !== 0) {
713 finlen = tailadd(finlen, -adxtail, cdytail, bdz, bdztail);
714 }
715 }
716 if (bdxtail !== 0) {
717 if (cdytail !== 0) {
718 finlen = tailadd(finlen, bdxtail, cdytail, adz, adztail);
719 }
720 if (adytail !== 0) {
721 finlen = tailadd(finlen, -bdxtail, adytail, cdz, cdztail);
722 }
723 }
724 if (cdxtail !== 0) {
725 if (adytail !== 0) {
726 finlen = tailadd(finlen, cdxtail, adytail, bdz, bdztail);
727 }
728 if (bdytail !== 0) {
729 finlen = tailadd(finlen, -cdxtail, bdytail, adz, adztail);
730 }
731 }
732
733 return fin$2[finlen - 1];
734}
735
736function orient3d(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz) {
737 const adx = ax - dx;
738 const bdx = bx - dx;
739 const cdx = cx - dx;
740 const ady = ay - dy;
741 const bdy = by - dy;
742 const cdy = cy - dy;
743 const adz = az - dz;
744 const bdz = bz - dz;
745 const cdz = cz - dz;
746
747 const bdxcdy = bdx * cdy;
748 const cdxbdy = cdx * bdy;
749
750 const cdxady = cdx * ady;
751 const adxcdy = adx * cdy;
752
753 const adxbdy = adx * bdy;
754 const bdxady = bdx * ady;
755
756 const det =
757 adz * (bdxcdy - cdxbdy) +
758 bdz * (cdxady - adxcdy) +
759 cdz * (adxbdy - bdxady);
760
761 const permanent =
762 (Math.abs(bdxcdy) + Math.abs(cdxbdy)) * Math.abs(adz) +
763 (Math.abs(cdxady) + Math.abs(adxcdy)) * Math.abs(bdz) +
764 (Math.abs(adxbdy) + Math.abs(bdxady)) * Math.abs(cdz);
765
766 const errbound = o3derrboundA * permanent;
767 if (det > errbound || -det > errbound) {
768 return det;
769 }
770
771 return orient3dadapt(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, permanent);
772}
773
774function orient3dfast(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz) {
775 const adx = ax - dx;
776 const bdx = bx - dx;
777 const cdx = cx - dx;
778 const ady = ay - dy;
779 const bdy = by - dy;
780 const cdy = cy - dy;
781 const adz = az - dz;
782 const bdz = bz - dz;
783 const cdz = cz - dz;
784
785 return adx * (bdy * cdz - bdz * cdy) +
786 bdx * (cdy * adz - cdz * ady) +
787 cdx * (ady * bdz - adz * bdy);
788}
789
790const iccerrboundA = (10 + 96 * epsilon) * epsilon;
791const iccerrboundB = (4 + 48 * epsilon) * epsilon;
792const iccerrboundC = (44 + 576 * epsilon) * epsilon * epsilon;
793
794const bc$1 = vec(4);
795const ca = vec(4);
796const ab$1 = vec(4);
797const aa = vec(4);
798const bb = vec(4);
799const cc = vec(4);
800const u = vec(4);
801const v = vec(4);
802const axtbc = vec(8);
803const aytbc = vec(8);
804const bxtca = vec(8);
805const bytca = vec(8);
806const cxtab = vec(8);
807const cytab = vec(8);
808const abt = vec(8);
809const bct = vec(8);
810const cat = vec(8);
811const abtt = vec(4);
812const bctt = vec(4);
813const catt = vec(4);
814
815const _8$1 = vec(8);
816const _16$1 = vec(16);
817const _16b = vec(16);
818const _16c = vec(16);
819const _32 = vec(32);
820const _32b = vec(32);
821const _48$1 = vec(48);
822const _64 = vec(64);
823
824let fin$1 = vec(1152);
825let fin2 = vec(1152);
826
827function finadd(finlen, a, alen) {
828 finlen = sum(finlen, fin$1, a, alen, fin2);
829 const tmp = fin$1; fin$1 = fin2; fin2 = tmp;
830 return finlen;
831}
832
833function incircleadapt(ax, ay, bx, by, cx, cy, dx, dy, permanent) {
834 let finlen;
835 let adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
836 let axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
837 let abtlen, bctlen, catlen;
838 let abttlen, bcttlen, cattlen;
839 let n1, n0;
840
841 let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0, u3;
842
843 const adx = ax - dx;
844 const bdx = bx - dx;
845 const cdx = cx - dx;
846 const ady = ay - dy;
847 const bdy = by - dy;
848 const cdy = cy - dy;
849
850 s1 = bdx * cdy;
851 c = splitter * bdx;
852 ahi = c - (c - bdx);
853 alo = bdx - ahi;
854 c = splitter * cdy;
855 bhi = c - (c - cdy);
856 blo = cdy - bhi;
857 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
858 t1 = cdx * bdy;
859 c = splitter * cdx;
860 ahi = c - (c - cdx);
861 alo = cdx - ahi;
862 c = splitter * bdy;
863 bhi = c - (c - bdy);
864 blo = bdy - bhi;
865 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
866 _i = s0 - t0;
867 bvirt = s0 - _i;
868 bc$1[0] = s0 - (_i + bvirt) + (bvirt - t0);
869 _j = s1 + _i;
870 bvirt = _j - s1;
871 _0 = s1 - (_j - bvirt) + (_i - bvirt);
872 _i = _0 - t1;
873 bvirt = _0 - _i;
874 bc$1[1] = _0 - (_i + bvirt) + (bvirt - t1);
875 u3 = _j + _i;
876 bvirt = u3 - _j;
877 bc$1[2] = _j - (u3 - bvirt) + (_i - bvirt);
878 bc$1[3] = u3;
879 s1 = cdx * ady;
880 c = splitter * cdx;
881 ahi = c - (c - cdx);
882 alo = cdx - ahi;
883 c = splitter * ady;
884 bhi = c - (c - ady);
885 blo = ady - bhi;
886 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
887 t1 = adx * cdy;
888 c = splitter * adx;
889 ahi = c - (c - adx);
890 alo = adx - ahi;
891 c = splitter * cdy;
892 bhi = c - (c - cdy);
893 blo = cdy - bhi;
894 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
895 _i = s0 - t0;
896 bvirt = s0 - _i;
897 ca[0] = s0 - (_i + bvirt) + (bvirt - t0);
898 _j = s1 + _i;
899 bvirt = _j - s1;
900 _0 = s1 - (_j - bvirt) + (_i - bvirt);
901 _i = _0 - t1;
902 bvirt = _0 - _i;
903 ca[1] = _0 - (_i + bvirt) + (bvirt - t1);
904 u3 = _j + _i;
905 bvirt = u3 - _j;
906 ca[2] = _j - (u3 - bvirt) + (_i - bvirt);
907 ca[3] = u3;
908 s1 = adx * bdy;
909 c = splitter * adx;
910 ahi = c - (c - adx);
911 alo = adx - ahi;
912 c = splitter * bdy;
913 bhi = c - (c - bdy);
914 blo = bdy - bhi;
915 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
916 t1 = bdx * ady;
917 c = splitter * bdx;
918 ahi = c - (c - bdx);
919 alo = bdx - ahi;
920 c = splitter * ady;
921 bhi = c - (c - ady);
922 blo = ady - bhi;
923 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
924 _i = s0 - t0;
925 bvirt = s0 - _i;
926 ab$1[0] = s0 - (_i + bvirt) + (bvirt - t0);
927 _j = s1 + _i;
928 bvirt = _j - s1;
929 _0 = s1 - (_j - bvirt) + (_i - bvirt);
930 _i = _0 - t1;
931 bvirt = _0 - _i;
932 ab$1[1] = _0 - (_i + bvirt) + (bvirt - t1);
933 u3 = _j + _i;
934 bvirt = u3 - _j;
935 ab$1[2] = _j - (u3 - bvirt) + (_i - bvirt);
936 ab$1[3] = u3;
937
938 finlen = sum(
939 sum(
940 sum(
941 scale(scale(4, bc$1, adx, _8$1), _8$1, adx, _16$1), _16$1,
942 scale(scale(4, bc$1, ady, _8$1), _8$1, ady, _16b), _16b, _32), _32,
943 sum(
944 scale(scale(4, ca, bdx, _8$1), _8$1, bdx, _16$1), _16$1,
945 scale(scale(4, ca, bdy, _8$1), _8$1, bdy, _16b), _16b, _32b), _32b, _64), _64,
946 sum(
947 scale(scale(4, ab$1, cdx, _8$1), _8$1, cdx, _16$1), _16$1,
948 scale(scale(4, ab$1, cdy, _8$1), _8$1, cdy, _16b), _16b, _32), _32, fin$1);
949
950 let det = estimate(finlen, fin$1);
951 let errbound = iccerrboundB * permanent;
952 if (det >= errbound || -det >= errbound) {
953 return det;
954 }
955
956 bvirt = ax - adx;
957 adxtail = ax - (adx + bvirt) + (bvirt - dx);
958 bvirt = ay - ady;
959 adytail = ay - (ady + bvirt) + (bvirt - dy);
960 bvirt = bx - bdx;
961 bdxtail = bx - (bdx + bvirt) + (bvirt - dx);
962 bvirt = by - bdy;
963 bdytail = by - (bdy + bvirt) + (bvirt - dy);
964 bvirt = cx - cdx;
965 cdxtail = cx - (cdx + bvirt) + (bvirt - dx);
966 bvirt = cy - cdy;
967 cdytail = cy - (cdy + bvirt) + (bvirt - dy);
968 if (adxtail === 0 && bdxtail === 0 && cdxtail === 0 && adytail === 0 && bdytail === 0 && cdytail === 0) {
969 return det;
970 }
971
972 errbound = iccerrboundC * permanent + resulterrbound * Math.abs(det);
973 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail)) +
974 2 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) +
975 ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail)) +
976 2 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) +
977 ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail)) +
978 2 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
979
980 if (det >= errbound || -det >= errbound) {
981 return det;
982 }
983
984 if (bdxtail !== 0 || bdytail !== 0 || cdxtail !== 0 || cdytail !== 0) {
985 s1 = adx * adx;
986 c = splitter * adx;
987 ahi = c - (c - adx);
988 alo = adx - ahi;
989 s0 = alo * alo - (s1 - ahi * ahi - (ahi + ahi) * alo);
990 t1 = ady * ady;
991 c = splitter * ady;
992 ahi = c - (c - ady);
993 alo = ady - ahi;
994 t0 = alo * alo - (t1 - ahi * ahi - (ahi + ahi) * alo);
995 _i = s0 + t0;
996 bvirt = _i - s0;
997 aa[0] = s0 - (_i - bvirt) + (t0 - bvirt);
998 _j = s1 + _i;
999 bvirt = _j - s1;
1000 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1001 _i = _0 + t1;
1002 bvirt = _i - _0;
1003 aa[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1004 u3 = _j + _i;
1005 bvirt = u3 - _j;
1006 aa[2] = _j - (u3 - bvirt) + (_i - bvirt);
1007 aa[3] = u3;
1008 }
1009 if (cdxtail !== 0 || cdytail !== 0 || adxtail !== 0 || adytail !== 0) {
1010 s1 = bdx * bdx;
1011 c = splitter * bdx;
1012 ahi = c - (c - bdx);
1013 alo = bdx - ahi;
1014 s0 = alo * alo - (s1 - ahi * ahi - (ahi + ahi) * alo);
1015 t1 = bdy * bdy;
1016 c = splitter * bdy;
1017 ahi = c - (c - bdy);
1018 alo = bdy - ahi;
1019 t0 = alo * alo - (t1 - ahi * ahi - (ahi + ahi) * alo);
1020 _i = s0 + t0;
1021 bvirt = _i - s0;
1022 bb[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1023 _j = s1 + _i;
1024 bvirt = _j - s1;
1025 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1026 _i = _0 + t1;
1027 bvirt = _i - _0;
1028 bb[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1029 u3 = _j + _i;
1030 bvirt = u3 - _j;
1031 bb[2] = _j - (u3 - bvirt) + (_i - bvirt);
1032 bb[3] = u3;
1033 }
1034 if (adxtail !== 0 || adytail !== 0 || bdxtail !== 0 || bdytail !== 0) {
1035 s1 = cdx * cdx;
1036 c = splitter * cdx;
1037 ahi = c - (c - cdx);
1038 alo = cdx - ahi;
1039 s0 = alo * alo - (s1 - ahi * ahi - (ahi + ahi) * alo);
1040 t1 = cdy * cdy;
1041 c = splitter * cdy;
1042 ahi = c - (c - cdy);
1043 alo = cdy - ahi;
1044 t0 = alo * alo - (t1 - ahi * ahi - (ahi + ahi) * alo);
1045 _i = s0 + t0;
1046 bvirt = _i - s0;
1047 cc[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1048 _j = s1 + _i;
1049 bvirt = _j - s1;
1050 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1051 _i = _0 + t1;
1052 bvirt = _i - _0;
1053 cc[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1054 u3 = _j + _i;
1055 bvirt = u3 - _j;
1056 cc[2] = _j - (u3 - bvirt) + (_i - bvirt);
1057 cc[3] = u3;
1058 }
1059
1060 if (adxtail !== 0) {
1061 axtbclen = scale(4, bc$1, adxtail, axtbc);
1062 finlen = finadd(finlen, sum_three(
1063 scale(axtbclen, axtbc, 2 * adx, _16$1), _16$1,
1064 scale(scale(4, cc, adxtail, _8$1), _8$1, bdy, _16b), _16b,
1065 scale(scale(4, bb, adxtail, _8$1), _8$1, -cdy, _16c), _16c, _32, _48$1), _48$1);
1066 }
1067 if (adytail !== 0) {
1068 aytbclen = scale(4, bc$1, adytail, aytbc);
1069 finlen = finadd(finlen, sum_three(
1070 scale(aytbclen, aytbc, 2 * ady, _16$1), _16$1,
1071 scale(scale(4, bb, adytail, _8$1), _8$1, cdx, _16b), _16b,
1072 scale(scale(4, cc, adytail, _8$1), _8$1, -bdx, _16c), _16c, _32, _48$1), _48$1);
1073 }
1074 if (bdxtail !== 0) {
1075 bxtcalen = scale(4, ca, bdxtail, bxtca);
1076 finlen = finadd(finlen, sum_three(
1077 scale(bxtcalen, bxtca, 2 * bdx, _16$1), _16$1,
1078 scale(scale(4, aa, bdxtail, _8$1), _8$1, cdy, _16b), _16b,
1079 scale(scale(4, cc, bdxtail, _8$1), _8$1, -ady, _16c), _16c, _32, _48$1), _48$1);
1080 }
1081 if (bdytail !== 0) {
1082 bytcalen = scale(4, ca, bdytail, bytca);
1083 finlen = finadd(finlen, sum_three(
1084 scale(bytcalen, bytca, 2 * bdy, _16$1), _16$1,
1085 scale(scale(4, cc, bdytail, _8$1), _8$1, adx, _16b), _16b,
1086 scale(scale(4, aa, bdytail, _8$1), _8$1, -cdx, _16c), _16c, _32, _48$1), _48$1);
1087 }
1088 if (cdxtail !== 0) {
1089 cxtablen = scale(4, ab$1, cdxtail, cxtab);
1090 finlen = finadd(finlen, sum_three(
1091 scale(cxtablen, cxtab, 2 * cdx, _16$1), _16$1,
1092 scale(scale(4, bb, cdxtail, _8$1), _8$1, ady, _16b), _16b,
1093 scale(scale(4, aa, cdxtail, _8$1), _8$1, -bdy, _16c), _16c, _32, _48$1), _48$1);
1094 }
1095 if (cdytail !== 0) {
1096 cytablen = scale(4, ab$1, cdytail, cytab);
1097 finlen = finadd(finlen, sum_three(
1098 scale(cytablen, cytab, 2 * cdy, _16$1), _16$1,
1099 scale(scale(4, aa, cdytail, _8$1), _8$1, bdx, _16b), _16b,
1100 scale(scale(4, bb, cdytail, _8$1), _8$1, -adx, _16c), _16c, _32, _48$1), _48$1);
1101 }
1102
1103 if (adxtail !== 0 || adytail !== 0) {
1104 if (bdxtail !== 0 || bdytail !== 0 || cdxtail !== 0 || cdytail !== 0) {
1105 s1 = bdxtail * cdy;
1106 c = splitter * bdxtail;
1107 ahi = c - (c - bdxtail);
1108 alo = bdxtail - ahi;
1109 c = splitter * cdy;
1110 bhi = c - (c - cdy);
1111 blo = cdy - bhi;
1112 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1113 t1 = bdx * cdytail;
1114 c = splitter * bdx;
1115 ahi = c - (c - bdx);
1116 alo = bdx - ahi;
1117 c = splitter * cdytail;
1118 bhi = c - (c - cdytail);
1119 blo = cdytail - bhi;
1120 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1121 _i = s0 + t0;
1122 bvirt = _i - s0;
1123 u[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1124 _j = s1 + _i;
1125 bvirt = _j - s1;
1126 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1127 _i = _0 + t1;
1128 bvirt = _i - _0;
1129 u[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1130 u3 = _j + _i;
1131 bvirt = u3 - _j;
1132 u[2] = _j - (u3 - bvirt) + (_i - bvirt);
1133 u[3] = u3;
1134 s1 = cdxtail * -bdy;
1135 c = splitter * cdxtail;
1136 ahi = c - (c - cdxtail);
1137 alo = cdxtail - ahi;
1138 c = splitter * -bdy;
1139 bhi = c - (c - -bdy);
1140 blo = -bdy - bhi;
1141 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1142 t1 = cdx * -bdytail;
1143 c = splitter * cdx;
1144 ahi = c - (c - cdx);
1145 alo = cdx - ahi;
1146 c = splitter * -bdytail;
1147 bhi = c - (c - -bdytail);
1148 blo = -bdytail - bhi;
1149 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1150 _i = s0 + t0;
1151 bvirt = _i - s0;
1152 v[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1153 _j = s1 + _i;
1154 bvirt = _j - s1;
1155 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1156 _i = _0 + t1;
1157 bvirt = _i - _0;
1158 v[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1159 u3 = _j + _i;
1160 bvirt = u3 - _j;
1161 v[2] = _j - (u3 - bvirt) + (_i - bvirt);
1162 v[3] = u3;
1163 bctlen = sum(4, u, 4, v, bct);
1164 s1 = bdxtail * cdytail;
1165 c = splitter * bdxtail;
1166 ahi = c - (c - bdxtail);
1167 alo = bdxtail - ahi;
1168 c = splitter * cdytail;
1169 bhi = c - (c - cdytail);
1170 blo = cdytail - bhi;
1171 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1172 t1 = cdxtail * bdytail;
1173 c = splitter * cdxtail;
1174 ahi = c - (c - cdxtail);
1175 alo = cdxtail - ahi;
1176 c = splitter * bdytail;
1177 bhi = c - (c - bdytail);
1178 blo = bdytail - bhi;
1179 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1180 _i = s0 - t0;
1181 bvirt = s0 - _i;
1182 bctt[0] = s0 - (_i + bvirt) + (bvirt - t0);
1183 _j = s1 + _i;
1184 bvirt = _j - s1;
1185 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1186 _i = _0 - t1;
1187 bvirt = _0 - _i;
1188 bctt[1] = _0 - (_i + bvirt) + (bvirt - t1);
1189 u3 = _j + _i;
1190 bvirt = u3 - _j;
1191 bctt[2] = _j - (u3 - bvirt) + (_i - bvirt);
1192 bctt[3] = u3;
1193 bcttlen = 4;
1194 } else {
1195 bct[0] = 0;
1196 bctlen = 1;
1197 bctt[0] = 0;
1198 bcttlen = 1;
1199 }
1200 if (adxtail !== 0) {
1201 const len = scale(bctlen, bct, adxtail, _16c);
1202 finlen = finadd(finlen, sum(
1203 scale(axtbclen, axtbc, adxtail, _16$1), _16$1,
1204 scale(len, _16c, 2 * adx, _32), _32, _48$1), _48$1);
1205
1206 const len2 = scale(bcttlen, bctt, adxtail, _8$1);
1207 finlen = finadd(finlen, sum_three(
1208 scale(len2, _8$1, 2 * adx, _16$1), _16$1,
1209 scale(len2, _8$1, adxtail, _16b), _16b,
1210 scale(len, _16c, adxtail, _32), _32, _32b, _64), _64);
1211
1212 if (bdytail !== 0) {
1213 finlen = finadd(finlen, scale(scale(4, cc, adxtail, _8$1), _8$1, bdytail, _16$1), _16$1);
1214 }
1215 if (cdytail !== 0) {
1216 finlen = finadd(finlen, scale(scale(4, bb, -adxtail, _8$1), _8$1, cdytail, _16$1), _16$1);
1217 }
1218 }
1219 if (adytail !== 0) {
1220 const len = scale(bctlen, bct, adytail, _16c);
1221 finlen = finadd(finlen, sum(
1222 scale(aytbclen, aytbc, adytail, _16$1), _16$1,
1223 scale(len, _16c, 2 * ady, _32), _32, _48$1), _48$1);
1224
1225 const len2 = scale(bcttlen, bctt, adytail, _8$1);
1226 finlen = finadd(finlen, sum_three(
1227 scale(len2, _8$1, 2 * ady, _16$1), _16$1,
1228 scale(len2, _8$1, adytail, _16b), _16b,
1229 scale(len, _16c, adytail, _32), _32, _32b, _64), _64);
1230 }
1231 }
1232 if (bdxtail !== 0 || bdytail !== 0) {
1233 if (cdxtail !== 0 || cdytail !== 0 || adxtail !== 0 || adytail !== 0) {
1234 s1 = cdxtail * ady;
1235 c = splitter * cdxtail;
1236 ahi = c - (c - cdxtail);
1237 alo = cdxtail - ahi;
1238 c = splitter * ady;
1239 bhi = c - (c - ady);
1240 blo = ady - bhi;
1241 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1242 t1 = cdx * adytail;
1243 c = splitter * cdx;
1244 ahi = c - (c - cdx);
1245 alo = cdx - ahi;
1246 c = splitter * adytail;
1247 bhi = c - (c - adytail);
1248 blo = adytail - bhi;
1249 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1250 _i = s0 + t0;
1251 bvirt = _i - s0;
1252 u[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1253 _j = s1 + _i;
1254 bvirt = _j - s1;
1255 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1256 _i = _0 + t1;
1257 bvirt = _i - _0;
1258 u[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1259 u3 = _j + _i;
1260 bvirt = u3 - _j;
1261 u[2] = _j - (u3 - bvirt) + (_i - bvirt);
1262 u[3] = u3;
1263 n1 = -cdy;
1264 n0 = -cdytail;
1265 s1 = adxtail * n1;
1266 c = splitter * adxtail;
1267 ahi = c - (c - adxtail);
1268 alo = adxtail - ahi;
1269 c = splitter * n1;
1270 bhi = c - (c - n1);
1271 blo = n1 - bhi;
1272 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1273 t1 = adx * n0;
1274 c = splitter * adx;
1275 ahi = c - (c - adx);
1276 alo = adx - ahi;
1277 c = splitter * n0;
1278 bhi = c - (c - n0);
1279 blo = n0 - bhi;
1280 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1281 _i = s0 + t0;
1282 bvirt = _i - s0;
1283 v[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1284 _j = s1 + _i;
1285 bvirt = _j - s1;
1286 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1287 _i = _0 + t1;
1288 bvirt = _i - _0;
1289 v[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1290 u3 = _j + _i;
1291 bvirt = u3 - _j;
1292 v[2] = _j - (u3 - bvirt) + (_i - bvirt);
1293 v[3] = u3;
1294 catlen = sum(4, u, 4, v, cat);
1295 s1 = cdxtail * adytail;
1296 c = splitter * cdxtail;
1297 ahi = c - (c - cdxtail);
1298 alo = cdxtail - ahi;
1299 c = splitter * adytail;
1300 bhi = c - (c - adytail);
1301 blo = adytail - bhi;
1302 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1303 t1 = adxtail * cdytail;
1304 c = splitter * adxtail;
1305 ahi = c - (c - adxtail);
1306 alo = adxtail - ahi;
1307 c = splitter * cdytail;
1308 bhi = c - (c - cdytail);
1309 blo = cdytail - bhi;
1310 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1311 _i = s0 - t0;
1312 bvirt = s0 - _i;
1313 catt[0] = s0 - (_i + bvirt) + (bvirt - t0);
1314 _j = s1 + _i;
1315 bvirt = _j - s1;
1316 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1317 _i = _0 - t1;
1318 bvirt = _0 - _i;
1319 catt[1] = _0 - (_i + bvirt) + (bvirt - t1);
1320 u3 = _j + _i;
1321 bvirt = u3 - _j;
1322 catt[2] = _j - (u3 - bvirt) + (_i - bvirt);
1323 catt[3] = u3;
1324 cattlen = 4;
1325 } else {
1326 cat[0] = 0;
1327 catlen = 1;
1328 catt[0] = 0;
1329 cattlen = 1;
1330 }
1331 if (bdxtail !== 0) {
1332 const len = scale(catlen, cat, bdxtail, _16c);
1333 finlen = finadd(finlen, sum(
1334 scale(bxtcalen, bxtca, bdxtail, _16$1), _16$1,
1335 scale(len, _16c, 2 * bdx, _32), _32, _48$1), _48$1);
1336
1337 const len2 = scale(cattlen, catt, bdxtail, _8$1);
1338 finlen = finadd(finlen, sum_three(
1339 scale(len2, _8$1, 2 * bdx, _16$1), _16$1,
1340 scale(len2, _8$1, bdxtail, _16b), _16b,
1341 scale(len, _16c, bdxtail, _32), _32, _32b, _64), _64);
1342
1343 if (cdytail !== 0) {
1344 finlen = finadd(finlen, scale(scale(4, aa, bdxtail, _8$1), _8$1, cdytail, _16$1), _16$1);
1345 }
1346 if (adytail !== 0) {
1347 finlen = finadd(finlen, scale(scale(4, cc, -bdxtail, _8$1), _8$1, adytail, _16$1), _16$1);
1348 }
1349 }
1350 if (bdytail !== 0) {
1351 const len = scale(catlen, cat, bdytail, _16c);
1352 finlen = finadd(finlen, sum(
1353 scale(bytcalen, bytca, bdytail, _16$1), _16$1,
1354 scale(len, _16c, 2 * bdy, _32), _32, _48$1), _48$1);
1355
1356 const len2 = scale(cattlen, catt, bdytail, _8$1);
1357 finlen = finadd(finlen, sum_three(
1358 scale(len2, _8$1, 2 * bdy, _16$1), _16$1,
1359 scale(len2, _8$1, bdytail, _16b), _16b,
1360 scale(len, _16c, bdytail, _32), _32, _32b, _64), _64);
1361 }
1362 }
1363 if (cdxtail !== 0 || cdytail !== 0) {
1364 if (adxtail !== 0 || adytail !== 0 || bdxtail !== 0 || bdytail !== 0) {
1365 s1 = adxtail * bdy;
1366 c = splitter * adxtail;
1367 ahi = c - (c - adxtail);
1368 alo = adxtail - ahi;
1369 c = splitter * bdy;
1370 bhi = c - (c - bdy);
1371 blo = bdy - bhi;
1372 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1373 t1 = adx * bdytail;
1374 c = splitter * adx;
1375 ahi = c - (c - adx);
1376 alo = adx - ahi;
1377 c = splitter * bdytail;
1378 bhi = c - (c - bdytail);
1379 blo = bdytail - bhi;
1380 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1381 _i = s0 + t0;
1382 bvirt = _i - s0;
1383 u[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1384 _j = s1 + _i;
1385 bvirt = _j - s1;
1386 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1387 _i = _0 + t1;
1388 bvirt = _i - _0;
1389 u[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1390 u3 = _j + _i;
1391 bvirt = u3 - _j;
1392 u[2] = _j - (u3 - bvirt) + (_i - bvirt);
1393 u[3] = u3;
1394 n1 = -ady;
1395 n0 = -adytail;
1396 s1 = bdxtail * n1;
1397 c = splitter * bdxtail;
1398 ahi = c - (c - bdxtail);
1399 alo = bdxtail - ahi;
1400 c = splitter * n1;
1401 bhi = c - (c - n1);
1402 blo = n1 - bhi;
1403 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1404 t1 = bdx * n0;
1405 c = splitter * bdx;
1406 ahi = c - (c - bdx);
1407 alo = bdx - ahi;
1408 c = splitter * n0;
1409 bhi = c - (c - n0);
1410 blo = n0 - bhi;
1411 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1412 _i = s0 + t0;
1413 bvirt = _i - s0;
1414 v[0] = s0 - (_i - bvirt) + (t0 - bvirt);
1415 _j = s1 + _i;
1416 bvirt = _j - s1;
1417 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1418 _i = _0 + t1;
1419 bvirt = _i - _0;
1420 v[1] = _0 - (_i - bvirt) + (t1 - bvirt);
1421 u3 = _j + _i;
1422 bvirt = u3 - _j;
1423 v[2] = _j - (u3 - bvirt) + (_i - bvirt);
1424 v[3] = u3;
1425 abtlen = sum(4, u, 4, v, abt);
1426 s1 = adxtail * bdytail;
1427 c = splitter * adxtail;
1428 ahi = c - (c - adxtail);
1429 alo = adxtail - ahi;
1430 c = splitter * bdytail;
1431 bhi = c - (c - bdytail);
1432 blo = bdytail - bhi;
1433 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1434 t1 = bdxtail * adytail;
1435 c = splitter * bdxtail;
1436 ahi = c - (c - bdxtail);
1437 alo = bdxtail - ahi;
1438 c = splitter * adytail;
1439 bhi = c - (c - adytail);
1440 blo = adytail - bhi;
1441 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1442 _i = s0 - t0;
1443 bvirt = s0 - _i;
1444 abtt[0] = s0 - (_i + bvirt) + (bvirt - t0);
1445 _j = s1 + _i;
1446 bvirt = _j - s1;
1447 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1448 _i = _0 - t1;
1449 bvirt = _0 - _i;
1450 abtt[1] = _0 - (_i + bvirt) + (bvirt - t1);
1451 u3 = _j + _i;
1452 bvirt = u3 - _j;
1453 abtt[2] = _j - (u3 - bvirt) + (_i - bvirt);
1454 abtt[3] = u3;
1455 abttlen = 4;
1456 } else {
1457 abt[0] = 0;
1458 abtlen = 1;
1459 abtt[0] = 0;
1460 abttlen = 1;
1461 }
1462 if (cdxtail !== 0) {
1463 const len = scale(abtlen, abt, cdxtail, _16c);
1464 finlen = finadd(finlen, sum(
1465 scale(cxtablen, cxtab, cdxtail, _16$1), _16$1,
1466 scale(len, _16c, 2 * cdx, _32), _32, _48$1), _48$1);
1467
1468 const len2 = scale(abttlen, abtt, cdxtail, _8$1);
1469 finlen = finadd(finlen, sum_three(
1470 scale(len2, _8$1, 2 * cdx, _16$1), _16$1,
1471 scale(len2, _8$1, cdxtail, _16b), _16b,
1472 scale(len, _16c, cdxtail, _32), _32, _32b, _64), _64);
1473
1474 if (adytail !== 0) {
1475 finlen = finadd(finlen, scale(scale(4, bb, cdxtail, _8$1), _8$1, adytail, _16$1), _16$1);
1476 }
1477 if (bdytail !== 0) {
1478 finlen = finadd(finlen, scale(scale(4, aa, -cdxtail, _8$1), _8$1, bdytail, _16$1), _16$1);
1479 }
1480 }
1481 if (cdytail !== 0) {
1482 const len = scale(abtlen, abt, cdytail, _16c);
1483 finlen = finadd(finlen, sum(
1484 scale(cytablen, cytab, cdytail, _16$1), _16$1,
1485 scale(len, _16c, 2 * cdy, _32), _32, _48$1), _48$1);
1486
1487 const len2 = scale(abttlen, abtt, cdytail, _8$1);
1488 finlen = finadd(finlen, sum_three(
1489 scale(len2, _8$1, 2 * cdy, _16$1), _16$1,
1490 scale(len2, _8$1, cdytail, _16b), _16b,
1491 scale(len, _16c, cdytail, _32), _32, _32b, _64), _64);
1492 }
1493 }
1494
1495 return fin$1[finlen - 1];
1496}
1497
1498function incircle(ax, ay, bx, by, cx, cy, dx, dy) {
1499 const adx = ax - dx;
1500 const bdx = bx - dx;
1501 const cdx = cx - dx;
1502 const ady = ay - dy;
1503 const bdy = by - dy;
1504 const cdy = cy - dy;
1505
1506 const bdxcdy = bdx * cdy;
1507 const cdxbdy = cdx * bdy;
1508 const alift = adx * adx + ady * ady;
1509
1510 const cdxady = cdx * ady;
1511 const adxcdy = adx * cdy;
1512 const blift = bdx * bdx + bdy * bdy;
1513
1514 const adxbdy = adx * bdy;
1515 const bdxady = bdx * ady;
1516 const clift = cdx * cdx + cdy * cdy;
1517
1518 const det =
1519 alift * (bdxcdy - cdxbdy) +
1520 blift * (cdxady - adxcdy) +
1521 clift * (adxbdy - bdxady);
1522
1523 const permanent =
1524 (Math.abs(bdxcdy) + Math.abs(cdxbdy)) * alift +
1525 (Math.abs(cdxady) + Math.abs(adxcdy)) * blift +
1526 (Math.abs(adxbdy) + Math.abs(bdxady)) * clift;
1527
1528 const errbound = iccerrboundA * permanent;
1529
1530 if (det > errbound || -det > errbound) {
1531 return det;
1532 }
1533 return incircleadapt(ax, ay, bx, by, cx, cy, dx, dy, permanent);
1534}
1535
1536function incirclefast(ax, ay, bx, by, cx, cy, dx, dy) {
1537 const adx = ax - dx;
1538 const ady = ay - dy;
1539 const bdx = bx - dx;
1540 const bdy = by - dy;
1541 const cdx = cx - dx;
1542 const cdy = cy - dy;
1543
1544 const abdet = adx * bdy - bdx * ady;
1545 const bcdet = bdx * cdy - cdx * bdy;
1546 const cadet = cdx * ady - adx * cdy;
1547 const alift = adx * adx + ady * ady;
1548 const blift = bdx * bdx + bdy * bdy;
1549 const clift = cdx * cdx + cdy * cdy;
1550
1551 return alift * bcdet + blift * cadet + clift * abdet;
1552}
1553
1554const isperrboundA = (16 + 224 * epsilon) * epsilon;
1555const isperrboundB = (5 + 72 * epsilon) * epsilon;
1556const isperrboundC = (71 + 1408 * epsilon) * epsilon * epsilon;
1557
1558const ab = vec(4);
1559const bc = vec(4);
1560const cd = vec(4);
1561const de = vec(4);
1562const ea = vec(4);
1563const ac = vec(4);
1564const bd = vec(4);
1565const ce = vec(4);
1566const da = vec(4);
1567const eb = vec(4);
1568
1569const abc = vec(24);
1570const bcd = vec(24);
1571const cde = vec(24);
1572const dea = vec(24);
1573const eab = vec(24);
1574const abd = vec(24);
1575const bce = vec(24);
1576const cda = vec(24);
1577const deb = vec(24);
1578const eac = vec(24);
1579
1580const adet = vec(1152);
1581const bdet = vec(1152);
1582const cdet = vec(1152);
1583const ddet = vec(1152);
1584const edet = vec(1152);
1585const abdet = vec(2304);
1586const cddet = vec(2304);
1587const cdedet = vec(3456);
1588const deter = vec(5760);
1589
1590const _8 = vec(8);
1591const _8b = vec(8);
1592const _8c = vec(8);
1593const _16 = vec(16);
1594const _24 = vec(24);
1595const _48 = vec(48);
1596const _48b = vec(48);
1597const _96 = vec(96);
1598const _192 = vec(192);
1599const _384x = vec(384);
1600const _384y = vec(384);
1601const _384z = vec(384);
1602const _768 = vec(768);
1603
1604function sum_three_scale(a, b, c, az, bz, cz, out) {
1605 return sum_three(
1606 scale(4, a, az, _8), _8,
1607 scale(4, b, bz, _8b), _8b,
1608 scale(4, c, cz, _8c), _8c, _16, out);
1609}
1610
1611function liftexact(alen, a, blen, b, clen, c, dlen, d, x, y, z, out) {
1612 const len = sum(
1613 sum(alen, a, blen, b, _48), _48,
1614 negate(sum(clen, c, dlen, d, _48b), _48b), _48b, _96);
1615
1616 return sum_three(
1617 scale(scale(len, _96, x, _192), _192, x, _384x), _384x,
1618 scale(scale(len, _96, y, _192), _192, y, _384y), _384y,
1619 scale(scale(len, _96, z, _192), _192, z, _384z), _384z, _768, out);
1620}
1621
1622function insphereexact(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez) {
1623 let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0, u3;
1624
1625 s1 = ax * by;
1626 c = splitter * ax;
1627 ahi = c - (c - ax);
1628 alo = ax - ahi;
1629 c = splitter * by;
1630 bhi = c - (c - by);
1631 blo = by - bhi;
1632 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1633 t1 = bx * ay;
1634 c = splitter * bx;
1635 ahi = c - (c - bx);
1636 alo = bx - ahi;
1637 c = splitter * ay;
1638 bhi = c - (c - ay);
1639 blo = ay - bhi;
1640 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1641 _i = s0 - t0;
1642 bvirt = s0 - _i;
1643 ab[0] = s0 - (_i + bvirt) + (bvirt - t0);
1644 _j = s1 + _i;
1645 bvirt = _j - s1;
1646 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1647 _i = _0 - t1;
1648 bvirt = _0 - _i;
1649 ab[1] = _0 - (_i + bvirt) + (bvirt - t1);
1650 u3 = _j + _i;
1651 bvirt = u3 - _j;
1652 ab[2] = _j - (u3 - bvirt) + (_i - bvirt);
1653 ab[3] = u3;
1654 s1 = bx * cy;
1655 c = splitter * bx;
1656 ahi = c - (c - bx);
1657 alo = bx - ahi;
1658 c = splitter * cy;
1659 bhi = c - (c - cy);
1660 blo = cy - bhi;
1661 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1662 t1 = cx * by;
1663 c = splitter * cx;
1664 ahi = c - (c - cx);
1665 alo = cx - ahi;
1666 c = splitter * by;
1667 bhi = c - (c - by);
1668 blo = by - bhi;
1669 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1670 _i = s0 - t0;
1671 bvirt = s0 - _i;
1672 bc[0] = s0 - (_i + bvirt) + (bvirt - t0);
1673 _j = s1 + _i;
1674 bvirt = _j - s1;
1675 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1676 _i = _0 - t1;
1677 bvirt = _0 - _i;
1678 bc[1] = _0 - (_i + bvirt) + (bvirt - t1);
1679 u3 = _j + _i;
1680 bvirt = u3 - _j;
1681 bc[2] = _j - (u3 - bvirt) + (_i - bvirt);
1682 bc[3] = u3;
1683 s1 = cx * dy;
1684 c = splitter * cx;
1685 ahi = c - (c - cx);
1686 alo = cx - ahi;
1687 c = splitter * dy;
1688 bhi = c - (c - dy);
1689 blo = dy - bhi;
1690 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1691 t1 = dx * cy;
1692 c = splitter * dx;
1693 ahi = c - (c - dx);
1694 alo = dx - ahi;
1695 c = splitter * cy;
1696 bhi = c - (c - cy);
1697 blo = cy - bhi;
1698 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1699 _i = s0 - t0;
1700 bvirt = s0 - _i;
1701 cd[0] = s0 - (_i + bvirt) + (bvirt - t0);
1702 _j = s1 + _i;
1703 bvirt = _j - s1;
1704 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1705 _i = _0 - t1;
1706 bvirt = _0 - _i;
1707 cd[1] = _0 - (_i + bvirt) + (bvirt - t1);
1708 u3 = _j + _i;
1709 bvirt = u3 - _j;
1710 cd[2] = _j - (u3 - bvirt) + (_i - bvirt);
1711 cd[3] = u3;
1712 s1 = dx * ey;
1713 c = splitter * dx;
1714 ahi = c - (c - dx);
1715 alo = dx - ahi;
1716 c = splitter * ey;
1717 bhi = c - (c - ey);
1718 blo = ey - bhi;
1719 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1720 t1 = ex * dy;
1721 c = splitter * ex;
1722 ahi = c - (c - ex);
1723 alo = ex - ahi;
1724 c = splitter * dy;
1725 bhi = c - (c - dy);
1726 blo = dy - bhi;
1727 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1728 _i = s0 - t0;
1729 bvirt = s0 - _i;
1730 de[0] = s0 - (_i + bvirt) + (bvirt - t0);
1731 _j = s1 + _i;
1732 bvirt = _j - s1;
1733 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1734 _i = _0 - t1;
1735 bvirt = _0 - _i;
1736 de[1] = _0 - (_i + bvirt) + (bvirt - t1);
1737 u3 = _j + _i;
1738 bvirt = u3 - _j;
1739 de[2] = _j - (u3 - bvirt) + (_i - bvirt);
1740 de[3] = u3;
1741 s1 = ex * ay;
1742 c = splitter * ex;
1743 ahi = c - (c - ex);
1744 alo = ex - ahi;
1745 c = splitter * ay;
1746 bhi = c - (c - ay);
1747 blo = ay - bhi;
1748 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1749 t1 = ax * ey;
1750 c = splitter * ax;
1751 ahi = c - (c - ax);
1752 alo = ax - ahi;
1753 c = splitter * ey;
1754 bhi = c - (c - ey);
1755 blo = ey - bhi;
1756 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1757 _i = s0 - t0;
1758 bvirt = s0 - _i;
1759 ea[0] = s0 - (_i + bvirt) + (bvirt - t0);
1760 _j = s1 + _i;
1761 bvirt = _j - s1;
1762 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1763 _i = _0 - t1;
1764 bvirt = _0 - _i;
1765 ea[1] = _0 - (_i + bvirt) + (bvirt - t1);
1766 u3 = _j + _i;
1767 bvirt = u3 - _j;
1768 ea[2] = _j - (u3 - bvirt) + (_i - bvirt);
1769 ea[3] = u3;
1770 s1 = ax * cy;
1771 c = splitter * ax;
1772 ahi = c - (c - ax);
1773 alo = ax - ahi;
1774 c = splitter * cy;
1775 bhi = c - (c - cy);
1776 blo = cy - bhi;
1777 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1778 t1 = cx * ay;
1779 c = splitter * cx;
1780 ahi = c - (c - cx);
1781 alo = cx - ahi;
1782 c = splitter * ay;
1783 bhi = c - (c - ay);
1784 blo = ay - bhi;
1785 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1786 _i = s0 - t0;
1787 bvirt = s0 - _i;
1788 ac[0] = s0 - (_i + bvirt) + (bvirt - t0);
1789 _j = s1 + _i;
1790 bvirt = _j - s1;
1791 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1792 _i = _0 - t1;
1793 bvirt = _0 - _i;
1794 ac[1] = _0 - (_i + bvirt) + (bvirt - t1);
1795 u3 = _j + _i;
1796 bvirt = u3 - _j;
1797 ac[2] = _j - (u3 - bvirt) + (_i - bvirt);
1798 ac[3] = u3;
1799 s1 = bx * dy;
1800 c = splitter * bx;
1801 ahi = c - (c - bx);
1802 alo = bx - ahi;
1803 c = splitter * dy;
1804 bhi = c - (c - dy);
1805 blo = dy - bhi;
1806 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1807 t1 = dx * by;
1808 c = splitter * dx;
1809 ahi = c - (c - dx);
1810 alo = dx - ahi;
1811 c = splitter * by;
1812 bhi = c - (c - by);
1813 blo = by - bhi;
1814 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1815 _i = s0 - t0;
1816 bvirt = s0 - _i;
1817 bd[0] = s0 - (_i + bvirt) + (bvirt - t0);
1818 _j = s1 + _i;
1819 bvirt = _j - s1;
1820 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1821 _i = _0 - t1;
1822 bvirt = _0 - _i;
1823 bd[1] = _0 - (_i + bvirt) + (bvirt - t1);
1824 u3 = _j + _i;
1825 bvirt = u3 - _j;
1826 bd[2] = _j - (u3 - bvirt) + (_i - bvirt);
1827 bd[3] = u3;
1828 s1 = cx * ey;
1829 c = splitter * cx;
1830 ahi = c - (c - cx);
1831 alo = cx - ahi;
1832 c = splitter * ey;
1833 bhi = c - (c - ey);
1834 blo = ey - bhi;
1835 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1836 t1 = ex * cy;
1837 c = splitter * ex;
1838 ahi = c - (c - ex);
1839 alo = ex - ahi;
1840 c = splitter * cy;
1841 bhi = c - (c - cy);
1842 blo = cy - bhi;
1843 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1844 _i = s0 - t0;
1845 bvirt = s0 - _i;
1846 ce[0] = s0 - (_i + bvirt) + (bvirt - t0);
1847 _j = s1 + _i;
1848 bvirt = _j - s1;
1849 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1850 _i = _0 - t1;
1851 bvirt = _0 - _i;
1852 ce[1] = _0 - (_i + bvirt) + (bvirt - t1);
1853 u3 = _j + _i;
1854 bvirt = u3 - _j;
1855 ce[2] = _j - (u3 - bvirt) + (_i - bvirt);
1856 ce[3] = u3;
1857 s1 = dx * ay;
1858 c = splitter * dx;
1859 ahi = c - (c - dx);
1860 alo = dx - ahi;
1861 c = splitter * ay;
1862 bhi = c - (c - ay);
1863 blo = ay - bhi;
1864 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1865 t1 = ax * dy;
1866 c = splitter * ax;
1867 ahi = c - (c - ax);
1868 alo = ax - ahi;
1869 c = splitter * dy;
1870 bhi = c - (c - dy);
1871 blo = dy - bhi;
1872 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1873 _i = s0 - t0;
1874 bvirt = s0 - _i;
1875 da[0] = s0 - (_i + bvirt) + (bvirt - t0);
1876 _j = s1 + _i;
1877 bvirt = _j - s1;
1878 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1879 _i = _0 - t1;
1880 bvirt = _0 - _i;
1881 da[1] = _0 - (_i + bvirt) + (bvirt - t1);
1882 u3 = _j + _i;
1883 bvirt = u3 - _j;
1884 da[2] = _j - (u3 - bvirt) + (_i - bvirt);
1885 da[3] = u3;
1886 s1 = ex * by;
1887 c = splitter * ex;
1888 ahi = c - (c - ex);
1889 alo = ex - ahi;
1890 c = splitter * by;
1891 bhi = c - (c - by);
1892 blo = by - bhi;
1893 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1894 t1 = bx * ey;
1895 c = splitter * bx;
1896 ahi = c - (c - bx);
1897 alo = bx - ahi;
1898 c = splitter * ey;
1899 bhi = c - (c - ey);
1900 blo = ey - bhi;
1901 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1902 _i = s0 - t0;
1903 bvirt = s0 - _i;
1904 eb[0] = s0 - (_i + bvirt) + (bvirt - t0);
1905 _j = s1 + _i;
1906 bvirt = _j - s1;
1907 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1908 _i = _0 - t1;
1909 bvirt = _0 - _i;
1910 eb[1] = _0 - (_i + bvirt) + (bvirt - t1);
1911 u3 = _j + _i;
1912 bvirt = u3 - _j;
1913 eb[2] = _j - (u3 - bvirt) + (_i - bvirt);
1914 eb[3] = u3;
1915
1916 const abclen = sum_three_scale(ab, bc, ac, cz, az, -bz, abc);
1917 const bcdlen = sum_three_scale(bc, cd, bd, dz, bz, -cz, bcd);
1918 const cdelen = sum_three_scale(cd, de, ce, ez, cz, -dz, cde);
1919 const dealen = sum_three_scale(de, ea, da, az, dz, -ez, dea);
1920 const eablen = sum_three_scale(ea, ab, eb, bz, ez, -az, eab);
1921 const abdlen = sum_three_scale(ab, bd, da, dz, az, bz, abd);
1922 const bcelen = sum_three_scale(bc, ce, eb, ez, bz, cz, bce);
1923 const cdalen = sum_three_scale(cd, da, ac, az, cz, dz, cda);
1924 const deblen = sum_three_scale(de, eb, bd, bz, dz, ez, deb);
1925 const eaclen = sum_three_scale(ea, ac, ce, cz, ez, az, eac);
1926
1927 const deterlen = sum_three(
1928 liftexact(cdelen, cde, bcelen, bce, deblen, deb, bcdlen, bcd, ax, ay, az, adet), adet,
1929 liftexact(dealen, dea, cdalen, cda, eaclen, eac, cdelen, cde, bx, by, bz, bdet), bdet,
1930 sum_three(
1931 liftexact(eablen, eab, deblen, deb, abdlen, abd, dealen, dea, cx, cy, cz, cdet), cdet,
1932 liftexact(abclen, abc, eaclen, eac, bcelen, bce, eablen, eab, dx, dy, dz, ddet), ddet,
1933 liftexact(bcdlen, bcd, abdlen, abd, cdalen, cda, abclen, abc, ex, ey, ez, edet), edet, cddet, cdedet), cdedet, abdet, deter);
1934
1935 return deter[deterlen - 1];
1936}
1937
1938const xdet = vec(96);
1939const ydet = vec(96);
1940const zdet = vec(96);
1941const fin = vec(1152);
1942
1943function liftadapt(a, b, c, az, bz, cz, x, y, z, out) {
1944 const len = sum_three_scale(a, b, c, az, bz, cz, _24);
1945 return sum_three(
1946 scale(scale(len, _24, x, _48), _48, x, xdet), xdet,
1947 scale(scale(len, _24, y, _48), _48, y, ydet), ydet,
1948 scale(scale(len, _24, z, _48), _48, z, zdet), zdet, _192, out);
1949}
1950
1951function insphereadapt(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez, permanent) {
1952 let ab3, bc3, cd3, da3, ac3, bd3;
1953
1954 let aextail, bextail, cextail, dextail;
1955 let aeytail, beytail, ceytail, deytail;
1956 let aeztail, beztail, ceztail, deztail;
1957
1958 let bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0;
1959
1960 const aex = ax - ex;
1961 const bex = bx - ex;
1962 const cex = cx - ex;
1963 const dex = dx - ex;
1964 const aey = ay - ey;
1965 const bey = by - ey;
1966 const cey = cy - ey;
1967 const dey = dy - ey;
1968 const aez = az - ez;
1969 const bez = bz - ez;
1970 const cez = cz - ez;
1971 const dez = dz - ez;
1972
1973 s1 = aex * bey;
1974 c = splitter * aex;
1975 ahi = c - (c - aex);
1976 alo = aex - ahi;
1977 c = splitter * bey;
1978 bhi = c - (c - bey);
1979 blo = bey - bhi;
1980 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
1981 t1 = bex * aey;
1982 c = splitter * bex;
1983 ahi = c - (c - bex);
1984 alo = bex - ahi;
1985 c = splitter * aey;
1986 bhi = c - (c - aey);
1987 blo = aey - bhi;
1988 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
1989 _i = s0 - t0;
1990 bvirt = s0 - _i;
1991 ab[0] = s0 - (_i + bvirt) + (bvirt - t0);
1992 _j = s1 + _i;
1993 bvirt = _j - s1;
1994 _0 = s1 - (_j - bvirt) + (_i - bvirt);
1995 _i = _0 - t1;
1996 bvirt = _0 - _i;
1997 ab[1] = _0 - (_i + bvirt) + (bvirt - t1);
1998 ab3 = _j + _i;
1999 bvirt = ab3 - _j;
2000 ab[2] = _j - (ab3 - bvirt) + (_i - bvirt);
2001 ab[3] = ab3;
2002 s1 = bex * cey;
2003 c = splitter * bex;
2004 ahi = c - (c - bex);
2005 alo = bex - ahi;
2006 c = splitter * cey;
2007 bhi = c - (c - cey);
2008 blo = cey - bhi;
2009 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
2010 t1 = cex * bey;
2011 c = splitter * cex;
2012 ahi = c - (c - cex);
2013 alo = cex - ahi;
2014 c = splitter * bey;
2015 bhi = c - (c - bey);
2016 blo = bey - bhi;
2017 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
2018 _i = s0 - t0;
2019 bvirt = s0 - _i;
2020 bc[0] = s0 - (_i + bvirt) + (bvirt - t0);
2021 _j = s1 + _i;
2022 bvirt = _j - s1;
2023 _0 = s1 - (_j - bvirt) + (_i - bvirt);
2024 _i = _0 - t1;
2025 bvirt = _0 - _i;
2026 bc[1] = _0 - (_i + bvirt) + (bvirt - t1);
2027 bc3 = _j + _i;
2028 bvirt = bc3 - _j;
2029 bc[2] = _j - (bc3 - bvirt) + (_i - bvirt);
2030 bc[3] = bc3;
2031 s1 = cex * dey;
2032 c = splitter * cex;
2033 ahi = c - (c - cex);
2034 alo = cex - ahi;
2035 c = splitter * dey;
2036 bhi = c - (c - dey);
2037 blo = dey - bhi;
2038 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
2039 t1 = dex * cey;
2040 c = splitter * dex;
2041 ahi = c - (c - dex);
2042 alo = dex - ahi;
2043 c = splitter * cey;
2044 bhi = c - (c - cey);
2045 blo = cey - bhi;
2046 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
2047 _i = s0 - t0;
2048 bvirt = s0 - _i;
2049 cd[0] = s0 - (_i + bvirt) + (bvirt - t0);
2050 _j = s1 + _i;
2051 bvirt = _j - s1;
2052 _0 = s1 - (_j - bvirt) + (_i - bvirt);
2053 _i = _0 - t1;
2054 bvirt = _0 - _i;
2055 cd[1] = _0 - (_i + bvirt) + (bvirt - t1);
2056 cd3 = _j + _i;
2057 bvirt = cd3 - _j;
2058 cd[2] = _j - (cd3 - bvirt) + (_i - bvirt);
2059 cd[3] = cd3;
2060 s1 = dex * aey;
2061 c = splitter * dex;
2062 ahi = c - (c - dex);
2063 alo = dex - ahi;
2064 c = splitter * aey;
2065 bhi = c - (c - aey);
2066 blo = aey - bhi;
2067 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
2068 t1 = aex * dey;
2069 c = splitter * aex;
2070 ahi = c - (c - aex);
2071 alo = aex - ahi;
2072 c = splitter * dey;
2073 bhi = c - (c - dey);
2074 blo = dey - bhi;
2075 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
2076 _i = s0 - t0;
2077 bvirt = s0 - _i;
2078 da[0] = s0 - (_i + bvirt) + (bvirt - t0);
2079 _j = s1 + _i;
2080 bvirt = _j - s1;
2081 _0 = s1 - (_j - bvirt) + (_i - bvirt);
2082 _i = _0 - t1;
2083 bvirt = _0 - _i;
2084 da[1] = _0 - (_i + bvirt) + (bvirt - t1);
2085 da3 = _j + _i;
2086 bvirt = da3 - _j;
2087 da[2] = _j - (da3 - bvirt) + (_i - bvirt);
2088 da[3] = da3;
2089 s1 = aex * cey;
2090 c = splitter * aex;
2091 ahi = c - (c - aex);
2092 alo = aex - ahi;
2093 c = splitter * cey;
2094 bhi = c - (c - cey);
2095 blo = cey - bhi;
2096 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
2097 t1 = cex * aey;
2098 c = splitter * cex;
2099 ahi = c - (c - cex);
2100 alo = cex - ahi;
2101 c = splitter * aey;
2102 bhi = c - (c - aey);
2103 blo = aey - bhi;
2104 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
2105 _i = s0 - t0;
2106 bvirt = s0 - _i;
2107 ac[0] = s0 - (_i + bvirt) + (bvirt - t0);
2108 _j = s1 + _i;
2109 bvirt = _j - s1;
2110 _0 = s1 - (_j - bvirt) + (_i - bvirt);
2111 _i = _0 - t1;
2112 bvirt = _0 - _i;
2113 ac[1] = _0 - (_i + bvirt) + (bvirt - t1);
2114 ac3 = _j + _i;
2115 bvirt = ac3 - _j;
2116 ac[2] = _j - (ac3 - bvirt) + (_i - bvirt);
2117 ac[3] = ac3;
2118 s1 = bex * dey;
2119 c = splitter * bex;
2120 ahi = c - (c - bex);
2121 alo = bex - ahi;
2122 c = splitter * dey;
2123 bhi = c - (c - dey);
2124 blo = dey - bhi;
2125 s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
2126 t1 = dex * bey;
2127 c = splitter * dex;
2128 ahi = c - (c - dex);
2129 alo = dex - ahi;
2130 c = splitter * bey;
2131 bhi = c - (c - bey);
2132 blo = bey - bhi;
2133 t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
2134 _i = s0 - t0;
2135 bvirt = s0 - _i;
2136 bd[0] = s0 - (_i + bvirt) + (bvirt - t0);
2137 _j = s1 + _i;
2138 bvirt = _j - s1;
2139 _0 = s1 - (_j - bvirt) + (_i - bvirt);
2140 _i = _0 - t1;
2141 bvirt = _0 - _i;
2142 bd[1] = _0 - (_i + bvirt) + (bvirt - t1);
2143 bd3 = _j + _i;
2144 bvirt = bd3 - _j;
2145 bd[2] = _j - (bd3 - bvirt) + (_i - bvirt);
2146 bd[3] = bd3;
2147
2148 const finlen = sum(
2149 sum(
2150 negate(liftadapt(bc, cd, bd, dez, bez, -cez, aex, aey, aez, adet), adet), adet,
2151 liftadapt(cd, da, ac, aez, cez, dez, bex, bey, bez, bdet), bdet, abdet), abdet,
2152 sum(
2153 negate(liftadapt(da, ab, bd, bez, dez, aez, cex, cey, cez, cdet), cdet), cdet,
2154 liftadapt(ab, bc, ac, cez, aez, -bez, dex, dey, dez, ddet), ddet, cddet), cddet, fin);
2155
2156 let det = estimate(finlen, fin);
2157 let errbound = isperrboundB * permanent;
2158 if (det >= errbound || -det >= errbound) {
2159 return det;
2160 }
2161
2162 bvirt = ax - aex;
2163 aextail = ax - (aex + bvirt) + (bvirt - ex);
2164 bvirt = ay - aey;
2165 aeytail = ay - (aey + bvirt) + (bvirt - ey);
2166 bvirt = az - aez;
2167 aeztail = az - (aez + bvirt) + (bvirt - ez);
2168 bvirt = bx - bex;
2169 bextail = bx - (bex + bvirt) + (bvirt - ex);
2170 bvirt = by - bey;
2171 beytail = by - (bey + bvirt) + (bvirt - ey);
2172 bvirt = bz - bez;
2173 beztail = bz - (bez + bvirt) + (bvirt - ez);
2174 bvirt = cx - cex;
2175 cextail = cx - (cex + bvirt) + (bvirt - ex);
2176 bvirt = cy - cey;
2177 ceytail = cy - (cey + bvirt) + (bvirt - ey);
2178 bvirt = cz - cez;
2179 ceztail = cz - (cez + bvirt) + (bvirt - ez);
2180 bvirt = dx - dex;
2181 dextail = dx - (dex + bvirt) + (bvirt - ex);
2182 bvirt = dy - dey;
2183 deytail = dy - (dey + bvirt) + (bvirt - ey);
2184 bvirt = dz - dez;
2185 deztail = dz - (dez + bvirt) + (bvirt - ez);
2186 if (aextail === 0 && aeytail === 0 && aeztail === 0 &&
2187 bextail === 0 && beytail === 0 && beztail === 0 &&
2188 cextail === 0 && ceytail === 0 && ceztail === 0 &&
2189 dextail === 0 && deytail === 0 && deztail === 0) {
2190 return det;
2191 }
2192
2193 errbound = isperrboundC * permanent + resulterrbound * Math.abs(det);
2194
2195 const abeps = (aex * beytail + bey * aextail) - (aey * bextail + bex * aeytail);
2196 const bceps = (bex * ceytail + cey * bextail) - (bey * cextail + cex * beytail);
2197 const cdeps = (cex * deytail + dey * cextail) - (cey * dextail + dex * ceytail);
2198 const daeps = (dex * aeytail + aey * dextail) - (dey * aextail + aex * deytail);
2199 const aceps = (aex * ceytail + cey * aextail) - (aey * cextail + cex * aeytail);
2200 const bdeps = (bex * deytail + dey * bextail) - (bey * dextail + dex * beytail);
2201 det +=
2202 (((bex * bex + bey * bey + bez * bez) * ((cez * daeps + dez * aceps + aez * cdeps) +
2203 (ceztail * da3 + deztail * ac3 + aeztail * cd3)) + (dex * dex + dey * dey + dez * dez) *
2204 ((aez * bceps - bez * aceps + cez * abeps) + (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) -
2205 ((aex * aex + aey * aey + aez * aez) * ((bez * cdeps - cez * bdeps + dez * bceps) +
2206 (beztail * cd3 - ceztail * bd3 + deztail * bc3)) + (cex * cex + cey * cey + cez * cez) *
2207 ((dez * abeps + aez * bdeps + bez * daeps) + (deztail * ab3 + aeztail * bd3 + beztail * da3)))) +
2208 2 * (((bex * bextail + bey * beytail + bez * beztail) * (cez * da3 + dez * ac3 + aez * cd3) +
2209 (dex * dextail + dey * deytail + dez * deztail) * (aez * bc3 - bez * ac3 + cez * ab3)) -
2210 ((aex * aextail + aey * aeytail + aez * aeztail) * (bez * cd3 - cez * bd3 + dez * bc3) +
2211 (cex * cextail + cey * ceytail + cez * ceztail) * (dez * ab3 + aez * bd3 + bez * da3)));
2212
2213 if (det >= errbound || -det >= errbound) {
2214 return det;
2215 }
2216
2217 return insphereexact(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez);
2218}
2219
2220function insphere(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez) {
2221 const aex = ax - ex;
2222 const bex = bx - ex;
2223 const cex = cx - ex;
2224 const dex = dx - ex;
2225 const aey = ay - ey;
2226 const bey = by - ey;
2227 const cey = cy - ey;
2228 const dey = dy - ey;
2229 const aez = az - ez;
2230 const bez = bz - ez;
2231 const cez = cz - ez;
2232 const dez = dz - ez;
2233
2234 const aexbey = aex * bey;
2235 const bexaey = bex * aey;
2236 const ab = aexbey - bexaey;
2237 const bexcey = bex * cey;
2238 const cexbey = cex * bey;
2239 const bc = bexcey - cexbey;
2240 const cexdey = cex * dey;
2241 const dexcey = dex * cey;
2242 const cd = cexdey - dexcey;
2243 const dexaey = dex * aey;
2244 const aexdey = aex * dey;
2245 const da = dexaey - aexdey;
2246 const aexcey = aex * cey;
2247 const cexaey = cex * aey;
2248 const ac = aexcey - cexaey;
2249 const bexdey = bex * dey;
2250 const dexbey = dex * bey;
2251 const bd = bexdey - dexbey;
2252
2253 const alift = aex * aex + aey * aey + aez * aez;
2254 const blift = bex * bex + bey * bey + bez * bez;
2255 const clift = cex * cex + cey * cey + cez * cez;
2256 const dlift = dex * dex + dey * dey + dez * dez;
2257
2258 const det =
2259 (clift * (dez * ab + aez * bd + bez * da) - dlift * (aez * bc - bez * ac + cez * ab)) +
2260 (alift * (bez * cd - cez * bd + dez * bc) - blift * (cez * da + dez * ac + aez * cd));
2261
2262 const aezplus = Math.abs(aez);
2263 const bezplus = Math.abs(bez);
2264 const cezplus = Math.abs(cez);
2265 const dezplus = Math.abs(dez);
2266 const aexbeyplus = Math.abs(aexbey) + Math.abs(bexaey);
2267 const bexceyplus = Math.abs(bexcey) + Math.abs(cexbey);
2268 const cexdeyplus = Math.abs(cexdey) + Math.abs(dexcey);
2269 const dexaeyplus = Math.abs(dexaey) + Math.abs(aexdey);
2270 const aexceyplus = Math.abs(aexcey) + Math.abs(cexaey);
2271 const bexdeyplus = Math.abs(bexdey) + Math.abs(dexbey);
2272 const permanent =
2273 (cexdeyplus * bezplus + bexdeyplus * cezplus + bexceyplus * dezplus) * alift +
2274 (dexaeyplus * cezplus + aexceyplus * dezplus + cexdeyplus * aezplus) * blift +
2275 (aexbeyplus * dezplus + bexdeyplus * aezplus + dexaeyplus * bezplus) * clift +
2276 (bexceyplus * aezplus + aexceyplus * bezplus + aexbeyplus * cezplus) * dlift;
2277
2278 const errbound = isperrboundA * permanent;
2279 if (det > errbound || -det > errbound) {
2280 return det;
2281 }
2282 return -insphereadapt(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, ex, ey, ez, permanent);
2283}
2284
2285function inspherefast(pax, pay, paz, pbx, pby, pbz, pcx, pcy, pcz, pdx, pdy, pdz, pex, pey, pez) {
2286 const aex = pax - pex;
2287 const bex = pbx - pex;
2288 const cex = pcx - pex;
2289 const dex = pdx - pex;
2290 const aey = pay - pey;
2291 const bey = pby - pey;
2292 const cey = pcy - pey;
2293 const dey = pdy - pey;
2294 const aez = paz - pez;
2295 const bez = pbz - pez;
2296 const cez = pcz - pez;
2297 const dez = pdz - pez;
2298
2299 const ab = aex * bey - bex * aey;
2300 const bc = bex * cey - cex * bey;
2301 const cd = cex * dey - dex * cey;
2302 const da = dex * aey - aex * dey;
2303 const ac = aex * cey - cex * aey;
2304 const bd = bex * dey - dex * bey;
2305
2306 const abc = aez * bc - bez * ac + cez * ab;
2307 const bcd = bez * cd - cez * bd + dez * bc;
2308 const cda = cez * da + dez * ac + aez * cd;
2309 const dab = dez * ab + aez * bd + bez * da;
2310
2311 const alift = aex * aex + aey * aey + aez * aez;
2312 const blift = bex * bex + bey * bey + bez * bez;
2313 const clift = cex * cex + cey * cey + cez * cez;
2314 const dlift = dex * dex + dey * dey + dez * dez;
2315
2316 return (clift * dab - dlift * abc) + (alift * bcd - blift * cda);
2317}
2318
2319exports.incircle = incircle;
2320exports.incirclefast = incirclefast;
2321exports.insphere = insphere;
2322exports.inspherefast = inspherefast;
2323exports.orient2d = orient2d;
2324exports.orient2dfast = orient2dfast;
2325exports.orient3d = orient3d;
2326exports.orient3dfast = orient3dfast;
2327
2328}));
Note: See TracBrowser for help on using the repository browser.