1 /*****************************************************************************/
2 /* */
3 /* Routines for Arbitrary Precision Floating-point Arithmetic */
4 /* and Fast Robust Geometric Predicates */
5 /* (predicates.cpp) */
6 /* */
7 /* May 18, 1996 */
8 /* */
9 /* Placed in the public domain by */
10 /* Jonathan Richard Shewchuk */
11 /* School of Computer Science */
12 /* Carnegie Mellon University */
13 /* 5000 Forbes Avenue */
14 /* Pittsburgh, Pennsylvania 15213-3891 */
15 /* jrs@cs.cmu.edu */
16 /* */
17 /* This file contains C implementation of algorithms for exact addition */
18 /* and multiplication of floating-point numbers, and predicates for */
19 /* robustly performing the orientation and incircle tests used in */
20 /* computational geometry. The algorithms and underlying theory are */
21 /* described in Jonathan Richard Shewchuk. "Adaptive Precision Floating- */
22 /* Point Arithmetic and Fast Robust Geometric Predicates." Technical */
23 /* Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon */
24 /* University, Pittsburgh, Pennsylvania, May 1996. (Submitted to */
25 /* Discrete & Computational Geometry.) */
26 /* */
27 /* This file, the paper listed above, and other information are available */
28 /* from the Web page http://www.cs.cmu.edu/~quake/robust.html . */
29 /* */
30 /*****************************************************************************/
31
32 /*****************************************************************************/
33 /* */
34 /* Using this code: */
35 /* */
36 /* First, read the short or long version of the paper (from the Web page */
37 /* above). */
38 /* */
39 /* Be sure to call exactinit() once, before calling any of the arithmetic */
40 /* functions or geometric predicates. Also be sure to turn on the */
41 /* optimizer when compiling this file. */
42 /* */
43 /* */
44 /* Several geometric predicates are defined. Their parameters are all */
45 /* points. Each point is an array of two or three floating-point */
46 /* numbers. The geometric predicates, described in the papers, are */
47 /* */
48 /* incircle(pa, pb, pc, pd) */
49 /* incirclefast(pa, pb, pc, pd) */
50 /* */
51 /* Those with suffix "fast" are approximate, non-robust versions. Those */
52 /* without the suffix are adaptive precision, robust versions. There */
53 /* are also versions with the suffices "exact" and "slow", which are */
54 /* non-adaptive, exact arithmetic versions, which I use only for timings */
55 /* in my arithmetic papers. */
56 /* */
57 /* */
58 /* An expansion is represented by an array of floating-point numbers, */
59 /* sorted from smallest to largest magnitude (possibly with interspersed */
60 /* zeros). The length of each expansion is stored as a separate integer, */
61 /* and each arithmetic function returns an integer which is the length */
62 /* of the expansion it created. */
63 /* */
64 /* Several arithmetic functions are defined. Their parameters are */
65 /* */
66 /* e, f Input expansions */
67 /* elen, flen Lengths of input expansions (must be >= 1) */
68 /* h Output expansion */
69 /* b Input scalar */
70 /* */
71 /* The arithmetic functions are */
72 /* */
73 /* expansion_sum(elen, e, flen, f, h) */
74 /* expansion_sum_zeroelim1(elen, e, flen, f, h) */
75 /* expansion_sum_zeroelim2(elen, e, flen, f, h) */
76 /* fast_expansion_sum(elen, e, flen, f, h) */
77 /* fast_expansion_sum_zeroelim(elen, e, flen, f, h) */
78 /* linear_expansion_sum(elen, e, flen, f, h) */
79 /* linear_expansion_sum_zeroelim(elen, e, flen, f, h) */
80 /* scale_expansion(elen, e, b, h) */
81 /* scale_expansion_zeroelim(elen, e, b, h) */
82 /* compress(elen, e, h) */
83 /* */
84 /* All of these are described in the long version of the paper; some are */
85 /* described in the short version. All return an integer that is the */
86 /* length of h. Those with suffix _zeroelim perform zero elimination, */
87 /* and are recommended over their counterparts. The procedure */
88 /* fast_expansion_sum_zeroelim() (or linear_expansion_sum_zeroelim() on */
89 /* processors that do not use the round-to-even tiebreaking rule) is */
90 /* recommended over expansion_sum_zeroelim(). Each procedure has a */
91 /* little note next to it (in the code below) that tells you whether or */
92 /* not the output expansion may be the same array as one of the input */
93 /* expansions. */
94 /* */
95 /* */
96 /* If you look around below, you'll also find macros for a bunch of */
97 /* simple unrolled arithmetic operations, and procedures for printing */
98 /* expansions (commented out because they don't work with all C */
99 /* compilers) and for generating random floating-point numbers whose */
100 /* significand bits are all random. Most of the macros have undocumented */
101 /* requirements that certain of their parameters should not be the same */
102 /* variable; for safety, better to make sure all the parameters are */
103 /* distinct variables. Feel free to send email to jrs@cs.cmu.edu if you */
104 /* have questions. */
105 /* */
106 /*****************************************************************************/
107
108 #include <QtGlobal>
109
110 #include "../vmisc/diagnostic.h"
111
112 QT_WARNING_PUSH
113 QT_WARNING_DISABLE_MSVC(4701)
114 QT_WARNING_DISABLE_GCC("-Wold-style-cast")
115 QT_WARNING_DISABLE_GCC("-Wfloat-equal")
116 #if defined(__GNUC__) && (__GNUC__ * 100 + __GNUC_MINOR__) >= 408
117 QT_WARNING_DISABLE_GCC("-Wmaybe-uninitialized")
118 #else
119 QT_WARNING_DISABLE_GCC("-Wuninitialized")
120 #endif
121 QT_WARNING_DISABLE_CLANG("-Wold-style-cast")
122 QT_WARNING_DISABLE_CLANG("-Wmissing-variable-declarations")
123 QT_WARNING_DISABLE_CLANG("-Wfloat-equal")
124 QT_WARNING_DISABLE_CLANG("-Wmissing-prototypes")
125 QT_WARNING_DISABLE_CLANG("-Wconditional-uninitialized")
126
127 /* On some machines, the exact arithmetic routines might be defeated by the */
128 /* use of internal extended precision floating-point registers. Sometimes */
129 /* this problem can be fixed by defining certain values to be volatile, */
130 /* thus forcing them to be stored to memory and rounded off. This isn't */
131 /* a great solution, though, as it slows the arithmetic down. */
132 /* */
133 /* To try this out, write "#define INEXACT volatile" below. Normally, */
134 /* however, INEXACT should be defined to be nothing. ("#define INEXACT".) */
135
136 #define INEXACT /* Nothing */
137 /* #define INEXACT volatile */
138
139 #define REALPRINT doubleprint
140 #define REALRAND doublerand
141 #define NARROWRAND narrowdoublerand
142 #define UNIFORMRAND uniformdoublerand
143
144 /* Which of the following two methods of finding the absolute values is */
145 /* fastest is compiler-dependent. A few compilers can inline and optimize */
146 /* the fabs() call; but most will incur the overhead of a function call, */
147 /* which is disastrously slow. A faster way on IEEE machines might be to */
148 /* mask the appropriate bit, but that's difficult to do in C. */
149
150 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
151 /* #define Absolute(a) fabs(a) */
152
153
154 /* Many of the operations are broken up into two pieces, a main part that */
155 /* performs an approximate operation, and a "tail" that computes the */
156 /* roundoff error of that operation. */
157 /* */
158 /* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), */
159 /* Split(), and Two_Product() are all implemented as described in the */
160 /* reference. Each of these macros requires certain variables to be */
161 /* defined in the calling routine. The variables `bvirt', `c', `abig', */
162 /* `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because */
163 /* they store the result of an operation that may incur roundoff error. */
164 /* The input parameter `x' (or the highest numbered `x_' parameter) must */
165 /* also be declared `INEXACT'. */
166
167 #define Fast_Two_Sum_Tail(a, b, x, y) \
168 bvirt = x - a; \
169 y = b - bvirt
170
171 #define Fast_Two_Sum(a, b, x, y) \
172 x = (qreal) (a + b); \
173 Fast_Two_Sum_Tail(a, b, x, y)
174
175 #define Two_Sum_Tail(a, b, x, y) \
176 bvirt = (qreal) (x - a); \
177 avirt = x - bvirt; \
178 bround = b - bvirt; \
179 around = a - avirt; \
180 y = around + bround
181
182 #define Two_Sum(a, b, x, y) \
183 x = (qreal) (a + b); \
184 Two_Sum_Tail(a, b, x, y)
185
186 #define Two_Diff_Tail(a, b, x, y) \
187 bvirt = (qreal) (a - x); \
188 avirt = x + bvirt; \
189 bround = bvirt - b; \
190 around = a - avirt; \
191 y = around + bround
192
193 #define Two_Diff(a, b, x, y) \
194 x = (qreal) (a - b); \
195 Two_Diff_Tail(a, b, x, y)
196
197 #define Split(a, ahi, alo) \
198 c = (qreal) (splitter * a); \
199 abig = (qreal) (c - a); \
200 ahi = c - abig; \
201 alo = a - ahi
202
203 #define Two_Product_Tail(a, b, x, y) \
204 Split(a, ahi, alo); \
205 Split(b, bhi, blo); \
206 err1 = x - (ahi * bhi); \
207 err2 = err1 - (alo * bhi); \
208 err3 = err2 - (ahi * blo); \
209 y = (alo * blo) - err3
210
211 #define Two_Product(a, b, x, y) \
212 x = (qreal) (a * b); \
213 Two_Product_Tail(a, b, x, y)
214
215 /* Two_Product_Presplit() is Two_Product() where one of the inputs has */
216 /* already been split. Avoids redundant splitting. */
217
218 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
219 x = (qreal) (a * b); \
220 Split(a, ahi, alo); \
221 err1 = x - (ahi * bhi); \
222 err2 = err1 - (alo * bhi); \
223 err3 = err2 - (ahi * blo); \
224 y = (alo * blo) - err3
225
226 /* Square() can be done more quickly than Two_Product(). */
227
228 #define Square_Tail(a, x, y) \
229 Split(a, ahi, alo); \
230 err1 = x - (ahi * ahi); \
231 err3 = err1 - ((ahi + ahi) * alo); \
232 y = (alo * alo) - err3
233
234 #define Square(a, x, y) \
235 x = (qreal) (a * a); \
236 Square_Tail(a, x, y)
237
238 /* Macros for summing expansions of various fixed lengths. These are all */
239 /* unrolled versions of Expansion_Sum(). */
240
241 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
242 Two_Sum(a0, b , _i, x0); \
243 Two_Sum(a1, _i, x2, x1)
244
245 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
246 Two_Diff(a0, b , _i, x0); \
247 Two_Sum( a1, _i, x2, x1)
248
249 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
250 Two_One_Sum(a1, a0, b0, _j, _0, x0); \
251 Two_One_Sum(_j, _0, b1, x3, x2, x1)
252
253 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
254 Two_One_Diff(a1, a0, b0, _j, _0, x0); \
255 Two_One_Diff(_j, _0, b1, x3, x2, x1)
256
257 qreal splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */
258 qreal epsilon; /* = 2^(-p). Used to estimate roundoff errors. */
259 /* A set of coefficients used to calculate maximum roundoff errors. */
260 qreal resulterrbound;
261 qreal ccwerrboundA, ccwerrboundB, ccwerrboundC;
262 qreal o3derrboundA, o3derrboundB, o3derrboundC;
263 qreal iccerrboundA, iccerrboundB, iccerrboundC;
264 qreal isperrboundA, isperrboundB, isperrboundC;
265
266 /*****************************************************************************/
267 /* */
268 /* exactinit() Initialize the variables used for exact arithmetic. */
269 /* */
270 /* `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in */
271 /* floating-point arithmetic. `epsilon' bounds the relative roundoff */
272 /* error. It is used for floating-point error analysis. */
273 /* */
274 /* `splitter' is used to split floating-point numbers into two half- */
275 /* length significands for exact multiplication. */
276 /* */
277 /* I imagine that a highly optimizing compiler might be too smart for its */
278 /* own good, and somehow cause this routine to fail, if it pretends that */
279 /* floating-point arithmetic is too much like real arithmetic. */
280 /* */
281 /* Don't change this routine unless you fully understand it. */
282 /* */
283 /*****************************************************************************/
284
exactinit()285 void exactinit()
286 {
287 qreal half;
288 qreal check, lastcheck;
289 int every_other;
290
291 every_other = 1;
292 half = 0.5;
293 epsilon = 1.0;
294 splitter = 1.0;
295 check = 1.0;
296 /* Repeatedly divide `epsilon' by two until it is too small to add to */
297 /* one without causing roundoff. (Also check if the sum is equal to */
298 /* the previous sum, for machines that round up instead of using exact */
299 /* rounding. Not that this library will work on such machines anyway. */
300 do
301 {
302 lastcheck = check;
303 epsilon *= half;
304 if (every_other)
305 {
306 splitter *= 2.0;
307 }
308 every_other = !every_other;
309 check = 1.0 + epsilon;
310 } while ((check != 1.0) && (check != lastcheck));
311 splitter += 1.0;
312
313 /* Error bounds for orientation and incircle tests. */
314 resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
315 ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
316 ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
317 ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
318 o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
319 o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
320 o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
321 iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
322 iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
323 iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
324 isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
325 isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
326 isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
327 }
328
329 /*****************************************************************************/
330 /* */
331 /* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */
332 /* components from the output expansion. */
333 /* */
334 /* Sets h = e + f. See the long version of my paper for details. */
335 /* */
336 /* If round-to-even is used (as with IEEE 754), maintains the strongly */
337 /* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */
338 /* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */
339 /* properties. */
340 /* */
341 /*****************************************************************************/
342
fast_expansion_sum_zeroelim(int elen,qreal * e,int flen,qreal * f,qreal * h)343 int fast_expansion_sum_zeroelim(int elen, qreal *e, int flen, qreal *f, qreal *h) /* h cannot be e or f. */
344 {
345 qreal Q;
346 INEXACT qreal Qnew;
347 INEXACT qreal hh;
348 INEXACT qreal bvirt;
349 qreal avirt, bround, around;
350 int eindex, findex, hindex;
351 qreal enow, fnow;
352
353 enow = e[0];
354 fnow = f[0];
355 eindex = findex = 0;
356 if ((fnow > enow) == (fnow > -enow))
357 {
358 Q = enow;
359 enow = e[++eindex];
360 }
361 else
362 {
363 Q = fnow;
364 fnow = f[++findex];
365 }
366 hindex = 0;
367 if ((eindex < elen) && (findex < flen))
368 {
369 if ((fnow > enow) == (fnow > -enow))
370 {
371 Fast_Two_Sum(enow, Q, Qnew, hh);
372 enow = e[++eindex];
373 }
374 else
375 {
376 Fast_Two_Sum(fnow, Q, Qnew, hh);
377 fnow = f[++findex];
378 }
379 Q = Qnew;
380 if (hh != 0.0)
381 {
382 h[hindex++] = hh;
383 }
384 while ((eindex < elen) && (findex < flen))
385 {
386 if ((fnow > enow) == (fnow > -enow))
387 {
388 Two_Sum(Q, enow, Qnew, hh);
389 enow = e[++eindex];
390 }
391 else
392 {
393 Two_Sum(Q, fnow, Qnew, hh);
394 fnow = f[++findex];
395 }
396 Q = Qnew;
397 if (hh != 0.0)
398 {
399 h[hindex++] = hh;
400 }
401 }
402 }
403 while (eindex < elen)
404 {
405 Two_Sum(Q, enow, Qnew, hh);
406 enow = e[++eindex];
407 Q = Qnew;
408 if (hh != 0.0)
409 {
410 h[hindex++] = hh;
411 }
412 }
413 while (findex < flen)
414 {
415 Two_Sum(Q, fnow, Qnew, hh);
416 fnow = f[++findex];
417 Q = Qnew;
418 if (hh != 0.0)
419 {
420 h[hindex++] = hh;
421 }
422 }
423 if ((Q != 0.0) || (hindex == 0))
424 {
425 h[hindex++] = Q;
426 }
427 return hindex;
428 }
429
430 /*****************************************************************************/
431 /* */
432 /* scale_expansion_zeroelim() Multiply an expansion by a scalar, */
433 /* eliminating zero components from the */
434 /* output expansion. */
435 /* */
436 /* Sets h = be. See either version of my paper for details. */
437 /* */
438 /* Maintains the nonoverlapping property. If round-to-even is used (as */
439 /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
440 /* properties as well. (That is, if e has one of these properties, so */
441 /* will h.) */
442 /* */
443 /*****************************************************************************/
444
scale_expansion_zeroelim(int elen,qreal * e,qreal b,qreal * h)445 int scale_expansion_zeroelim(int elen, qreal *e, qreal b, qreal *h) /* e and h cannot be the same. */
446 {
447 INEXACT qreal Q, sum;
448 qreal hh;
449 INEXACT qreal product1;
450 qreal product0;
451 int eindex, hindex;
452 qreal enow;
453 INEXACT qreal bvirt;
454 qreal avirt, bround, around;
455 INEXACT qreal c;
456 INEXACT qreal abig;
457 qreal ahi, alo, bhi, blo;
458 qreal err1, err2, err3;
459
460 Split(b, bhi, blo);
461 Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
462 hindex = 0;
463 if (hh != 0) {
464 h[hindex++] = hh;
465 }
466 for (eindex = 1; eindex < elen; eindex++)
467 {
468 enow = e[eindex];
469 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
470 Two_Sum(Q, product0, sum, hh);
471 if (hh != 0)
472 {
473 h[hindex++] = hh;
474 }
475 Fast_Two_Sum(product1, sum, Q, hh);
476 if (hh != 0)
477 {
478 h[hindex++] = hh;
479 }
480 }
481 if ((Q != 0.0) || (hindex == 0))
482 {
483 h[hindex++] = Q;
484 }
485 return hindex;
486 }
487
488 /*****************************************************************************/
489 /* */
490 /* estimate() Produce a one-word estimate of an expansion's value. */
491 /* */
492 /* See either version of my paper for details. */
493 /* */
494 /*****************************************************************************/
495
estimate(int elen,qreal * e)496 qreal estimate(int elen, qreal *e)
497 {
498 qreal Q;
499 int eindex;
500
501 Q = e[0];
502 for (eindex = 1; eindex < elen; eindex++)
503 {
504 Q += e[eindex];
505 }
506 return Q;
507 }
508
incircleadapt(qreal * pa,qreal * pb,qreal * pc,qreal * pd,qreal permanent)509 qreal incircleadapt(qreal *pa, qreal *pb, qreal *pc, qreal *pd, qreal permanent)
510 {
511 INEXACT qreal adx, bdx, cdx, ady, bdy, cdy;
512 qreal det, errbound;
513
514 INEXACT qreal bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
515 qreal bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
516 qreal bc[4], ca[4], ab[4];
517 INEXACT qreal bc3, ca3, ab3;
518 qreal axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
519 int axbclen, axxbclen, aybclen, ayybclen, alen;
520 qreal bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
521 int bxcalen, bxxcalen, bycalen, byycalen, blen;
522 qreal cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
523 int cxablen, cxxablen, cyablen, cyyablen, clen;
524 qreal abdet[64];
525 int ablen;
526 qreal fin1[1152], fin2[1152];
527 qreal *finnow, *finother, *finswap;
528 int finlength;
529
530 qreal adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
531 INEXACT qreal adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
532 qreal adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
533 qreal aa[4], bb[4], cc[4];
534 INEXACT qreal aa3, bb3, cc3;
535 INEXACT qreal ti1, tj1;
536 qreal ti0, tj0;
537 qreal u[4], v[4];
538 INEXACT qreal u3, v3;
539 qreal temp8[8], temp16a[16], temp16b[16], temp16c[16];
540 qreal temp32a[32], temp32b[32], temp48[48], temp64[64];
541 int temp8len, temp16alen, temp16blen, temp16clen;
542 int temp32alen, temp32blen, temp48len, temp64len;
543 qreal axtbb[8], axtcc[8], aytbb[8], aytcc[8];
544 // cppcheck-suppress variableScope
545 int axtbblen, axtcclen, aytbblen, aytcclen;
546 qreal bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
547 // cppcheck-suppress variableScope
548 int bxtaalen, bxtcclen, bytaalen, bytcclen;
549 qreal cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
550 // cppcheck-suppress variableScope
551 int cxtaalen, cxtbblen, cytaalen, cytbblen;
552 qreal axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
553 int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
554 qreal axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
555 // cppcheck-suppress variableScope
556 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
557 qreal axtbctt[8], aytbctt[8], bxtcatt[8];
558 qreal bytcatt[8], cxtabtt[8], cytabtt[8];
559 // cppcheck-suppress variableScope
560 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
561 qreal abt[8], bct[8], cat[8];
562 // cppcheck-suppress variableScope
563 int abtlen, bctlen, catlen;
564 qreal abtt[4], bctt[4], catt[4];
565 // cppcheck-suppress variableScope
566 int abttlen, bcttlen, cattlen;
567 INEXACT qreal abtt3, bctt3, catt3;
568 qreal negate;
569
570 INEXACT qreal bvirt;
571 qreal avirt, bround, around;
572 INEXACT qreal c;
573 INEXACT qreal abig;
574 qreal ahi, alo, bhi, blo;
575 qreal err1, err2, err3;
576 INEXACT qreal _i, _j;
577 qreal _0;
578
579 adx = (qreal) (pa[0] - pd[0]);
580 bdx = (qreal) (pb[0] - pd[0]);
581 cdx = (qreal) (pc[0] - pd[0]);
582 ady = (qreal) (pa[1] - pd[1]);
583 bdy = (qreal) (pb[1] - pd[1]);
584 cdy = (qreal) (pc[1] - pd[1]);
585
586 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
587 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
588 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
589 bc[3] = bc3;
590 axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
591 axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
592 aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
593 ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
594 alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
595
596 Two_Product(cdx, ady, cdxady1, cdxady0);
597 Two_Product(adx, cdy, adxcdy1, adxcdy0);
598 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
599 ca[3] = ca3;
600 bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
601 bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
602 bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
603 byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
604 blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
605
606 Two_Product(adx, bdy, adxbdy1, adxbdy0);
607 Two_Product(bdx, ady, bdxady1, bdxady0);
608 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
609 ab[3] = ab3;
610 cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
611 cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
612 cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
613 cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
614 clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
615
616 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
617 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
618
619 det = estimate(finlength, fin1);
620 errbound = iccerrboundB * permanent;
621 if ((det >= errbound) || (-det >= errbound))
622 {
623 return det;
624 }
625
626 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
627 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
628 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
629 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
630 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
631 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
632 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
633 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0))
634 {
635 return det;
636 }
637
638 errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
639 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail))
640 + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
641 + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail))
642 + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
643 + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail))
644 + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
645 if ((det >= errbound) || (-det >= errbound))
646 {
647 return det;
648 }
649
650 finnow = fin1;
651 finother = fin2;
652
653 if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0))
654 {
655 Square(adx, adxadx1, adxadx0);
656 Square(ady, adyady1, adyady0);
657 Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
658 aa[3] = aa3;
659 }
660 if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0))
661 {
662 Square(bdx, bdxbdx1, bdxbdx0);
663 Square(bdy, bdybdy1, bdybdy0);
664 Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
665 bb[3] = bb3;
666 }
667 if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0))
668 {
669 Square(cdx, cdxcdx1, cdxcdx0);
670 Square(cdy, cdycdy1, cdycdy0);
671 Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
672 cc[3] = cc3;
673 }
674
675 if (adxtail != 0.0)
676 {
677 axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
678 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx, temp16a);
679
680 axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
681 temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
682
683 axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
684 temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
685
686 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
687 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
688 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
689 finswap = finnow; finnow = finother; finother = finswap;
690 }
691 if (adytail != 0.0)
692 {
693 aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
694 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady, temp16a);
695
696 aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
697 temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
698
699 aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
700 temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
701
702 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
703 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
704 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
705 finswap = finnow; finnow = finother; finother = finswap;
706 }
707 if (bdxtail != 0.0)
708 {
709 bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
710 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx, temp16a);
711
712 bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
713 temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
714
715 bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
716 temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
717
718 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
719 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
720 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
721 finswap = finnow; finnow = finother; finother = finswap;
722 }
723 if (bdytail != 0.0)
724 {
725 bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
726 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy, temp16a);
727
728 bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
729 temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
730
731 bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
732 temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
733
734 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
735 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
736 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
737 finswap = finnow; finnow = finother; finother = finswap;
738 }
739 if (cdxtail != 0.0)
740 {
741 cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
742 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx, temp16a);
743
744 cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
745 temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
746
747 cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
748 temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
749
750 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
751 temp16blen, temp16b, temp32a);
752 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
753 temp32alen, temp32a, temp48);
754 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
755 temp48, finother);
756 finswap = finnow; finnow = finother; finother = finswap;
757 }
758 if (cdytail != 0.0)
759 {
760 cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
761 temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy, temp16a);
762
763 cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
764 temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
765
766 cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
767 temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
768
769 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
770 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
771 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
772 finswap = finnow; finnow = finother; finother = finswap;
773 }
774
775 if ((adxtail != 0.0) || (adytail != 0.0))
776 {
777 if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0))
778 {
779 Two_Product(bdxtail, cdy, ti1, ti0);
780 Two_Product(bdx, cdytail, tj1, tj0);
781 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
782 u[3] = u3;
783 negate = -bdy;
784 Two_Product(cdxtail, negate, ti1, ti0);
785 negate = -bdytail;
786 Two_Product(cdx, negate, tj1, tj0);
787 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
788 v[3] = v3;
789 bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
790
791 Two_Product(bdxtail, cdytail, ti1, ti0);
792 Two_Product(cdxtail, bdytail, tj1, tj0);
793 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
794 bctt[3] = bctt3;
795 bcttlen = 4;
796 }
797 else
798 {
799 bct[0] = 0.0;
800 bctlen = 1;
801 bctt[0] = 0.0;
802 bcttlen = 1;
803 }
804
805 if (adxtail != 0.0)
806 {
807 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
808 axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
809 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx, temp32a);
810 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
811 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
812 finswap = finnow; finnow = finother; finother = finswap;
813 if (bdytail != 0.0)
814 {
815 temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
816 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a);
817 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
818 finswap = finnow; finnow = finother; finother = finswap;
819 }
820 if (cdytail != 0.0)
821 {
822 temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
823 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a);
824 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
825 finswap = finnow; finnow = finother; finother = finswap;
826 }
827
828 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail, temp32a);
829 axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
830 temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx, temp16a);
831 temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail, temp16b);
832 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
833 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
834 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
835 finswap = finnow; finnow = finother; finother = finswap;
836 }
837 if (adytail != 0.0)
838 {
839 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
840 aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
841 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady, temp32a);
842 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
843 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
844 finswap = finnow; finnow = finother; finother = finswap;
845
846
847 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail, temp32a);
848 aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
849 temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady, temp16a);
850 temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail, temp16b);
851 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
852 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
853 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
854 finswap = finnow; finnow = finother; finother = finswap;
855 }
856 }
857 if ((bdxtail != 0.0) || (bdytail != 0.0))
858 {
859 if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0))
860 {
861 Two_Product(cdxtail, ady, ti1, ti0);
862 Two_Product(cdx, adytail, tj1, tj0);
863 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
864 u[3] = u3;
865 negate = -cdy;
866 Two_Product(adxtail, negate, ti1, ti0);
867 negate = -cdytail;
868 Two_Product(adx, negate, tj1, tj0);
869 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
870 v[3] = v3;
871 catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
872
873 Two_Product(cdxtail, adytail, ti1, ti0);
874 Two_Product(adxtail, cdytail, tj1, tj0);
875 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
876 catt[3] = catt3;
877 cattlen = 4;
878 }
879 else
880 {
881 cat[0] = 0.0;
882 catlen = 1;
883 catt[0] = 0.0;
884 cattlen = 1;
885 }
886
887 if (bdxtail != 0.0)
888 {
889 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
890 bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
891 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx, temp32a);
892 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
893 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
894 finswap = finnow;
895 finnow = finother;
896 finother = finswap;
897 if (cdytail != 0.0)
898 {
899 temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
900 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a);
901 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
902 finswap = finnow;
903 finnow = finother;
904 finother = finswap;
905 }
906 if (adytail != 0.0)
907 {
908 temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
909 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a);
910 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
911 finswap = finnow;
912 finnow = finother;
913 finother = finswap;
914 }
915
916 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail, temp32a);
917 bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
918 temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx, temp16a);
919 temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail, temp16b);
920 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
921 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
922 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
923 finswap = finnow; finnow = finother; finother = finswap;
924 }
925 if (bdytail != 0.0)
926 {
927 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
928 bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
929 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy, temp32a);
930 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
931 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
932 finswap = finnow;
933 finnow = finother;
934 finother = finswap;
935
936
937 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail, temp32a);
938 bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
939 temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy, temp16a);
940 temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail, temp16b);
941 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
942 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
943 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
944 finswap = finnow;
945 finnow = finother;
946 finother = finswap;
947 }
948 }
949 if ((cdxtail != 0.0) || (cdytail != 0.0))
950 {
951 if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0))
952 {
953 Two_Product(adxtail, bdy, ti1, ti0);
954 Two_Product(adx, bdytail, tj1, tj0);
955 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
956 u[3] = u3;
957 negate = -ady;
958 Two_Product(bdxtail, negate, ti1, ti0);
959 negate = -adytail;
960 Two_Product(bdx, negate, tj1, tj0);
961 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
962 v[3] = v3;
963 abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
964
965 Two_Product(adxtail, bdytail, ti1, ti0);
966 Two_Product(bdxtail, adytail, tj1, tj0);
967 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
968 abtt[3] = abtt3;
969 abttlen = 4;
970 }
971 else
972 {
973 abt[0] = 0.0;
974 abtlen = 1;
975 abtt[0] = 0.0;
976 abttlen = 1;
977 }
978
979 if (cdxtail != 0.0)
980 {
981 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
982 cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
983 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx, temp32a);
984 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
985 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
986 finswap = finnow;
987 finnow = finother;
988 finother = finswap;
989 if (adytail != 0.0)
990 {
991 temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
992 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a);
993 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
994 finswap = finnow;
995 finnow = finother;
996 finother = finswap;
997 }
998 if (bdytail != 0.0)
999 {
1000 temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
1001 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a);
1002 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
1003 finswap = finnow; finnow = finother; finother = finswap;
1004 }
1005
1006 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail, temp32a);
1007 cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
1008 temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx, temp16a);
1009 temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail, temp16b);
1010 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1011 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1012 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1013 finswap = finnow;
1014 finnow = finother;
1015 finother = finswap;
1016 }
1017 if (cdytail != 0.0)
1018 {
1019 temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
1020 cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
1021 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy, temp32a);
1022 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
1023 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1024 finswap = finnow;
1025 finnow = finother;
1026 finother = finswap;
1027
1028
1029 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail, temp32a);
1030 cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
1031 temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy, temp16a);
1032 temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail, temp16b);
1033 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1034 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1035 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1036 finswap = finnow; finnow = finother; finother = finswap;
1037 }
1038 }
1039
1040 return finnow[finlength - 1];
1041 }
1042
1043
incircle(qreal * pa,qreal * pb,qreal * pc,qreal * pd)1044 qreal incircle(qreal *pa, qreal *pb, qreal *pc, qreal *pd)
1045 {
1046 qreal adx, bdx, cdx, ady, bdy, cdy;
1047 qreal bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
1048 qreal alift, blift, clift;
1049 qreal det;
1050 qreal permanent, errbound;
1051
1052 adx = pa[0] - pd[0];
1053 bdx = pb[0] - pd[0];
1054 cdx = pc[0] - pd[0];
1055 ady = pa[1] - pd[1];
1056 bdy = pb[1] - pd[1];
1057 cdy = pc[1] - pd[1];
1058
1059 bdxcdy = bdx * cdy;
1060 cdxbdy = cdx * bdy;
1061 alift = adx * adx + ady * ady;
1062
1063 cdxady = cdx * ady;
1064 adxcdy = adx * cdy;
1065 blift = bdx * bdx + bdy * bdy;
1066
1067 adxbdy = adx * bdy;
1068 bdxady = bdx * ady;
1069 clift = cdx * cdx + cdy * cdy;
1070
1071 det = alift * (bdxcdy - cdxbdy)
1072 + blift * (cdxady - adxcdy)
1073 + clift * (adxbdy - bdxady);
1074
1075 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
1076 + (Absolute(cdxady) + Absolute(adxcdy)) * blift
1077 + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
1078 errbound = iccerrboundA * permanent;
1079 if ((det > errbound) || (-det > errbound))
1080 {
1081 return det;
1082 }
1083
1084 return incircleadapt(pa, pb, pc, pd, permanent);
1085 }
1086
1087 QT_WARNING_POP
1088