00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 #include <stdio.h>
00117 #include <stdlib.h>
00118 #include <math.h>
00119
00120 #include "ApproxMVBB/GeometryPredicates/Predicates.hpp"
00121
00123 #include "ApproxMVBB/GeometryPredicates/Rounding.hpp"
00124
00125
00126 namespace GeometryPredicates{
00127
00129 FPU_DECLARE
00130
00131
00133 #define USE_PREDICATES_INIT
00134
00135
00136 #ifdef USE_PREDICATES_INIT
00137 #include "ApproxMVBB/GeometryPredicates/PredicatesInit.hpp"
00138 #endif
00139
00140
00141
00142
00143
00144
00145
00146
00147 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 #define Fast_Two_Sum_Tail(a, b, x, y) \
00164 bvirt = x - a; \
00165 y = b - bvirt
00166
00167 #define Fast_Two_Sum(a, b, x, y) \
00168 x = (REAL) (a + b); \
00169 Fast_Two_Sum_Tail(a, b, x, y)
00170
00171 #define Fast_Two_Diff_Tail(a, b, x, y) \
00172 bvirt = a - x; \
00173 y = bvirt - b
00174
00175 #define Fast_Two_Diff(a, b, x, y) \
00176 x = (REAL) (a - b); \
00177 Fast_Two_Diff_Tail(a, b, x, y)
00178
00179 #define Two_Sum_Tail(a, b, x, y) \
00180 bvirt = (REAL) (x - a); \
00181 avirt = x - bvirt; \
00182 bround = b - bvirt; \
00183 around = a - avirt; \
00184 y = around + bround
00185
00186 #define Two_Sum(a, b, x, y) \
00187 x = (REAL) (a + b); \
00188 Two_Sum_Tail(a, b, x, y)
00189
00190 #define Two_Diff_Tail(a, b, x, y) \
00191 bvirt = (REAL) (a - x); \
00192 avirt = x + bvirt; \
00193 bround = bvirt - b; \
00194 around = a - avirt; \
00195 y = around + bround
00196
00197 #define Two_Diff(a, b, x, y) \
00198 x = (REAL) (a - b); \
00199 Two_Diff_Tail(a, b, x, y)
00200
00201 #define Split(a, ahi, alo) \
00202 c = (REAL) (splitter * a); \
00203 abig = (REAL) (c - a); \
00204 ahi = c - abig; \
00205 alo = a - ahi
00206
00207 #define Two_Product_Tail(a, b, x, y) \
00208 Split(a, ahi, alo); \
00209 Split(b, bhi, blo); \
00210 err1 = x - (ahi * bhi); \
00211 err2 = err1 - (alo * bhi); \
00212 err3 = err2 - (ahi * blo); \
00213 y = (alo * blo) - err3
00214
00215 #define Two_Product(a, b, x, y) \
00216 x = (REAL) (a * b); \
00217 Two_Product_Tail(a, b, x, y)
00218
00219
00220
00221
00222 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
00223 x = (REAL) (a * b); \
00224 Split(a, ahi, alo); \
00225 err1 = x - (ahi * bhi); \
00226 err2 = err1 - (alo * bhi); \
00227 err3 = err2 - (ahi * blo); \
00228 y = (alo * blo) - err3
00229
00230
00231
00232
00233 #define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y) \
00234 x = (REAL) (a * b); \
00235 err1 = x - (ahi * bhi); \
00236 err2 = err1 - (alo * bhi); \
00237 err3 = err2 - (ahi * blo); \
00238 y = (alo * blo) - err3
00239
00240
00241
00242 #define Square_Tail(a, x, y) \
00243 Split(a, ahi, alo); \
00244 err1 = x - (ahi * ahi); \
00245 err3 = err1 - ((ahi + ahi) * alo); \
00246 y = (alo * alo) - err3
00247
00248 #define Square(a, x, y) \
00249 x = (REAL) (a * a); \
00250 Square_Tail(a, x, y)
00251
00252
00253
00254
00255 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
00256 Two_Sum(a0, b , _i, x0); \
00257 Two_Sum(a1, _i, x2, x1)
00258
00259 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
00260 Two_Diff(a0, b , _i, x0); \
00261 Two_Sum( a1, _i, x2, x1)
00262
00263 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
00264 Two_One_Sum(a1, a0, b0, _j, _0, x0); \
00265 Two_One_Sum(_j, _0, b1, x3, x2, x1)
00266
00267 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
00268 Two_One_Diff(a1, a0, b0, _j, _0, x0); \
00269 Two_One_Diff(_j, _0, b1, x3, x2, x1)
00270
00271 #define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \
00272 Two_One_Sum(a1, a0, b , _j, x1, x0); \
00273 Two_One_Sum(a3, a2, _j, x4, x3, x2)
00274
00275 #define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0) \
00276 Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0); \
00277 Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1)
00278
00279 #define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, \
00280 x1, x0) \
00281 Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0); \
00282 Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2)
00283
00284 #define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, \
00285 x3, x2, x1, x0) \
00286 Four_One_Sum(a3, a2, a1, a0, b , _j, x3, x2, x1, x0); \
00287 Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4)
00288
00289 #define Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, \
00290 x6, x5, x4, x3, x2, x1, x0) \
00291 Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, \
00292 _1, _0, x0); \
00293 Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, \
00294 x3, x2, x1)
00295
00296 #define Eight_Four_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b4, b3, b1, b0, x11, \
00297 x10, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0) \
00298 Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, \
00299 _2, _1, _0, x1, x0); \
00300 Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, \
00301 x7, x6, x5, x4, x3, x2)
00302
00303
00304
00305 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
00306 Split(b, bhi, blo); \
00307 Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
00308 Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
00309 Two_Sum(_i, _0, _k, x1); \
00310 Fast_Two_Sum(_j, _k, x3, x2)
00311
00312 #define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \
00313 Split(b, bhi, blo); \
00314 Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
00315 Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
00316 Two_Sum(_i, _0, _k, x1); \
00317 Fast_Two_Sum(_j, _k, _i, x2); \
00318 Two_Product_Presplit(a2, b, bhi, blo, _j, _0); \
00319 Two_Sum(_i, _0, _k, x3); \
00320 Fast_Two_Sum(_j, _k, _i, x4); \
00321 Two_Product_Presplit(a3, b, bhi, blo, _j, _0); \
00322 Two_Sum(_i, _0, _k, x5); \
00323 Fast_Two_Sum(_j, _k, x7, x6)
00324
00325 #define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \
00326 Split(a0, a0hi, a0lo); \
00327 Split(b0, bhi, blo); \
00328 Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0); \
00329 Split(a1, a1hi, a1lo); \
00330 Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0); \
00331 Two_Sum(_i, _0, _k, _1); \
00332 Fast_Two_Sum(_j, _k, _l, _2); \
00333 Split(b1, bhi, blo); \
00334 Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0); \
00335 Two_Sum(_1, _0, _k, x1); \
00336 Two_Sum(_2, _k, _j, _1); \
00337 Two_Sum(_l, _j, _m, _2); \
00338 Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0); \
00339 Two_Sum(_i, _0, _n, _0); \
00340 Two_Sum(_1, _0, _i, x2); \
00341 Two_Sum(_2, _i, _k, _1); \
00342 Two_Sum(_m, _k, _l, _2); \
00343 Two_Sum(_j, _n, _k, _0); \
00344 Two_Sum(_1, _0, _j, x3); \
00345 Two_Sum(_2, _j, _i, _1); \
00346 Two_Sum(_l, _i, _m, _2); \
00347 Two_Sum(_1, _k, _i, x4); \
00348 Two_Sum(_2, _i, _k, x5); \
00349 Two_Sum(_m, _k, x7, x6)
00350
00351
00352
00353
00354
00355 #define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \
00356 Square(a0, _j, x0); \
00357 _0 = a0 + a0; \
00358 Two_Product(a1, _0, _k, _1); \
00359 Two_One_Sum(_k, _1, _j, _l, _2, x1); \
00360 Square(a1, _j, _1); \
00361 Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2)
00362
00363 #ifndef USE_PREDICATES_INIT
00364
00365 static REAL splitter;
00366
00367 static REAL resulterrbound;
00368 static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
00369 static REAL o3derrboundA, o3derrboundB, o3derrboundC;
00370 static REAL iccerrboundA, iccerrboundB, iccerrboundC;
00371 static REAL isperrboundA, isperrboundB, isperrboundC;
00372
00373 #endif
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655 static int fast_expansion_sum_zeroelim(int elen, REAL *e,
00656 int flen, REAL *f, REAL *h)
00657
00658 {
00659 REAL Q;
00660 INEXACT REAL Qnew;
00661 INEXACT REAL hh;
00662 INEXACT REAL bvirt;
00663 REAL avirt, bround, around;
00664 int eindex, findex, hindex;
00665 REAL enow, fnow;
00666
00667 enow = e[0];
00668 fnow = f[0];
00669 eindex = findex = 0;
00670 if ((fnow > enow) == (fnow > -enow)) {
00671 Q = enow;
00672 enow = e[++eindex];
00673 } else {
00674 Q = fnow;
00675 fnow = f[++findex];
00676 }
00677 hindex = 0;
00678 if ((eindex < elen-1) && (findex < flen-1)) {
00679 if ((fnow > enow) == (fnow > -enow)) {
00680 Fast_Two_Sum(enow, Q, Qnew, hh);
00681 enow = e[++eindex];
00682 } else {
00683 Fast_Two_Sum(fnow, Q, Qnew, hh);
00684 fnow = f[++findex];
00685 }
00686 Q = Qnew;
00687 if (hh != 0.0) {
00688 h[hindex++] = hh;
00689 }
00690 while ((eindex < elen-1) && (findex < flen-1)) {
00691 if ((fnow > enow) == (fnow > -enow)) {
00692 Two_Sum(Q, enow, Qnew, hh);
00693 enow = e[++eindex];
00694 } else {
00695 Two_Sum(Q, fnow, Qnew, hh);
00696 fnow = f[++findex];
00697 }
00698 Q = Qnew;
00699 if (hh != 0.0) {
00700 h[hindex++] = hh;
00701 }
00702 }
00703 }
00704 while (eindex < elen-1) {
00705 Two_Sum(Q, enow, Qnew, hh);
00706 enow = e[++eindex];
00707 Q = Qnew;
00708 if (hh != 0.0) {
00709 h[hindex++] = hh;
00710 }
00711 }
00712 while (findex < flen-1) {
00713 Two_Sum(Q, fnow, Qnew, hh);
00714 fnow = f[++findex];
00715 Q = Qnew;
00716 if (hh != 0.0) {
00717 h[hindex++] = hh;
00718 }
00719 }
00720 if ((Q != 0.0) || (hindex == 0)) {
00721 h[hindex++] = Q;
00722 }
00723 return hindex;
00724 }
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741 static int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
00742
00743 {
00744 INEXACT REAL Q, sum;
00745 REAL hh;
00746 INEXACT REAL product1;
00747 REAL product0;
00748 int eindex, hindex;
00749 REAL enow;
00750 INEXACT REAL bvirt;
00751 REAL avirt, bround, around;
00752 INEXACT REAL c;
00753 INEXACT REAL abig;
00754 REAL ahi, alo, bhi, blo;
00755 REAL err1, err2, err3;
00756
00757 Split(b, bhi, blo);
00758 Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
00759 hindex = 0;
00760 if (hh != 0) {
00761 h[hindex++] = hh;
00762 }
00763 for (eindex = 1; eindex < elen; eindex++) {
00764 enow = e[eindex];
00765 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
00766 Two_Sum(Q, product0, sum, hh);
00767 if (hh != 0) {
00768 h[hindex++] = hh;
00769 }
00770 Fast_Two_Sum(product1, sum, Q, hh);
00771 if (hh != 0) {
00772 h[hindex++] = hh;
00773 }
00774 }
00775 if ((Q != 0.0) || (hindex == 0)) {
00776 h[hindex++] = Q;
00777 }
00778 return hindex;
00779 }
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789 static REAL estimate(int elen, REAL *e)
00790 {
00791 REAL Q;
00792 int eindex;
00793
00794 Q = e[0];
00795 for (eindex = 1; eindex < elen; eindex++) {
00796 Q += e[eindex];
00797 }
00798 return Q;
00799 }
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827 static REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum)
00828 {
00829 INEXACT REAL acx, acy, bcx, bcy;
00830 REAL acxtail, acytail, bcxtail, bcytail;
00831 INEXACT REAL detleft, detright;
00832 REAL detlefttail, detrighttail;
00833 REAL det, errbound;
00834 REAL B[4], C1[8], C2[12], D[16];
00835 INEXACT REAL B3;
00836 int C1length, C2length, Dlength;
00837 REAL u[4];
00838 INEXACT REAL u3;
00839 INEXACT REAL s1, t1;
00840 REAL s0, t0;
00841
00842 INEXACT REAL bvirt;
00843 REAL avirt, bround, around;
00844 INEXACT REAL c;
00845 INEXACT REAL abig;
00846 REAL ahi, alo, bhi, blo;
00847 REAL err1, err2, err3;
00848 INEXACT REAL _i, _j;
00849 REAL _0;
00850
00851 acx = (REAL) (pa[0] - pc[0]);
00852 bcx = (REAL) (pb[0] - pc[0]);
00853 acy = (REAL) (pa[1] - pc[1]);
00854 bcy = (REAL) (pb[1] - pc[1]);
00855
00856 Two_Product(acx, bcy, detleft, detlefttail);
00857 Two_Product(acy, bcx, detright, detrighttail);
00858
00859 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
00860 B3, B[2], B[1], B[0]);
00861 B[3] = B3;
00862
00863 det = estimate(4, B);
00864 errbound = ccwerrboundB * detsum;
00865 if ((det >= errbound) || (-det >= errbound)) {
00866 return det;
00867 }
00868
00869 Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
00870 Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
00871 Two_Diff_Tail(pa[1], pc[1], acy, acytail);
00872 Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
00873
00874 if ((acxtail == 0.0) && (acytail == 0.0)
00875 && (bcxtail == 0.0) && (bcytail == 0.0)) {
00876 return det;
00877 }
00878
00879 errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
00880 det += (acx * bcytail + bcy * acxtail)
00881 - (acy * bcxtail + bcx * acytail);
00882 if ((det >= errbound) || (-det >= errbound)) {
00883 return det;
00884 }
00885
00886 Two_Product(acxtail, bcy, s1, s0);
00887 Two_Product(acytail, bcx, t1, t0);
00888 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
00889 u[3] = u3;
00890 C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
00891
00892 Two_Product(acx, bcytail, s1, s0);
00893 Two_Product(acy, bcxtail, t1, t0);
00894 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
00895 u[3] = u3;
00896 C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
00897
00898 Two_Product(acxtail, bcytail, s1, s0);
00899 Two_Product(acytail, bcxtail, t1, t0);
00900 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
00901 u[3] = u3;
00902 Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
00903
00904 return(D[Dlength - 1]);
00905 }
00906
00907 REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
00908 {
00909 REAL detleft, detright, det;
00910 REAL detsum, errbound;
00911 REAL orient;
00912
00913 FPU_ROUND_DOUBLE;
00914
00915 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
00916 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
00917 det = detleft - detright;
00918
00919 if (detleft > 0.0) {
00920 if (detright <= 0.0) {
00921 FPU_RESTORE;
00922 return det;
00923 } else {
00924 detsum = detleft + detright;
00925 }
00926 } else if (detleft < 0.0) {
00927 if (detright >= 0.0) {
00928 FPU_RESTORE;
00929 return det;
00930 } else {
00931 detsum = -detleft - detright;
00932 }
00933 } else {
00934 FPU_RESTORE;
00935 return det;
00936 }
00937
00938 errbound = ccwerrboundA * detsum;
00939 if ((det >= errbound) || (-det >= errbound)) {
00940 FPU_RESTORE;
00941 return det;
00942 }
00943
00944 orient = orient2dadapt(pa, pb, pc, detsum);
00945 FPU_RESTORE;
00946 return orient;
00947 }
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978 static REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd,
00979 REAL permanent)
00980 {
00981 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
00982 REAL det, errbound;
00983
00984 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
00985 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
00986 REAL bc[4], ca[4], ab[4];
00987 INEXACT REAL bc3, ca3, ab3;
00988 REAL adet[8], bdet[8], cdet[8];
00989 int alen, blen, clen;
00990 REAL abdet[16];
00991 int ablen;
00992 REAL *finnow, *finother, *finswap;
00993 REAL fin1[192], fin2[192];
00994 int finlength;
00995
00996 REAL adxtail, bdxtail, cdxtail;
00997 REAL adytail, bdytail, cdytail;
00998 REAL adztail, bdztail, cdztail;
00999 INEXACT REAL at_blarge, at_clarge;
01000 INEXACT REAL bt_clarge, bt_alarge;
01001 INEXACT REAL ct_alarge, ct_blarge;
01002 REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
01003 int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
01004 INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
01005 INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1;
01006 REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
01007 REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
01008 INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
01009 INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1;
01010 REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
01011 REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
01012 REAL bct[8], cat[8], abt[8];
01013 int bctlen, catlen, abtlen;
01014 INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
01015 INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
01016 REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
01017 REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
01018 REAL u[4], v[12], w[16];
01019 INEXACT REAL u3;
01020 int vlength, wlength;
01021 REAL negate;
01022
01023 INEXACT REAL bvirt;
01024 REAL avirt, bround, around;
01025 INEXACT REAL c;
01026 INEXACT REAL abig;
01027 REAL ahi, alo, bhi, blo;
01028 REAL err1, err2, err3;
01029 INEXACT REAL _i, _j, _k;
01030 REAL _0;
01031
01032 adx = (REAL) (pa[0] - pd[0]);
01033 bdx = (REAL) (pb[0] - pd[0]);
01034 cdx = (REAL) (pc[0] - pd[0]);
01035 ady = (REAL) (pa[1] - pd[1]);
01036 bdy = (REAL) (pb[1] - pd[1]);
01037 cdy = (REAL) (pc[1] - pd[1]);
01038 adz = (REAL) (pa[2] - pd[2]);
01039 bdz = (REAL) (pb[2] - pd[2]);
01040 cdz = (REAL) (pc[2] - pd[2]);
01041
01042 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
01043 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
01044 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
01045 bc[3] = bc3;
01046 alen = scale_expansion_zeroelim(4, bc, adz, adet);
01047
01048 Two_Product(cdx, ady, cdxady1, cdxady0);
01049 Two_Product(adx, cdy, adxcdy1, adxcdy0);
01050 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
01051 ca[3] = ca3;
01052 blen = scale_expansion_zeroelim(4, ca, bdz, bdet);
01053
01054 Two_Product(adx, bdy, adxbdy1, adxbdy0);
01055 Two_Product(bdx, ady, bdxady1, bdxady0);
01056 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
01057 ab[3] = ab3;
01058 clen = scale_expansion_zeroelim(4, ab, cdz, cdet);
01059
01060 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
01061 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
01062
01063 det = estimate(finlength, fin1);
01064 errbound = o3derrboundB * permanent;
01065 if ((det >= errbound) || (-det >= errbound)) {
01066 return det;
01067 }
01068
01069 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
01070 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
01071 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
01072 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
01073 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
01074 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
01075 Two_Diff_Tail(pa[2], pd[2], adz, adztail);
01076 Two_Diff_Tail(pb[2], pd[2], bdz, bdztail);
01077 Two_Diff_Tail(pc[2], pd[2], cdz, cdztail);
01078
01079 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
01080 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)
01081 && (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) {
01082 return det;
01083 }
01084
01085 errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
01086 det += (adz * ((bdx * cdytail + cdy * bdxtail)
01087 - (bdy * cdxtail + cdx * bdytail))
01088 + adztail * (bdx * cdy - bdy * cdx))
01089 + (bdz * ((cdx * adytail + ady * cdxtail)
01090 - (cdy * adxtail + adx * cdytail))
01091 + bdztail * (cdx * ady - cdy * adx))
01092 + (cdz * ((adx * bdytail + bdy * adxtail)
01093 - (ady * bdxtail + bdx * adytail))
01094 + cdztail * (adx * bdy - ady * bdx));
01095 if ((det >= errbound) || (-det >= errbound)) {
01096 return det;
01097 }
01098
01099 finnow = fin1;
01100 finother = fin2;
01101
01102 if (adxtail == 0.0) {
01103 if (adytail == 0.0) {
01104 at_b[0] = 0.0;
01105 at_blen = 1;
01106 at_c[0] = 0.0;
01107 at_clen = 1;
01108 } else {
01109 negate = -adytail;
01110 Two_Product(negate, bdx, at_blarge, at_b[0]);
01111 at_b[1] = at_blarge;
01112 at_blen = 2;
01113 Two_Product(adytail, cdx, at_clarge, at_c[0]);
01114 at_c[1] = at_clarge;
01115 at_clen = 2;
01116 }
01117 } else {
01118 if (adytail == 0.0) {
01119 Two_Product(adxtail, bdy, at_blarge, at_b[0]);
01120 at_b[1] = at_blarge;
01121 at_blen = 2;
01122 negate = -adxtail;
01123 Two_Product(negate, cdy, at_clarge, at_c[0]);
01124 at_c[1] = at_clarge;
01125 at_clen = 2;
01126 } else {
01127 Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
01128 Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
01129 Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
01130 at_blarge, at_b[2], at_b[1], at_b[0]);
01131 at_b[3] = at_blarge;
01132 at_blen = 4;
01133 Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
01134 Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
01135 Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
01136 at_clarge, at_c[2], at_c[1], at_c[0]);
01137 at_c[3] = at_clarge;
01138 at_clen = 4;
01139 }
01140 }
01141 if (bdxtail == 0.0) {
01142 if (bdytail == 0.0) {
01143 bt_c[0] = 0.0;
01144 bt_clen = 1;
01145 bt_a[0] = 0.0;
01146 bt_alen = 1;
01147 } else {
01148 negate = -bdytail;
01149 Two_Product(negate, cdx, bt_clarge, bt_c[0]);
01150 bt_c[1] = bt_clarge;
01151 bt_clen = 2;
01152 Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
01153 bt_a[1] = bt_alarge;
01154 bt_alen = 2;
01155 }
01156 } else {
01157 if (bdytail == 0.0) {
01158 Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
01159 bt_c[1] = bt_clarge;
01160 bt_clen = 2;
01161 negate = -bdxtail;
01162 Two_Product(negate, ady, bt_alarge, bt_a[0]);
01163 bt_a[1] = bt_alarge;
01164 bt_alen = 2;
01165 } else {
01166 Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
01167 Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
01168 Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
01169 bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
01170 bt_c[3] = bt_clarge;
01171 bt_clen = 4;
01172 Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
01173 Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
01174 Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
01175 bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
01176 bt_a[3] = bt_alarge;
01177 bt_alen = 4;
01178 }
01179 }
01180 if (cdxtail == 0.0) {
01181 if (cdytail == 0.0) {
01182 ct_a[0] = 0.0;
01183 ct_alen = 1;
01184 ct_b[0] = 0.0;
01185 ct_blen = 1;
01186 } else {
01187 negate = -cdytail;
01188 Two_Product(negate, adx, ct_alarge, ct_a[0]);
01189 ct_a[1] = ct_alarge;
01190 ct_alen = 2;
01191 Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
01192 ct_b[1] = ct_blarge;
01193 ct_blen = 2;
01194 }
01195 } else {
01196 if (cdytail == 0.0) {
01197 Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
01198 ct_a[1] = ct_alarge;
01199 ct_alen = 2;
01200 negate = -cdxtail;
01201 Two_Product(negate, bdy, ct_blarge, ct_b[0]);
01202 ct_b[1] = ct_blarge;
01203 ct_blen = 2;
01204 } else {
01205 Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
01206 Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
01207 Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
01208 ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
01209 ct_a[3] = ct_alarge;
01210 ct_alen = 4;
01211 Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
01212 Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
01213 Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
01214 ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
01215 ct_b[3] = ct_blarge;
01216 ct_blen = 4;
01217 }
01218 }
01219
01220 bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
01221 wlength = scale_expansion_zeroelim(bctlen, bct, adz, w);
01222 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
01223 finother);
01224 finswap = finnow; finnow = finother; finother = finswap;
01225
01226 catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
01227 wlength = scale_expansion_zeroelim(catlen, cat, bdz, w);
01228 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
01229 finother);
01230 finswap = finnow; finnow = finother; finother = finswap;
01231
01232 abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
01233 wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w);
01234 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
01235 finother);
01236 finswap = finnow; finnow = finother; finother = finswap;
01237
01238 if (adztail != 0.0) {
01239 vlength = scale_expansion_zeroelim(4, bc, adztail, v);
01240 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
01241 finother);
01242 finswap = finnow; finnow = finother; finother = finswap;
01243 }
01244 if (bdztail != 0.0) {
01245 vlength = scale_expansion_zeroelim(4, ca, bdztail, v);
01246 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
01247 finother);
01248 finswap = finnow; finnow = finother; finother = finswap;
01249 }
01250 if (cdztail != 0.0) {
01251 vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
01252 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
01253 finother);
01254 finswap = finnow; finnow = finother; finother = finswap;
01255 }
01256
01257 if (adxtail != 0.0) {
01258 if (bdytail != 0.0) {
01259 Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
01260 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]);
01261 u[3] = u3;
01262 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01263 finother);
01264 finswap = finnow; finnow = finother; finother = finswap;
01265 if (cdztail != 0.0) {
01266 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]);
01267 u[3] = u3;
01268 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01269 finother);
01270 finswap = finnow; finnow = finother; finother = finswap;
01271 }
01272 }
01273 if (cdytail != 0.0) {
01274 negate = -adxtail;
01275 Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
01276 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]);
01277 u[3] = u3;
01278 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01279 finother);
01280 finswap = finnow; finnow = finother; finother = finswap;
01281 if (bdztail != 0.0) {
01282 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]);
01283 u[3] = u3;
01284 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01285 finother);
01286 finswap = finnow; finnow = finother; finother = finswap;
01287 }
01288 }
01289 }
01290 if (bdxtail != 0.0) {
01291 if (cdytail != 0.0) {
01292 Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
01293 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]);
01294 u[3] = u3;
01295 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01296 finother);
01297 finswap = finnow; finnow = finother; finother = finswap;
01298 if (adztail != 0.0) {
01299 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]);
01300 u[3] = u3;
01301 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01302 finother);
01303 finswap = finnow; finnow = finother; finother = finswap;
01304 }
01305 }
01306 if (adytail != 0.0) {
01307 negate = -bdxtail;
01308 Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
01309 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]);
01310 u[3] = u3;
01311 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01312 finother);
01313 finswap = finnow; finnow = finother; finother = finswap;
01314 if (cdztail != 0.0) {
01315 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]);
01316 u[3] = u3;
01317 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01318 finother);
01319 finswap = finnow; finnow = finother; finother = finswap;
01320 }
01321 }
01322 }
01323 if (cdxtail != 0.0) {
01324 if (adytail != 0.0) {
01325 Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
01326 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]);
01327 u[3] = u3;
01328 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01329 finother);
01330 finswap = finnow; finnow = finother; finother = finswap;
01331 if (bdztail != 0.0) {
01332 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]);
01333 u[3] = u3;
01334 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01335 finother);
01336 finswap = finnow; finnow = finother; finother = finswap;
01337 }
01338 }
01339 if (bdytail != 0.0) {
01340 negate = -cdxtail;
01341 Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
01342 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]);
01343 u[3] = u3;
01344 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01345 finother);
01346 finswap = finnow; finnow = finother; finother = finswap;
01347 if (adztail != 0.0) {
01348 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]);
01349 u[3] = u3;
01350 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
01351 finother);
01352 finswap = finnow; finnow = finother; finother = finswap;
01353 }
01354 }
01355 }
01356
01357 if (adztail != 0.0) {
01358 wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w);
01359 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
01360 finother);
01361 finswap = finnow; finnow = finother; finother = finswap;
01362 }
01363 if (bdztail != 0.0) {
01364 wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w);
01365 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
01366 finother);
01367 finswap = finnow; finnow = finother; finother = finswap;
01368 }
01369 if (cdztail != 0.0) {
01370 wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
01371 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
01372 finother);
01373 finswap = finnow; finnow = finother; finother = finswap;
01374 }
01375
01376 return finnow[finlength - 1];
01377 }
01378
01379 REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
01380 {
01381 REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
01382 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
01383 REAL det;
01384 REAL permanent, errbound;
01385 REAL orient;
01386
01387 FPU_ROUND_DOUBLE;
01388
01389 adx = pa[0] - pd[0];
01390 bdx = pb[0] - pd[0];
01391 cdx = pc[0] - pd[0];
01392 ady = pa[1] - pd[1];
01393 bdy = pb[1] - pd[1];
01394 cdy = pc[1] - pd[1];
01395 adz = pa[2] - pd[2];
01396 bdz = pb[2] - pd[2];
01397 cdz = pc[2] - pd[2];
01398
01399 bdxcdy = bdx * cdy;
01400 cdxbdy = cdx * bdy;
01401
01402 cdxady = cdx * ady;
01403 adxcdy = adx * cdy;
01404
01405 adxbdy = adx * bdy;
01406 bdxady = bdx * ady;
01407
01408 det = adz * (bdxcdy - cdxbdy)
01409 + bdz * (cdxady - adxcdy)
01410 + cdz * (adxbdy - bdxady);
01411
01412 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
01413 + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
01414 + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
01415 errbound = o3derrboundA * permanent;
01416 if ((det > errbound) || (-det > errbound)) {
01417 FPU_RESTORE;
01418 return det;
01419 }
01420
01421 orient = orient3dadapt(pa, pb, pc, pd, permanent);
01422 FPU_RESTORE;
01423 return orient;
01424 }
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452 static REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd,
01453 REAL permanent)
01454 {
01455 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
01456 REAL det, errbound;
01457
01458 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
01459 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
01460 REAL bc[4], ca[4], ab[4];
01461 INEXACT REAL bc3, ca3, ab3;
01462 REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
01463 int axbclen, axxbclen, aybclen, ayybclen, alen;
01464 REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
01465 int bxcalen, bxxcalen, bycalen, byycalen, blen;
01466 REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
01467 int cxablen, cxxablen, cyablen, cyyablen, clen;
01468 REAL abdet[64];
01469 int ablen;
01470 REAL fin1[1152], fin2[1152];
01471 REAL *finnow, *finother, *finswap;
01472 int finlength;
01473
01474 REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
01475 INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
01476 REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
01477 REAL aa[4], bb[4], cc[4];
01478 INEXACT REAL aa3, bb3, cc3;
01479 INEXACT REAL ti1, tj1;
01480 REAL ti0, tj0;
01481 REAL u[4], v[4];
01482 INEXACT REAL u3, v3;
01483 REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
01484 REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
01485 int temp8len, temp16alen, temp16blen, temp16clen;
01486 int temp32alen, temp32blen, temp48len, temp64len;
01487 REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
01488 int axtbblen, axtcclen, aytbblen, aytcclen;
01489 REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
01490 int bxtaalen, bxtcclen, bytaalen, bytcclen;
01491 REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
01492 int cxtaalen, cxtbblen, cytaalen, cytbblen;
01493 REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
01494 int axtbclen = 0, aytbclen = 0;
01495 int bxtcalen = 0, bytcalen = 0;
01496 int cxtablen = 0, cytablen = 0;
01497 REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
01498 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
01499 REAL axtbctt[8], aytbctt[8], bxtcatt[8];
01500 REAL bytcatt[8], cxtabtt[8], cytabtt[8];
01501 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
01502 REAL abt[8], bct[8], cat[8];
01503 int abtlen, bctlen, catlen;
01504 REAL abtt[4], bctt[4], catt[4];
01505 int abttlen, bcttlen, cattlen;
01506 INEXACT REAL abtt3, bctt3, catt3;
01507 REAL negate;
01508
01509 INEXACT REAL bvirt;
01510 REAL avirt, bround, around;
01511 INEXACT REAL c;
01512 INEXACT REAL abig;
01513 REAL ahi, alo, bhi, blo;
01514 REAL err1, err2, err3;
01515 INEXACT REAL _i, _j;
01516 REAL _0;
01517
01518 adx = (REAL) (pa[0] - pd[0]);
01519 bdx = (REAL) (pb[0] - pd[0]);
01520 cdx = (REAL) (pc[0] - pd[0]);
01521 ady = (REAL) (pa[1] - pd[1]);
01522 bdy = (REAL) (pb[1] - pd[1]);
01523 cdy = (REAL) (pc[1] - pd[1]);
01524
01525 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
01526 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
01527 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
01528 bc[3] = bc3;
01529 axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
01530 axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
01531 aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
01532 ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
01533 alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
01534
01535 Two_Product(cdx, ady, cdxady1, cdxady0);
01536 Two_Product(adx, cdy, adxcdy1, adxcdy0);
01537 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
01538 ca[3] = ca3;
01539 bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
01540 bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
01541 bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
01542 byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
01543 blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
01544
01545 Two_Product(adx, bdy, adxbdy1, adxbdy0);
01546 Two_Product(bdx, ady, bdxady1, bdxady0);
01547 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
01548 ab[3] = ab3;
01549 cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
01550 cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
01551 cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
01552 cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
01553 clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
01554
01555 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
01556 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
01557
01558 det = estimate(finlength, fin1);
01559 errbound = iccerrboundB * permanent;
01560 if ((det >= errbound) || (-det >= errbound)) {
01561 return det;
01562 }
01563
01564 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
01565 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
01566 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
01567 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
01568 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
01569 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
01570 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
01571 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
01572 return det;
01573 }
01574
01575 errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
01576 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
01577 - (bdy * cdxtail + cdx * bdytail))
01578 + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
01579 + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
01580 - (cdy * adxtail + adx * cdytail))
01581 + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
01582 + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
01583 - (ady * bdxtail + bdx * adytail))
01584 + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
01585 if ((det >= errbound) || (-det >= errbound)) {
01586 return det;
01587 }
01588
01589 finnow = fin1;
01590 finother = fin2;
01591
01592 if ((bdxtail != 0.0) || (bdytail != 0.0)
01593 || (cdxtail != 0.0) || (cdytail != 0.0)) {
01594 Square(adx, adxadx1, adxadx0);
01595 Square(ady, adyady1, adyady0);
01596 Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
01597 aa[3] = aa3;
01598 }
01599 if ((cdxtail != 0.0) || (cdytail != 0.0)
01600 || (adxtail != 0.0) || (adytail != 0.0)) {
01601 Square(bdx, bdxbdx1, bdxbdx0);
01602 Square(bdy, bdybdy1, bdybdy0);
01603 Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
01604 bb[3] = bb3;
01605 }
01606 if ((adxtail != 0.0) || (adytail != 0.0)
01607 || (bdxtail != 0.0) || (bdytail != 0.0)) {
01608 Square(cdx, cdxcdx1, cdxcdx0);
01609 Square(cdy, cdycdy1, cdycdy0);
01610 Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
01611 cc[3] = cc3;
01612 }
01613
01614 if (adxtail != 0.0) {
01615 axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
01616 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
01617 temp16a);
01618
01619 axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
01620 temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
01621
01622 axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
01623 temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
01624
01625 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01626 temp16blen, temp16b, temp32a);
01627 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
01628 temp32alen, temp32a, temp48);
01629 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01630 temp48, finother);
01631 finswap = finnow; finnow = finother; finother = finswap;
01632 }
01633 if (adytail != 0.0) {
01634 aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
01635 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
01636 temp16a);
01637
01638 aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
01639 temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
01640
01641 aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
01642 temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
01643
01644 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01645 temp16blen, temp16b, temp32a);
01646 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
01647 temp32alen, temp32a, temp48);
01648 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01649 temp48, finother);
01650 finswap = finnow; finnow = finother; finother = finswap;
01651 }
01652 if (bdxtail != 0.0) {
01653 bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
01654 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
01655 temp16a);
01656
01657 bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
01658 temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
01659
01660 bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
01661 temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
01662
01663 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01664 temp16blen, temp16b, temp32a);
01665 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
01666 temp32alen, temp32a, temp48);
01667 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01668 temp48, finother);
01669 finswap = finnow; finnow = finother; finother = finswap;
01670 }
01671 if (bdytail != 0.0) {
01672 bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
01673 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
01674 temp16a);
01675
01676 bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
01677 temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
01678
01679 bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
01680 temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
01681
01682 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01683 temp16blen, temp16b, temp32a);
01684 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
01685 temp32alen, temp32a, temp48);
01686 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01687 temp48, finother);
01688 finswap = finnow; finnow = finother; finother = finswap;
01689 }
01690 if (cdxtail != 0.0) {
01691 cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
01692 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
01693 temp16a);
01694
01695 cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
01696 temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
01697
01698 cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
01699 temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
01700
01701 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01702 temp16blen, temp16b, temp32a);
01703 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
01704 temp32alen, temp32a, temp48);
01705 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01706 temp48, finother);
01707 finswap = finnow; finnow = finother; finother = finswap;
01708 }
01709 if (cdytail != 0.0) {
01710 cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
01711 temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
01712 temp16a);
01713
01714 cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
01715 temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
01716
01717 cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
01718 temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
01719
01720 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01721 temp16blen, temp16b, temp32a);
01722 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
01723 temp32alen, temp32a, temp48);
01724 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01725 temp48, finother);
01726 finswap = finnow; finnow = finother; finother = finswap;
01727 }
01728
01729 if ((adxtail != 0.0) || (adytail != 0.0)) {
01730 if ((bdxtail != 0.0) || (bdytail != 0.0)
01731 || (cdxtail != 0.0) || (cdytail != 0.0)) {
01732 Two_Product(bdxtail, cdy, ti1, ti0);
01733 Two_Product(bdx, cdytail, tj1, tj0);
01734 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
01735 u[3] = u3;
01736 negate = -bdy;
01737 Two_Product(cdxtail, negate, ti1, ti0);
01738 negate = -bdytail;
01739 Two_Product(cdx, negate, tj1, tj0);
01740 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
01741 v[3] = v3;
01742 bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
01743
01744 Two_Product(bdxtail, cdytail, ti1, ti0);
01745 Two_Product(cdxtail, bdytail, tj1, tj0);
01746 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
01747 bctt[3] = bctt3;
01748 bcttlen = 4;
01749 } else {
01750 bct[0] = 0.0;
01751 bctlen = 1;
01752 bctt[0] = 0.0;
01753 bcttlen = 1;
01754 }
01755
01756 if (adxtail != 0.0) {
01757 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
01758 axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
01759 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
01760 temp32a);
01761 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01762 temp32alen, temp32a, temp48);
01763 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01764 temp48, finother);
01765 finswap = finnow; finnow = finother; finother = finswap;
01766 if (bdytail != 0.0) {
01767 temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
01768 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
01769 temp16a);
01770 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
01771 temp16a, finother);
01772 finswap = finnow; finnow = finother; finother = finswap;
01773 }
01774 if (cdytail != 0.0) {
01775 temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
01776 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
01777 temp16a);
01778 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
01779 temp16a, finother);
01780 finswap = finnow; finnow = finother; finother = finswap;
01781 }
01782
01783 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
01784 temp32a);
01785 axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
01786 temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
01787 temp16a);
01788 temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
01789 temp16b);
01790 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01791 temp16blen, temp16b, temp32b);
01792 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
01793 temp32blen, temp32b, temp64);
01794 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
01795 temp64, finother);
01796 finswap = finnow; finnow = finother; finother = finswap;
01797 }
01798 if (adytail != 0.0) {
01799 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
01800 aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
01801 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
01802 temp32a);
01803 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01804 temp32alen, temp32a, temp48);
01805 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01806 temp48, finother);
01807 finswap = finnow; finnow = finother; finother = finswap;
01808
01809
01810 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
01811 temp32a);
01812 aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
01813 temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
01814 temp16a);
01815 temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
01816 temp16b);
01817 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01818 temp16blen, temp16b, temp32b);
01819 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
01820 temp32blen, temp32b, temp64);
01821 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
01822 temp64, finother);
01823 finswap = finnow; finnow = finother; finother = finswap;
01824 }
01825 }
01826 if ((bdxtail != 0.0) || (bdytail != 0.0)) {
01827 if ((cdxtail != 0.0) || (cdytail != 0.0)
01828 || (adxtail != 0.0) || (adytail != 0.0)) {
01829 Two_Product(cdxtail, ady, ti1, ti0);
01830 Two_Product(cdx, adytail, tj1, tj0);
01831 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
01832 u[3] = u3;
01833 negate = -cdy;
01834 Two_Product(adxtail, negate, ti1, ti0);
01835 negate = -cdytail;
01836 Two_Product(adx, negate, tj1, tj0);
01837 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
01838 v[3] = v3;
01839 catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
01840
01841 Two_Product(cdxtail, adytail, ti1, ti0);
01842 Two_Product(adxtail, cdytail, tj1, tj0);
01843 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
01844 catt[3] = catt3;
01845 cattlen = 4;
01846 } else {
01847 cat[0] = 0.0;
01848 catlen = 1;
01849 catt[0] = 0.0;
01850 cattlen = 1;
01851 }
01852
01853 if (bdxtail != 0.0) {
01854 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
01855 bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
01856 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
01857 temp32a);
01858 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01859 temp32alen, temp32a, temp48);
01860 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01861 temp48, finother);
01862 finswap = finnow; finnow = finother; finother = finswap;
01863 if (cdytail != 0.0) {
01864 temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
01865 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
01866 temp16a);
01867 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
01868 temp16a, finother);
01869 finswap = finnow; finnow = finother; finother = finswap;
01870 }
01871 if (adytail != 0.0) {
01872 temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
01873 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
01874 temp16a);
01875 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
01876 temp16a, finother);
01877 finswap = finnow; finnow = finother; finother = finswap;
01878 }
01879
01880 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
01881 temp32a);
01882 bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
01883 temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
01884 temp16a);
01885 temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
01886 temp16b);
01887 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01888 temp16blen, temp16b, temp32b);
01889 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
01890 temp32blen, temp32b, temp64);
01891 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
01892 temp64, finother);
01893 finswap = finnow; finnow = finother; finother = finswap;
01894 }
01895 if (bdytail != 0.0) {
01896 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
01897 bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
01898 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
01899 temp32a);
01900 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01901 temp32alen, temp32a, temp48);
01902 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01903 temp48, finother);
01904 finswap = finnow; finnow = finother; finother = finswap;
01905
01906
01907 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
01908 temp32a);
01909 bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
01910 temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
01911 temp16a);
01912 temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
01913 temp16b);
01914 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01915 temp16blen, temp16b, temp32b);
01916 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
01917 temp32blen, temp32b, temp64);
01918 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
01919 temp64, finother);
01920 finswap = finnow; finnow = finother; finother = finswap;
01921 }
01922 }
01923 if ((cdxtail != 0.0) || (cdytail != 0.0)) {
01924 if ((adxtail != 0.0) || (adytail != 0.0)
01925 || (bdxtail != 0.0) || (bdytail != 0.0)) {
01926 Two_Product(adxtail, bdy, ti1, ti0);
01927 Two_Product(adx, bdytail, tj1, tj0);
01928 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
01929 u[3] = u3;
01930 negate = -ady;
01931 Two_Product(bdxtail, negate, ti1, ti0);
01932 negate = -adytail;
01933 Two_Product(bdx, negate, tj1, tj0);
01934 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
01935 v[3] = v3;
01936 abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
01937
01938 Two_Product(adxtail, bdytail, ti1, ti0);
01939 Two_Product(bdxtail, adytail, tj1, tj0);
01940 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
01941 abtt[3] = abtt3;
01942 abttlen = 4;
01943 } else {
01944 abt[0] = 0.0;
01945 abtlen = 1;
01946 abtt[0] = 0.0;
01947 abttlen = 1;
01948 }
01949
01950 if (cdxtail != 0.0) {
01951 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
01952 cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
01953 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
01954 temp32a);
01955 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01956 temp32alen, temp32a, temp48);
01957 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
01958 temp48, finother);
01959 finswap = finnow; finnow = finother; finother = finswap;
01960 if (adytail != 0.0) {
01961 temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
01962 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
01963 temp16a);
01964 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
01965 temp16a, finother);
01966 finswap = finnow; finnow = finother; finother = finswap;
01967 }
01968 if (bdytail != 0.0) {
01969 temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
01970 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
01971 temp16a);
01972 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
01973 temp16a, finother);
01974 finswap = finnow; finnow = finother; finother = finswap;
01975 }
01976
01977 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
01978 temp32a);
01979 cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
01980 temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
01981 temp16a);
01982 temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
01983 temp16b);
01984 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01985 temp16blen, temp16b, temp32b);
01986 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
01987 temp32blen, temp32b, temp64);
01988 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
01989 temp64, finother);
01990 finswap = finnow; finnow = finother; finother = finswap;
01991 }
01992 if (cdytail != 0.0) {
01993 temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
01994 cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
01995 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
01996 temp32a);
01997 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
01998 temp32alen, temp32a, temp48);
01999 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02000 temp48, finother);
02001 finswap = finnow; finnow = finother; finother = finswap;
02002
02003
02004 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
02005 temp32a);
02006 cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
02007 temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
02008 temp16a);
02009 temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
02010 temp16b);
02011 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02012 temp16blen, temp16b, temp32b);
02013 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
02014 temp32blen, temp32b, temp64);
02015 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
02016 temp64, finother);
02017 finswap = finnow; finnow = finother; finother = finswap;
02018 }
02019 }
02020
02021 return finnow[finlength - 1];
02022 }
02023
02024 REAL incircle(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
02025 {
02026 REAL adx, bdx, cdx, ady, bdy, cdy;
02027 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
02028 REAL alift, blift, clift;
02029 REAL det;
02030 REAL permanent, errbound;
02031 REAL inc;
02032
02033 FPU_ROUND_DOUBLE;
02034
02035 adx = pa[0] - pd[0];
02036 bdx = pb[0] - pd[0];
02037 cdx = pc[0] - pd[0];
02038 ady = pa[1] - pd[1];
02039 bdy = pb[1] - pd[1];
02040 cdy = pc[1] - pd[1];
02041
02042 bdxcdy = bdx * cdy;
02043 cdxbdy = cdx * bdy;
02044 alift = adx * adx + ady * ady;
02045
02046 cdxady = cdx * ady;
02047 adxcdy = adx * cdy;
02048 blift = bdx * bdx + bdy * bdy;
02049
02050 adxbdy = adx * bdy;
02051 bdxady = bdx * ady;
02052 clift = cdx * cdx + cdy * cdy;
02053
02054 det = alift * (bdxcdy - cdxbdy)
02055 + blift * (cdxady - adxcdy)
02056 + clift * (adxbdy - bdxady);
02057
02058 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
02059 + (Absolute(cdxady) + Absolute(adxcdy)) * blift
02060 + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
02061 errbound = iccerrboundA * permanent;
02062 if ((det > errbound) || (-det > errbound)) {
02063 FPU_RESTORE;
02064 return det;
02065 }
02066
02067 inc = incircleadapt(pa, pb, pc, pd, permanent);
02068 FPU_RESTORE;
02069 return inc;
02070 }
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099 static REAL insphereexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
02100 {
02101 INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1;
02102 INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1;
02103 INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1;
02104 INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1;
02105 REAL axby0, bxcy0, cxdy0, dxey0, exay0;
02106 REAL bxay0, cxby0, dxcy0, exdy0, axey0;
02107 REAL axcy0, bxdy0, cxey0, dxay0, exby0;
02108 REAL cxay0, dxby0, excy0, axdy0, bxey0;
02109 REAL ab[4], bc[4], cd[4], de[4], ea[4];
02110 REAL ac[4], bd[4], ce[4], da[4], eb[4];
02111 REAL temp8a[8], temp8b[8], temp16[16];
02112 int temp8alen, temp8blen, temp16len;
02113 REAL abc[24], bcd[24], cde[24], dea[24], eab[24];
02114 REAL abd[24], bce[24], cda[24], deb[24], eac[24];
02115 int abclen, bcdlen, cdelen, dealen, eablen;
02116 int abdlen, bcelen, cdalen, deblen, eaclen;
02117 REAL temp48a[48], temp48b[48];
02118 int temp48alen, temp48blen;
02119 REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96];
02120 int abcdlen, bcdelen, cdealen, deablen, eabclen;
02121 REAL temp192[192];
02122 REAL det384x[384], det384y[384], det384z[384];
02123 int xlen, ylen, zlen;
02124 REAL detxy[768];
02125 int xylen;
02126 REAL adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152];
02127 int alen, blen, clen, dlen, elen;
02128 REAL abdet[2304], cddet[2304], cdedet[3456];
02129 int ablen, cdlen;
02130 REAL deter[5760];
02131 int deterlen;
02132 int i;
02133
02134 INEXACT REAL bvirt;
02135 REAL avirt, bround, around;
02136 INEXACT REAL c;
02137 INEXACT REAL abig;
02138 REAL ahi, alo, bhi, blo;
02139 REAL err1, err2, err3;
02140 INEXACT REAL _i, _j;
02141 REAL _0;
02142
02143 Two_Product(pa[0], pb[1], axby1, axby0);
02144 Two_Product(pb[0], pa[1], bxay1, bxay0);
02145 Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
02146
02147 Two_Product(pb[0], pc[1], bxcy1, bxcy0);
02148 Two_Product(pc[0], pb[1], cxby1, cxby0);
02149 Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
02150
02151 Two_Product(pc[0], pd[1], cxdy1, cxdy0);
02152 Two_Product(pd[0], pc[1], dxcy1, dxcy0);
02153 Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
02154
02155 Two_Product(pd[0], pe[1], dxey1, dxey0);
02156 Two_Product(pe[0], pd[1], exdy1, exdy0);
02157 Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]);
02158
02159 Two_Product(pe[0], pa[1], exay1, exay0);
02160 Two_Product(pa[0], pe[1], axey1, axey0);
02161 Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]);
02162
02163 Two_Product(pa[0], pc[1], axcy1, axcy0);
02164 Two_Product(pc[0], pa[1], cxay1, cxay0);
02165 Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
02166
02167 Two_Product(pb[0], pd[1], bxdy1, bxdy0);
02168 Two_Product(pd[0], pb[1], dxby1, dxby0);
02169 Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
02170
02171 Two_Product(pc[0], pe[1], cxey1, cxey0);
02172 Two_Product(pe[0], pc[1], excy1, excy0);
02173 Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]);
02174
02175 Two_Product(pd[0], pa[1], dxay1, dxay0);
02176 Two_Product(pa[0], pd[1], axdy1, axdy0);
02177 Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
02178
02179 Two_Product(pe[0], pb[1], exby1, exby0);
02180 Two_Product(pb[0], pe[1], bxey1, bxey0);
02181 Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]);
02182
02183 temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a);
02184 temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b);
02185 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02186 temp16);
02187 temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a);
02188 abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02189 abc);
02190
02191 temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a);
02192 temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b);
02193 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02194 temp16);
02195 temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a);
02196 bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02197 bcd);
02198
02199 temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a);
02200 temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b);
02201 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02202 temp16);
02203 temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a);
02204 cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02205 cde);
02206
02207 temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a);
02208 temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b);
02209 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02210 temp16);
02211 temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a);
02212 dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02213 dea);
02214
02215 temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a);
02216 temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b);
02217 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02218 temp16);
02219 temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a);
02220 eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02221 eab);
02222
02223 temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a);
02224 temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b);
02225 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02226 temp16);
02227 temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a);
02228 abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02229 abd);
02230
02231 temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a);
02232 temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b);
02233 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02234 temp16);
02235 temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a);
02236 bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02237 bce);
02238
02239 temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a);
02240 temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b);
02241 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02242 temp16);
02243 temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a);
02244 cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02245 cda);
02246
02247 temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a);
02248 temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b);
02249 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02250 temp16);
02251 temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a);
02252 deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02253 deb);
02254
02255 temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a);
02256 temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b);
02257 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
02258 temp16);
02259 temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a);
02260 eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
02261 eac);
02262
02263 temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a);
02264 temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b);
02265 for (i = 0; i < temp48blen; i++) {
02266 temp48b[i] = -temp48b[i];
02267 }
02268 bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
02269 temp48blen, temp48b, bcde);
02270 xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192);
02271 xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x);
02272 ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192);
02273 ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y);
02274 zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192);
02275 zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z);
02276 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
02277 alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet);
02278
02279 temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a);
02280 temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b);
02281 for (i = 0; i < temp48blen; i++) {
02282 temp48b[i] = -temp48b[i];
02283 }
02284 cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
02285 temp48blen, temp48b, cdea);
02286 xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192);
02287 xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x);
02288 ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192);
02289 ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y);
02290 zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192);
02291 zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z);
02292 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
02293 blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet);
02294
02295 temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a);
02296 temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b);
02297 for (i = 0; i < temp48blen; i++) {
02298 temp48b[i] = -temp48b[i];
02299 }
02300 deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
02301 temp48blen, temp48b, deab);
02302 xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192);
02303 xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x);
02304 ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192);
02305 ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y);
02306 zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192);
02307 zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z);
02308 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
02309 clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet);
02310
02311 temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a);
02312 temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b);
02313 for (i = 0; i < temp48blen; i++) {
02314 temp48b[i] = -temp48b[i];
02315 }
02316 eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
02317 temp48blen, temp48b, eabc);
02318 xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192);
02319 xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x);
02320 ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192);
02321 ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y);
02322 zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192);
02323 zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z);
02324 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
02325 dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet);
02326
02327 temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a);
02328 temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b);
02329 for (i = 0; i < temp48blen; i++) {
02330 temp48b[i] = -temp48b[i];
02331 }
02332 abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
02333 temp48blen, temp48b, abcd);
02334 xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192);
02335 xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x);
02336 ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192);
02337 ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y);
02338 zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192);
02339 zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z);
02340 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
02341 elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet);
02342
02343 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
02344 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
02345 cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet);
02346 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter);
02347
02348 return deter[deterlen - 1];
02349 }
02350
02351 static REAL insphereadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe,
02352 REAL permanent)
02353 {
02354 INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
02355 REAL det, errbound;
02356
02357 INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1;
02358 INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1;
02359 INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1;
02360 REAL aexbey0, bexaey0, bexcey0, cexbey0;
02361 REAL cexdey0, dexcey0, dexaey0, aexdey0;
02362 REAL aexcey0, cexaey0, bexdey0, dexbey0;
02363 REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
02364 INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3;
02365 REAL abeps, bceps, cdeps, daeps, aceps, bdeps;
02366 REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48];
02367 int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len;
02368 REAL xdet[96], ydet[96], zdet[96], xydet[192];
02369 int xlen, ylen, zlen, xylen;
02370 REAL adet[288], bdet[288], cdet[288], ddet[288];
02371 int alen, blen, clen, dlen;
02372 REAL abdet[576], cddet[576];
02373 int ablen, cdlen;
02374 REAL fin1[1152];
02375 int finlength;
02376
02377 REAL aextail, bextail, cextail, dextail;
02378 REAL aeytail, beytail, ceytail, deytail;
02379 REAL aeztail, beztail, ceztail, deztail;
02380
02381 INEXACT REAL bvirt;
02382 REAL avirt, bround, around;
02383 INEXACT REAL c;
02384 INEXACT REAL abig;
02385 REAL ahi, alo, bhi, blo;
02386 REAL err1, err2, err3;
02387 INEXACT REAL _i, _j;
02388 REAL _0;
02389
02390 aex = (REAL) (pa[0] - pe[0]);
02391 bex = (REAL) (pb[0] - pe[0]);
02392 cex = (REAL) (pc[0] - pe[0]);
02393 dex = (REAL) (pd[0] - pe[0]);
02394 aey = (REAL) (pa[1] - pe[1]);
02395 bey = (REAL) (pb[1] - pe[1]);
02396 cey = (REAL) (pc[1] - pe[1]);
02397 dey = (REAL) (pd[1] - pe[1]);
02398 aez = (REAL) (pa[2] - pe[2]);
02399 bez = (REAL) (pb[2] - pe[2]);
02400 cez = (REAL) (pc[2] - pe[2]);
02401 dez = (REAL) (pd[2] - pe[2]);
02402
02403 Two_Product(aex, bey, aexbey1, aexbey0);
02404 Two_Product(bex, aey, bexaey1, bexaey0);
02405 Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]);
02406 ab[3] = ab3;
02407
02408 Two_Product(bex, cey, bexcey1, bexcey0);
02409 Two_Product(cex, bey, cexbey1, cexbey0);
02410 Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]);
02411 bc[3] = bc3;
02412
02413 Two_Product(cex, dey, cexdey1, cexdey0);
02414 Two_Product(dex, cey, dexcey1, dexcey0);
02415 Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]);
02416 cd[3] = cd3;
02417
02418 Two_Product(dex, aey, dexaey1, dexaey0);
02419 Two_Product(aex, dey, aexdey1, aexdey0);
02420 Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]);
02421 da[3] = da3;
02422
02423 Two_Product(aex, cey, aexcey1, aexcey0);
02424 Two_Product(cex, aey, cexaey1, cexaey0);
02425 Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]);
02426 ac[3] = ac3;
02427
02428 Two_Product(bex, dey, bexdey1, bexdey0);
02429 Two_Product(dex, bey, dexbey1, dexbey0);
02430 Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]);
02431 bd[3] = bd3;
02432
02433 temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a);
02434 temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b);
02435 temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c);
02436 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
02437 temp8blen, temp8b, temp16);
02438 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
02439 temp16len, temp16, temp24);
02440 temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48);
02441 xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet);
02442 temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48);
02443 ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet);
02444 temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48);
02445 zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet);
02446 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
02447 alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet);
02448
02449 temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a);
02450 temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b);
02451 temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c);
02452 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
02453 temp8blen, temp8b, temp16);
02454 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
02455 temp16len, temp16, temp24);
02456 temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48);
02457 xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet);
02458 temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48);
02459 ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet);
02460 temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48);
02461 zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet);
02462 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
02463 blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet);
02464
02465 temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a);
02466 temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b);
02467 temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c);
02468 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
02469 temp8blen, temp8b, temp16);
02470 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
02471 temp16len, temp16, temp24);
02472 temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48);
02473 xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet);
02474 temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48);
02475 ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet);
02476 temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48);
02477 zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet);
02478 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
02479 clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet);
02480
02481 temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a);
02482 temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b);
02483 temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c);
02484 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
02485 temp8blen, temp8b, temp16);
02486 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
02487 temp16len, temp16, temp24);
02488 temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48);
02489 xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet);
02490 temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48);
02491 ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet);
02492 temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48);
02493 zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet);
02494 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
02495 dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet);
02496
02497 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
02498 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
02499 finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1);
02500
02501 det = estimate(finlength, fin1);
02502 errbound = isperrboundB * permanent;
02503 if ((det >= errbound) || (-det >= errbound)) {
02504 return det;
02505 }
02506
02507 Two_Diff_Tail(pa[0], pe[0], aex, aextail);
02508 Two_Diff_Tail(pa[1], pe[1], aey, aeytail);
02509 Two_Diff_Tail(pa[2], pe[2], aez, aeztail);
02510 Two_Diff_Tail(pb[0], pe[0], bex, bextail);
02511 Two_Diff_Tail(pb[1], pe[1], bey, beytail);
02512 Two_Diff_Tail(pb[2], pe[2], bez, beztail);
02513 Two_Diff_Tail(pc[0], pe[0], cex, cextail);
02514 Two_Diff_Tail(pc[1], pe[1], cey, ceytail);
02515 Two_Diff_Tail(pc[2], pe[2], cez, ceztail);
02516 Two_Diff_Tail(pd[0], pe[0], dex, dextail);
02517 Two_Diff_Tail(pd[1], pe[1], dey, deytail);
02518 Two_Diff_Tail(pd[2], pe[2], dez, deztail);
02519 if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0)
02520 && (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0)
02521 && (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0)
02522 && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) {
02523 return det;
02524 }
02525
02526 errbound = isperrboundC * permanent + resulterrbound * Absolute(det);
02527 abeps = (aex * beytail + bey * aextail)
02528 - (aey * bextail + bex * aeytail);
02529 bceps = (bex * ceytail + cey * bextail)
02530 - (bey * cextail + cex * beytail);
02531 cdeps = (cex * deytail + dey * cextail)
02532 - (cey * dextail + dex * ceytail);
02533 daeps = (dex * aeytail + aey * dextail)
02534 - (dey * aextail + aex * deytail);
02535 aceps = (aex * ceytail + cey * aextail)
02536 - (aey * cextail + cex * aeytail);
02537 bdeps = (bex * deytail + dey * bextail)
02538 - (bey * dextail + dex * beytail);
02539 det += (((bex * bex + bey * bey + bez * bez)
02540 * ((cez * daeps + dez * aceps + aez * cdeps)
02541 + (ceztail * da3 + deztail * ac3 + aeztail * cd3))
02542 + (dex * dex + dey * dey + dez * dez)
02543 * ((aez * bceps - bez * aceps + cez * abeps)
02544 + (aeztail * bc3 - beztail * ac3 + ceztail * ab3)))
02545 - ((aex * aex + aey * aey + aez * aez)
02546 * ((bez * cdeps - cez * bdeps + dez * bceps)
02547 + (beztail * cd3 - ceztail * bd3 + deztail * bc3))
02548 + (cex * cex + cey * cey + cez * cez)
02549 * ((dez * abeps + aez * bdeps + bez * daeps)
02550 + (deztail * ab3 + aeztail * bd3 + beztail * da3))))
02551 + 2.0 * (((bex * bextail + bey * beytail + bez * beztail)
02552 * (cez * da3 + dez * ac3 + aez * cd3)
02553 + (dex * dextail + dey * deytail + dez * deztail)
02554 * (aez * bc3 - bez * ac3 + cez * ab3))
02555 - ((aex * aextail + aey * aeytail + aez * aeztail)
02556 * (bez * cd3 - cez * bd3 + dez * bc3)
02557 + (cex * cextail + cey * ceytail + cez * ceztail)
02558 * (dez * ab3 + aez * bd3 + bez * da3)));
02559 if ((det >= errbound) || (-det >= errbound)) {
02560 return det;
02561 }
02562
02563 return insphereexact(pa, pb, pc, pd, pe);
02564 }
02565
02566 REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
02567 {
02568 REAL aex, bex, cex, dex;
02569 REAL aey, bey, cey, dey;
02570 REAL aez, bez, cez, dez;
02571 REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
02572 REAL aexcey, cexaey, bexdey, dexbey;
02573 REAL alift, blift, clift, dlift;
02574 REAL ab, bc, cd, da, ac, bd;
02575 REAL abc, bcd, cda, dab;
02576 REAL aezplus, bezplus, cezplus, dezplus;
02577 REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus;
02578 REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus;
02579 REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus;
02580 REAL det;
02581 REAL permanent, errbound;
02582 REAL ins;
02583
02584 FPU_ROUND_DOUBLE;
02585
02586 aex = pa[0] - pe[0];
02587 bex = pb[0] - pe[0];
02588 cex = pc[0] - pe[0];
02589 dex = pd[0] - pe[0];
02590 aey = pa[1] - pe[1];
02591 bey = pb[1] - pe[1];
02592 cey = pc[1] - pe[1];
02593 dey = pd[1] - pe[1];
02594 aez = pa[2] - pe[2];
02595 bez = pb[2] - pe[2];
02596 cez = pc[2] - pe[2];
02597 dez = pd[2] - pe[2];
02598
02599 aexbey = aex * bey;
02600 bexaey = bex * aey;
02601 ab = aexbey - bexaey;
02602 bexcey = bex * cey;
02603 cexbey = cex * bey;
02604 bc = bexcey - cexbey;
02605 cexdey = cex * dey;
02606 dexcey = dex * cey;
02607 cd = cexdey - dexcey;
02608 dexaey = dex * aey;
02609 aexdey = aex * dey;
02610 da = dexaey - aexdey;
02611
02612 aexcey = aex * cey;
02613 cexaey = cex * aey;
02614 ac = aexcey - cexaey;
02615 bexdey = bex * dey;
02616 dexbey = dex * bey;
02617 bd = bexdey - dexbey;
02618
02619 abc = aez * bc - bez * ac + cez * ab;
02620 bcd = bez * cd - cez * bd + dez * bc;
02621 cda = cez * da + dez * ac + aez * cd;
02622 dab = dez * ab + aez * bd + bez * da;
02623
02624 alift = aex * aex + aey * aey + aez * aez;
02625 blift = bex * bex + bey * bey + bez * bez;
02626 clift = cex * cex + cey * cey + cez * cez;
02627 dlift = dex * dex + dey * dey + dez * dez;
02628
02629 det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
02630
02631 aezplus = Absolute(aez);
02632 bezplus = Absolute(bez);
02633 cezplus = Absolute(cez);
02634 dezplus = Absolute(dez);
02635 aexbeyplus = Absolute(aexbey);
02636 bexaeyplus = Absolute(bexaey);
02637 bexceyplus = Absolute(bexcey);
02638 cexbeyplus = Absolute(cexbey);
02639 cexdeyplus = Absolute(cexdey);
02640 dexceyplus = Absolute(dexcey);
02641 dexaeyplus = Absolute(dexaey);
02642 aexdeyplus = Absolute(aexdey);
02643 aexceyplus = Absolute(aexcey);
02644 cexaeyplus = Absolute(cexaey);
02645 bexdeyplus = Absolute(bexdey);
02646 dexbeyplus = Absolute(dexbey);
02647 permanent = ((cexdeyplus + dexceyplus) * bezplus
02648 + (dexbeyplus + bexdeyplus) * cezplus
02649 + (bexceyplus + cexbeyplus) * dezplus)
02650 * alift
02651 + ((dexaeyplus + aexdeyplus) * cezplus
02652 + (aexceyplus + cexaeyplus) * dezplus
02653 + (cexdeyplus + dexceyplus) * aezplus)
02654 * blift
02655 + ((aexbeyplus + bexaeyplus) * dezplus
02656 + (bexdeyplus + dexbeyplus) * aezplus
02657 + (dexaeyplus + aexdeyplus) * bezplus)
02658 * clift
02659 + ((bexceyplus + cexbeyplus) * aezplus
02660 + (cexaeyplus + aexceyplus) * bezplus
02661 + (aexbeyplus + bexaeyplus) * cezplus)
02662 * dlift;
02663 errbound = isperrboundA * permanent;
02664 if ((det > errbound) || (-det > errbound)) {
02665 FPU_RESTORE;
02666 return det;
02667 }
02668
02669 ins = insphereadapt(pa, pb, pc, pd, pe, permanent);
02670 FPU_RESTORE;
02671 return ins;
02672 }
02673 }