Bundle latest version of Carve library which shall resolve compilation issues with...
[blender.git] / extern / carve / include / carve / exact.hpp
1 // Begin License:
2 // Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
3 // All rights reserved.
4 //
5 // This file is part of the Carve CSG Library (http://carve-csg.com/)
6 //
7 // This file may be used under the terms of the GNU General Public
8 // License version 2.0 as published by the Free Software Foundation
9 // and appearing in the file LICENSE.GPL2 included in the packaging of
10 // this file.
11 //
12 // This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
13 // INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
14 // A PARTICULAR PURPOSE.
15 // End:
16
17
18 #pragma once
19
20 #include <carve/carve.hpp>
21
22 #include <vector>
23 #include <numeric>
24 #include <algorithm>
25
26
27 namespace carve {
28   namespace exact {
29
30     class exact_t : public std::vector<double> {
31       typedef std::vector<double> super;
32
33     public:
34       exact_t() : super() {
35       }
36
37       exact_t(double v, size_t sz = 1) : super(sz, v) {
38       }
39
40       template<typename iter_t>
41       exact_t(iter_t a, iter_t b) : super(a, b) {
42       }
43
44       exact_t(double a, double b) : super() {
45         reserve(2);
46         push_back(a);
47         push_back(b);
48       }
49
50       exact_t(double a, double b, double c) : super() {
51         reserve(3);
52         push_back(a);
53         push_back(b);
54         push_back(c);
55       }
56
57       exact_t(double a, double b, double c, double d) : super() {
58         reserve(4);
59         push_back(a);
60         push_back(b);
61         push_back(c);
62         push_back(d);
63       }
64
65       exact_t(double a, double b, double c, double d, double e) : super() {
66         reserve(5);
67         push_back(a);
68         push_back(b);
69         push_back(c);
70         push_back(d);
71         push_back(e);
72       }
73
74       exact_t(double a, double b, double c, double d, double e, double f) : super() {
75         reserve(6);
76         push_back(a);
77         push_back(b);
78         push_back(c);
79         push_back(d);
80         push_back(e);
81         push_back(f);
82       }
83
84       exact_t(double a, double b, double c, double d, double e, double f, double g) : super() {
85         reserve(7);
86         push_back(a);
87         push_back(b);
88         push_back(c);
89         push_back(d);
90         push_back(e);
91         push_back(f);
92         push_back(g);
93       }
94
95       exact_t(double a, double b, double c, double d, double e, double f, double g, double h) : super() {
96         reserve(8);
97         push_back(a);
98         push_back(b);
99         push_back(c);
100         push_back(d);
101         push_back(e);
102         push_back(f);
103         push_back(g);
104         push_back(h);
105       }
106
107       void compress();
108
109       exact_t compressed() const {
110         exact_t result(*this);
111         result.compress();
112         return result;
113       }
114
115       operator double() const {
116         return std::accumulate(begin(), end(), 0.0);
117       }
118
119       void removeZeroes() {
120         erase(std::remove(begin(), end(), 0.0), end());
121       }
122     };
123
124     inline std::ostream &operator<<(std::ostream &out, const exact_t &p) {
125       out << '{';
126       out << p[0];
127       for (size_t i = 1; i < p.size(); ++i) out << ';' << p[i];
128       out << '}';
129       return out;
130     }
131
132
133
134     namespace detail {
135       const struct constants_t {
136         double splitter;     /* = 2^ceiling(p / 2) + 1.  Used to split floats in half. */
137         double epsilon;                /* = 2^(-p).  Used to estimate roundoff errors. */
138         /* A set of coefficients used to calculate maximum roundoff errors.          */
139         double resulterrbound;
140         double ccwerrboundA, ccwerrboundB, ccwerrboundC;
141         double o3derrboundA, o3derrboundB, o3derrboundC;
142         double iccerrboundA, iccerrboundB, iccerrboundC;
143         double isperrboundA, isperrboundB, isperrboundC;
144
145         constants_t() {
146           double half;
147           double check, lastcheck;
148           int every_other;
149   
150           every_other = 1;
151           half = 0.5;
152           epsilon = 1.0;
153           splitter = 1.0;
154           check = 1.0;
155           /* Repeatedly divide `epsilon' by two until it is too small to add to    */
156           /*   one without causing roundoff.  (Also check if the sum is equal to   */
157           /*   the previous sum, for machines that round up instead of using exact */
158           /*   rounding.  Not that this library will work on such machines anyway. */
159           do {
160             lastcheck = check;
161             epsilon *= half;
162             if (every_other) {
163               splitter *= 2.0;
164             }
165             every_other = !every_other;
166             check = 1.0 + epsilon;
167           } while ((check != 1.0) && (check != lastcheck));
168           splitter += 1.0;
169   
170           /* Error bounds for orientation and incircle tests. */
171           resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
172           ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
173           ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
174           ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
175           o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
176           o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
177           o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
178           iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
179           iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
180           iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
181           isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
182           isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
183           isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
184         }
185       } constants;
186
187       template<unsigned U, unsigned V>
188       struct op {
189         enum {
190           Vlo = V / 2,
191           Vhi = V - Vlo
192         };
193
194         static inline void add(const double *a, const double *b, double *r) {
195           double t[U + Vlo];
196           op<U, Vlo>::add(a, b, t);
197           for (size_t i = 0; i < Vlo; ++i) r[i] = t[i];
198           op<U, Vhi>::add(t + Vlo, b + Vlo, r + Vlo);
199         }
200
201         static inline void sub(const double *a, const double *b, double *r) {
202           double t[U + Vlo];
203           op<U, Vlo>::sub(a, b, t);
204           for (size_t i = 0; i < Vlo; ++i) r[i] = t[i];
205           op<U, Vhi>::sub(t + Vlo, b + Vlo, r + Vlo);
206         }
207       };
208
209       template<unsigned U>
210       struct op<U, 1> {
211         enum {
212           Ulo = U / 2,
213           Uhi = U - Ulo
214         };
215         static void add(const double *a, const double *b, double *r) {
216           double t[Ulo + 1];
217           op<Ulo, 1>::add(a, b, t);
218           for (size_t i = 0; i < Ulo; ++i) r[i] = t[i];
219           op<Uhi, 1>::add(a + Ulo, t + Ulo, r + Ulo);
220         }
221
222         static void sub(const double *a, const double *b, double *r) {
223           double t[Ulo + 1];
224           op<Ulo, 1>::sub(a, b, t);
225           for (size_t i = 0; i < Ulo; ++i) r[i] = t[i];
226           op<Uhi, 1>::add(a + Ulo, t + Ulo, r + Ulo);
227         }
228       };
229
230       template<>
231       struct op<1, 1> {
232         static void add_fast(const double *a, const double *b, double *r) {
233           assert(fabs(a[0]) >= fabs(b[0]));
234           volatile double sum = a[0] + b[0];
235           volatile double bvirt = sum - a[0];
236           r[0] = b[0] - bvirt;
237           r[1] = sum;
238         }
239
240         static void sub_fast(const double *a, const double *b, double *r) {
241           assert(fabs(a[0]) >= fabs(b[0]));
242           volatile double diff = a[0] - b[0];
243           volatile double bvirt = a[0] - diff;
244           r[0] = bvirt - b[0];
245           r[1] = diff;
246         }
247
248         static void add(const double *a, const double *b, double *r) {
249           volatile double sum = a[0] + b[0];
250           volatile double bvirt = sum - a[0];
251           double avirt = sum - bvirt;
252           double bround = b[0] - bvirt;
253           double around = a[0] - avirt;
254           r[0] = around + bround;
255           r[1] = sum;
256         }
257
258         static void sub(const double *a, const double *b, double *r) {
259           volatile double diff = a[0] - b[0];
260           volatile double bvirt = a[0] - diff;
261           double avirt = diff + bvirt;
262           double bround = bvirt - b[0];
263           double around = a[0] - avirt;
264           r[0] = around + bround;
265           r[1] = diff;
266         }
267       };
268
269
270       template<unsigned U, unsigned V>
271       static exact_t add(const double *a, const double *b) {
272         exact_t result;
273         result.resize(U + V);
274         op<U,V>::add(a, b, &result[0]);
275         return result;
276       }
277
278
279       template<unsigned U, unsigned V>
280       static exact_t sub(const double *a, const double *b) {
281         exact_t result;
282         result.resize(U + V);
283         op<U,V>::sub(a, b, &result[0]);
284         return result;
285       }
286
287
288       template<unsigned U, unsigned V>
289       static exact_t add(const exact_t &a, const exact_t &b) {
290         assert(a.size() == U);
291         assert(b.size() == V);
292         exact_t result;
293         result.resize(U + V);
294         std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
295         op<U,V>::add(&a[0], &b[0], &result[0]);
296         return result;
297       }
298
299
300       template<unsigned U, unsigned V>
301       static exact_t add(const exact_t &a, const double *b) {
302         assert(a.size() == U);
303         exact_t result;
304         result.resize(U + V);
305         std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
306         op<U,V>::add(&a[0], b, &result[0]);
307         return result;
308       }
309
310
311       template<unsigned U, unsigned V>
312       static exact_t sub(const exact_t &a, const exact_t &b) {
313         assert(a.size() == U);
314         assert(b.size() == V);
315         exact_t result;
316         result.resize(U + V);
317         std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
318         op<U,V>::sub(&a[0], &b[0], &result[0]);
319         return result;
320       }
321
322
323       template<unsigned U, unsigned V>
324       static exact_t sub(const exact_t &a, const double *b) {
325         assert(a.size() == U);
326         exact_t result;
327         result.resize(U + V);
328         std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
329         op<U,V>::sub(&a[0], &b[0], &result[0]);
330         return result;
331       }
332
333
334       static inline void split(const double a, double *r) {
335         volatile double c = constants.splitter * a;
336         volatile double abig = c - a;
337         r[1] = c - abig;
338         r[0] = a - r[1];
339       }
340
341       static inline void prod_1_1(const double *a, const double *b, double *r) {
342         r[1] = a[0] * b[0];
343         double a_sp[2]; split(a[0], a_sp);
344         double b_sp[2]; split(b[0], b_sp);
345         double err1 = r[1] - a_sp[1] * b_sp[1];
346         double err2 = err1 - a_sp[0] * b_sp[1];
347         double err3 = err2 - a_sp[1] * b_sp[0];
348         r[0] = a_sp[0] * b_sp[0] - err3;
349       }
350
351       static inline void prod_1_1s(const double *a, const double *b, const double *b_sp, double *r) {
352         r[1] = a[0] * b[0];
353         double a_sp[2]; split(a[0], a_sp);
354         double err1 = r[1] - a_sp[1] * b_sp[1];
355         double err2 = err1 - a_sp[0] * b_sp[1];
356         double err3 = err2 - a_sp[1] * b_sp[0];
357         r[0] = a_sp[0] * b_sp[0] - err3;
358       }
359
360       static inline void prod_1s_1s(const double *a, const double *a_sp, const double *b, const double *b_sp, double *r) {
361         r[1] = a[0] * b[0];
362         double err1 = r[1] - a_sp[1] * b_sp[1];
363         double err2 = err1 - a_sp[0] * b_sp[1];
364         double err3 = err2 - a_sp[1] * b_sp[0];
365         r[0] = a_sp[0] * b_sp[0] - err3;
366       }
367
368       static inline void prod_2_1(const double *a, const double *b, double *r) {
369         double b_sp[2]; split(b[0], b_sp);
370         double t1[2]; prod_1_1s(a+0, b, b_sp, t1);
371         r[0] = t1[0];
372         double t2[2]; prod_1_1s(a+1, b, b_sp, t2);
373         double t3[2]; op<1,1>::add(t1+1, t2, t3);
374         r[1] = t3[0];
375         double t4[2]; op<1,1>::add_fast(t2+1, t3+1, r + 2);
376       }
377
378       static inline void prod_1_2(const double *a, const double *b, double *r) {
379         prod_2_1(b, a, r);
380       }
381
382       static inline void prod_4_1(const double *a, const double *b, double *r) {
383         double b_sp[2]; split(b[0], b_sp);
384         double t1[2]; prod_1_1s(a+0, b, b_sp, t1);
385         r[0] = t1[0];
386         double t2[2]; prod_1_1s(a+1, b, b_sp, t2);
387         double t3[2]; op<1,1>::add(t1+1, t2, t3);
388         r[1] = t3[0];
389         double t4[2]; op<1,1>::add_fast(t2+1, t3+1, t4);
390         r[2] = t4[0];
391         double t5[2]; prod_1_1s(a+2, b, b_sp, t5);
392         double t6[2]; op<1,1>::add(t4+1, t5, t6);
393         r[3] = t6[0];
394         double t7[2]; op<1,1>::add_fast(t5+1, t6+1, t7);
395         r[4] = t7[0];
396         double t8[2]; prod_1_1s(a+3, b, b_sp, t8);
397         double t9[2]; op<1,1>::add(t7+1, t8, t9);
398         r[5] = t9[0];
399         op<1,1>::add_fast(t8+1, t9+1, r + 6);
400       }
401
402       static inline void prod_1_4(const double *a, const double *b, double *r) {
403         prod_4_1(b, a, r);
404       }
405
406       static inline void prod_2_2(const double *a, const double *b, double *r) {
407         double a1_sp[2]; split(a[1], a1_sp);
408         double a0_sp[2]; split(a[0], a0_sp);
409         double b1_sp[2]; split(b[1], b1_sp);
410         double b0_sp[2]; split(b[0], b0_sp);
411
412         double t1[2]; prod_1s_1s(a+0, a0_sp, b+0, b0_sp, t1);
413         r[0] = t1[0];
414         double t2[2]; prod_1s_1s(a+1, a1_sp, b+0, b0_sp, t2);
415
416         double t3[2]; op<1,1>::add(t1+1, t2, t3);
417         double t4[2]; op<1,1>::add_fast(t2+1, t3+1, t4);
418
419         double t5[2]; prod_1s_1s(a+0, a0_sp, b+1, b1_sp, t5);
420
421         double t6[2]; op<1,1>::add(t3, t5, t6);
422         r[1] = t6[0];
423         double t7[2]; op<1,1>::add(t4, t6+1, t7);
424         double t8[2]; op<1,1>::add(t4+1, t7+1, t8);
425
426         double t9[2]; prod_1s_1s(a+1, a1_sp, b+1, b1_sp, t9);
427
428         double t10[2]; op<1,1>::add(t5+1, t9, t10);
429         double t11[2]; op<1,1>::add(t7, t10, t11);
430         r[2] = t11[0];
431         double t12[2]; op<1,1>::add(t8, t11+1, t12);
432         double t13[2]; op<1,1>::add(t8+1, t12+1, t13);
433         double t14[2]; op<1,1>::add(t9+1, t10+1, t14);
434         double t15[2]; op<1,1>::add(t12, t14, t15);
435         r[3] = t15[0];
436         double t16[2]; op<1,1>::add(t13, t15+1, t16);
437         double t17[2]; op<1,1>::add(t13+1, t16+1, t17);
438         double t18[2]; op<1,1>::add(t16, t14+1, t18);
439         r[4] = t18[0];
440         double t19[2]; op<1,1>::add(t17, t18+1, t19);
441         r[5] = t19[0];
442         double t20[2]; op<1,1>::add(t17+1, t19+1, t20);
443         r[6] = t20[0];
444         r[7] = t20[1];
445       }
446
447
448
449       static inline void square(const double a, double *r) {
450         r[1] = a * a;
451         double a_sp[2]; split(a, a_sp);
452         double err1 = r[1] - (a_sp[1] * a_sp[1]);
453         double err3 = err1 - ((a_sp[1] + a_sp[1]) * a_sp[0]);
454         r[0] = a_sp[0] * a_sp[0] - err3;
455       }
456
457       static inline void square_2(const double *a, double *r) {
458         double t1[2]; square(a[0], t1);
459         r[0] = t1[0];
460         double t2 = a[0] + a[0];
461         double t3[2]; prod_1_1(a+1, &t2, t3);
462         double t4[3]; op<2,1>::add(t3, t1 + 1, t4);
463         r[1] = t4[0];
464         double t5[2]; square(a[1], t5);
465         double t6[4]; op<2,2>::add(t5, t4 + 1, r + 2);
466       }
467     }
468
469
470
471     void exact_t::compress() {
472       double sum[2];
473
474       int j = size() - 1;
475       double Q = (*this)[j];
476       for (int i = (int)size()-2; i >= 0; --i) {
477         detail::op<1,1>::add_fast(&Q, &(*this)[i], sum);
478         if (sum[0] != 0) {
479           (*this)[j--] = sum[1];
480           Q = sum[0];
481         } else {
482           Q = sum[1];
483         }
484       }
485       int j2 = 0;
486       for (int i = j + 1; i < (int)size(); ++i) {
487         detail::op<1,1>::add_fast(&(*this)[i], &Q, sum);
488         if (sum[0] != 0) {
489           (*this)[j2++] = sum[0];
490         }
491         Q = sum[1];
492       }
493       (*this)[j2++] = Q;
494
495       erase(begin() + j2, end());
496     }
497
498     template<typename iter_t>
499     void negate(iter_t begin, iter_t end) {
500       while (begin != end) { *begin = -*begin; ++begin; }
501     }
502
503     void negate(exact_t &e) {
504       negate(&e[0], &e[e.size()]);
505     }
506
507     template<typename iter_t>
508     void scale_zeroelim(iter_t ebegin,
509                         iter_t eend,
510                         double b,
511                         exact_t &h) {
512       double Q;
513
514       h.clear();
515       double b_sp[2]; detail::split(b, b_sp);
516
517       double prod[2], sum[2];
518
519       detail::prod_1_1s((double *)ebegin++, &b, b_sp, prod);
520       Q = prod[1];
521       if (prod[0] != 0.0) {
522         h.push_back(prod[0]);
523       }
524       while (ebegin != eend) {
525         double enow = *ebegin++;
526         detail::prod_1_1s(&enow, &b, b_sp, prod);
527         detail::op<1,1>::add(&Q, prod, sum);
528         if (sum[0] != 0) {
529           h.push_back(sum[0]);
530         }
531         detail::op<1,1>::add_fast(prod+1, sum+1, sum);
532         Q = sum[1];
533         if (sum[0] != 0) {
534           h.push_back(sum[0]);
535         }
536       }
537       if ((Q != 0.0) || (h.size() == 0)) {
538         h.push_back(Q);
539       }
540     }
541
542     void scale_zeroelim(const exact_t &e,
543                         double b,
544                         exact_t &h) {
545       scale_zeroelim(&e[0], &e[e.size()], b, h);
546     }
547
548     template<typename iter_t>
549     void sum_zeroelim(iter_t ebegin,
550                       iter_t eend,
551                       iter_t fbegin,
552                       iter_t fend,
553                       exact_t &h) {
554       double Q;
555       double enow, fnow;
556
557       double sum[2];
558
559       enow = *ebegin;
560       fnow = *fbegin;
561
562       h.clear();
563
564       if ((fnow > enow) == (fnow > -enow)) {
565         Q = enow;
566         enow = *++ebegin;
567       } else {
568         Q = fnow;
569         fnow = *++fbegin;
570       }
571
572       if (ebegin != eend && fbegin != fend) {
573         if ((fnow > enow) == (fnow > -enow)) {
574           detail::op<1,1>::add_fast(&enow, &Q, sum);
575           enow = *++ebegin;
576         } else {
577           detail::op<1,1>::add_fast(&fnow, &Q, sum);
578           fnow = *++fbegin;
579         }
580         Q = sum[1];
581         if (sum[0] != 0.0) {
582           h.push_back(sum[0]);
583         }
584         while (ebegin != eend && fbegin != fend) {
585           if ((fnow > enow) == (fnow > -enow)) {
586             detail::op<1,1>::add(&Q, &enow, sum);
587             enow = *++ebegin;
588           } else {
589             detail::op<1,1>::add(&Q, &fnow, sum);
590             fnow = *++fbegin;
591           }
592           Q = sum[1];
593           if (sum[0] != 0.0) {
594             h.push_back(sum[0]);
595           }
596         }
597       }
598
599       while (ebegin != eend) {
600         detail::op<1,1>::add(&Q, &enow, sum);
601         enow = *++ebegin;
602         Q = sum[1];
603         if (sum[0] != 0.0) {
604           h.push_back(sum[0]);
605         }
606       }
607       while (fbegin != fend) {
608         detail::op<1,1>::add(&Q, &fnow, sum);
609         fnow = *++fbegin;
610         Q = sum[1];
611         if (sum[0] != 0.0) {
612           h.push_back(sum[0]);
613         }
614       }
615
616       if (Q != 0.0 || !h.size()) {
617         h.push_back(Q);
618       }
619     }
620
621     void sum_zeroelim(const exact_t &e,
622                       const exact_t &f,
623                       exact_t &h) {
624       sum_zeroelim(&e[0], &e[e.size()], &f[0], &f[f.size()], h);
625     }
626
627     void sum_zeroelim(const double *ebegin,
628                       const double *eend,
629                       const exact_t &f,
630                       exact_t &h) {
631       sum_zeroelim(ebegin, eend, &f[0], &f[f.size()], h);
632     }
633
634     void sum_zeroelim(const exact_t &e,
635                       const double *fbegin,
636                       const double *fend,
637                       exact_t &h) {
638       sum_zeroelim(&e[0], &e[e.size()], fbegin, fend, h);
639     }
640
641
642     exact_t operator+(const exact_t &a, const exact_t &b) {
643       exact_t r;
644       sum_zeroelim(a, b, r);
645       return r;
646     }
647
648
649
650     void diffprod(const double a, const double b, const double c, const double d, double *r) {
651       // return ab - cd;
652       double ab[2], cd[2];
653       detail::prod_1_1(&a, &b, ab);
654       detail::prod_1_1(&c, &d, cd);
655       detail::op<2,2>::sub(ab, cd, r);
656     }
657
658     double orient3dexact(const double *pa,
659                          const double *pb,
660                          const double *pc,
661                          const double *pd) {
662       using namespace detail;
663
664       double ab[4]; diffprod(pa[0], pb[1], pb[0], pa[1], ab);
665       double bc[4]; diffprod(pb[0], pc[1], pc[0], pb[1], bc);
666       double cd[4]; diffprod(pc[0], pd[1], pd[0], pc[1], cd);
667       double da[4]; diffprod(pd[0], pa[1], pa[0], pd[1], da);
668       double ac[4]; diffprod(pa[0], pc[1], pc[0], pa[1], ac);
669       double bd[4]; diffprod(pb[0], pd[1], pd[0], pb[1], bd);
670
671       exact_t temp;
672       exact_t cda, dab, abc, bcd;
673       exact_t adet, bdet, cdet, ddet, abdet, cddet, det;
674
675       sum_zeroelim(cd, cd + 4, da, da + 4, temp);
676       sum_zeroelim(temp, ac, ac + 4, cda);
677
678       sum_zeroelim(da, da + 4, ab, ab + 4, temp);
679       sum_zeroelim(temp, bd, bd + 4, dab);
680
681       negate(bd, bd + 4);
682       negate(ac, bd + 4);
683
684       sum_zeroelim(ab, ab + 4, bc, bc + 4, temp);
685       sum_zeroelim(temp, ac, ac + 4, abc);
686
687       sum_zeroelim(bc, bc + 4, cd, cd + 4, temp);
688       sum_zeroelim(temp, bd, bd + 4, bcd);
689
690       scale_zeroelim(bcd, +pa[2], adet);
691       scale_zeroelim(cda, -pb[2], bdet);
692       scale_zeroelim(dab, +pc[2], cdet);
693       scale_zeroelim(abc, -pd[2], ddet);
694
695       sum_zeroelim(adet, bdet, abdet);
696       sum_zeroelim(cdet, ddet, cddet);
697
698       sum_zeroelim(abdet, cddet, det);
699   
700       return det[det.size() - 1];
701     }
702
703   }
704 }