Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/extern/carve/include/carve/exact.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'extern/carve/include/carve/exact.hpp')
-rw-r--r--extern/carve/include/carve/exact.hpp704
1 files changed, 0 insertions, 704 deletions
diff --git a/extern/carve/include/carve/exact.hpp b/extern/carve/include/carve/exact.hpp
deleted file mode 100644
index a74fb77d2ec..00000000000
--- a/extern/carve/include/carve/exact.hpp
+++ /dev/null
@@ -1,704 +0,0 @@
-// Begin License:
-// Copyright (C) 2006-2014 Tobias Sargeant (tobias.sargeant@gmail.com).
-// All rights reserved.
-//
-// This file is part of the Carve CSG Library (http://carve-csg.com/)
-//
-// This file may be used under the terms of either the GNU General
-// Public License version 2 or 3 (at your option) as published by the
-// Free Software Foundation and appearing in the files LICENSE.GPL2
-// and LICENSE.GPL3 included in the packaging of this file.
-//
-// This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
-// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
-// A PARTICULAR PURPOSE.
-// End:
-
-
-#pragma once
-
-#include <carve/carve.hpp>
-
-#include <vector>
-#include <numeric>
-#include <algorithm>
-
-
-namespace carve {
- namespace exact {
-
- class exact_t : public std::vector<double> {
- typedef std::vector<double> super;
-
- public:
- exact_t() : super() {
- }
-
- exact_t(double v, size_t sz = 1) : super(sz, v) {
- }
-
- template<typename iter_t>
- exact_t(iter_t a, iter_t b) : super(a, b) {
- }
-
- exact_t(double a, double b) : super() {
- reserve(2);
- push_back(a);
- push_back(b);
- }
-
- exact_t(double a, double b, double c) : super() {
- reserve(3);
- push_back(a);
- push_back(b);
- push_back(c);
- }
-
- exact_t(double a, double b, double c, double d) : super() {
- reserve(4);
- push_back(a);
- push_back(b);
- push_back(c);
- push_back(d);
- }
-
- exact_t(double a, double b, double c, double d, double e) : super() {
- reserve(5);
- push_back(a);
- push_back(b);
- push_back(c);
- push_back(d);
- push_back(e);
- }
-
- exact_t(double a, double b, double c, double d, double e, double f) : super() {
- reserve(6);
- push_back(a);
- push_back(b);
- push_back(c);
- push_back(d);
- push_back(e);
- push_back(f);
- }
-
- exact_t(double a, double b, double c, double d, double e, double f, double g) : super() {
- reserve(7);
- push_back(a);
- push_back(b);
- push_back(c);
- push_back(d);
- push_back(e);
- push_back(f);
- push_back(g);
- }
-
- exact_t(double a, double b, double c, double d, double e, double f, double g, double h) : super() {
- reserve(8);
- push_back(a);
- push_back(b);
- push_back(c);
- push_back(d);
- push_back(e);
- push_back(f);
- push_back(g);
- push_back(h);
- }
-
- void compress();
-
- exact_t compressed() const {
- exact_t result(*this);
- result.compress();
- return result;
- }
-
- operator double() const {
- return std::accumulate(begin(), end(), 0.0);
- }
-
- void removeZeroes() {
- erase(std::remove(begin(), end(), 0.0), end());
- }
- };
-
- inline std::ostream &operator<<(std::ostream &out, const exact_t &p) {
- out << '{';
- out << p[0];
- for (size_t i = 1; i < p.size(); ++i) out << ';' << p[i];
- out << '}';
- return out;
- }
-
-
-
- namespace detail {
- const struct constants_t {
- double splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */
- double epsilon; /* = 2^(-p). Used to estimate roundoff errors. */
- /* A set of coefficients used to calculate maximum roundoff errors. */
- double resulterrbound;
- double ccwerrboundA, ccwerrboundB, ccwerrboundC;
- double o3derrboundA, o3derrboundB, o3derrboundC;
- double iccerrboundA, iccerrboundB, iccerrboundC;
- double isperrboundA, isperrboundB, isperrboundC;
-
- constants_t() {
- double half;
- double check, lastcheck;
- int every_other;
-
- every_other = 1;
- half = 0.5;
- epsilon = 1.0;
- splitter = 1.0;
- check = 1.0;
- /* Repeatedly divide `epsilon' by two until it is too small to add to */
- /* one without causing roundoff. (Also check if the sum is equal to */
- /* the previous sum, for machines that round up instead of using exact */
- /* rounding. Not that this library will work on such machines anyway. */
- do {
- lastcheck = check;
- epsilon *= half;
- if (every_other) {
- splitter *= 2.0;
- }
- every_other = !every_other;
- check = 1.0 + epsilon;
- } while ((check != 1.0) && (check != lastcheck));
- splitter += 1.0;
-
- /* Error bounds for orientation and incircle tests. */
- resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
- ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
- ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
- ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
- o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
- o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
- o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
- iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
- iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
- iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
- isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
- isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
- isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
- }
- } constants;
-
- template<unsigned U, unsigned V>
- struct op {
- enum {
- Vlo = V / 2,
- Vhi = V - Vlo
- };
-
- static inline void add(const double *a, const double *b, double *r) {
- double t[U + Vlo];
- op<U, Vlo>::add(a, b, t);
- for (size_t i = 0; i < Vlo; ++i) r[i] = t[i];
- op<U, Vhi>::add(t + Vlo, b + Vlo, r + Vlo);
- }
-
- static inline void sub(const double *a, const double *b, double *r) {
- double t[U + Vlo];
- op<U, Vlo>::sub(a, b, t);
- for (size_t i = 0; i < Vlo; ++i) r[i] = t[i];
- op<U, Vhi>::sub(t + Vlo, b + Vlo, r + Vlo);
- }
- };
-
- template<unsigned U>
- struct op<U, 1> {
- enum {
- Ulo = U / 2,
- Uhi = U - Ulo
- };
- static void add(const double *a, const double *b, double *r) {
- double t[Ulo + 1];
- op<Ulo, 1>::add(a, b, t);
- for (size_t i = 0; i < Ulo; ++i) r[i] = t[i];
- op<Uhi, 1>::add(a + Ulo, t + Ulo, r + Ulo);
- }
-
- static void sub(const double *a, const double *b, double *r) {
- double t[Ulo + 1];
- op<Ulo, 1>::sub(a, b, t);
- for (size_t i = 0; i < Ulo; ++i) r[i] = t[i];
- op<Uhi, 1>::add(a + Ulo, t + Ulo, r + Ulo);
- }
- };
-
- template<>
- struct op<1, 1> {
- static void add_fast(const double *a, const double *b, double *r) {
- assert(fabs(a[0]) >= fabs(b[0]));
- volatile double sum = a[0] + b[0];
- volatile double bvirt = sum - a[0];
- r[0] = b[0] - bvirt;
- r[1] = sum;
- }
-
- static void sub_fast(const double *a, const double *b, double *r) {
- assert(fabs(a[0]) >= fabs(b[0]));
- volatile double diff = a[0] - b[0];
- volatile double bvirt = a[0] - diff;
- r[0] = bvirt - b[0];
- r[1] = diff;
- }
-
- static void add(const double *a, const double *b, double *r) {
- volatile double sum = a[0] + b[0];
- volatile double bvirt = sum - a[0];
- double avirt = sum - bvirt;
- double bround = b[0] - bvirt;
- double around = a[0] - avirt;
- r[0] = around + bround;
- r[1] = sum;
- }
-
- static void sub(const double *a, const double *b, double *r) {
- volatile double diff = a[0] - b[0];
- volatile double bvirt = a[0] - diff;
- double avirt = diff + bvirt;
- double bround = bvirt - b[0];
- double around = a[0] - avirt;
- r[0] = around + bround;
- r[1] = diff;
- }
- };
-
-
- template<unsigned U, unsigned V>
- static exact_t add(const double *a, const double *b) {
- exact_t result;
- result.resize(U + V);
- op<U,V>::add(a, b, &result[0]);
- return result;
- }
-
-
- template<unsigned U, unsigned V>
- static exact_t sub(const double *a, const double *b) {
- exact_t result;
- result.resize(U + V);
- op<U,V>::sub(a, b, &result[0]);
- return result;
- }
-
-
- template<unsigned U, unsigned V>
- static exact_t add(const exact_t &a, const exact_t &b) {
- assert(a.size() == U);
- assert(b.size() == V);
- exact_t result;
- result.resize(U + V);
- std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
- op<U,V>::add(&a[0], &b[0], &result[0]);
- return result;
- }
-
-
- template<unsigned U, unsigned V>
- static exact_t add(const exact_t &a, const double *b) {
- assert(a.size() == U);
- exact_t result;
- result.resize(U + V);
- std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
- op<U,V>::add(&a[0], b, &result[0]);
- return result;
- }
-
-
- template<unsigned U, unsigned V>
- static exact_t sub(const exact_t &a, const exact_t &b) {
- assert(a.size() == U);
- assert(b.size() == V);
- exact_t result;
- result.resize(U + V);
- std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
- op<U,V>::sub(&a[0], &b[0], &result[0]);
- return result;
- }
-
-
- template<unsigned U, unsigned V>
- static exact_t sub(const exact_t &a, const double *b) {
- assert(a.size() == U);
- exact_t result;
- result.resize(U + V);
- std::fill(result.begin(), result.end(), std::numeric_limits<double>::quiet_NaN());
- op<U,V>::sub(&a[0], &b[0], &result[0]);
- return result;
- }
-
-
- static inline void split(const double a, double *r) {
- volatile double c = constants.splitter * a;
- volatile double abig = c - a;
- r[1] = c - abig;
- r[0] = a - r[1];
- }
-
- static inline void prod_1_1(const double *a, const double *b, double *r) {
- r[1] = a[0] * b[0];
- double a_sp[2]; split(a[0], a_sp);
- double b_sp[2]; split(b[0], b_sp);
- double err1 = r[1] - a_sp[1] * b_sp[1];
- double err2 = err1 - a_sp[0] * b_sp[1];
- double err3 = err2 - a_sp[1] * b_sp[0];
- r[0] = a_sp[0] * b_sp[0] - err3;
- }
-
- static inline void prod_1_1s(const double *a, const double *b, const double *b_sp, double *r) {
- r[1] = a[0] * b[0];
- double a_sp[2]; split(a[0], a_sp);
- double err1 = r[1] - a_sp[1] * b_sp[1];
- double err2 = err1 - a_sp[0] * b_sp[1];
- double err3 = err2 - a_sp[1] * b_sp[0];
- r[0] = a_sp[0] * b_sp[0] - err3;
- }
-
- static inline void prod_1s_1s(const double *a, const double *a_sp, const double *b, const double *b_sp, double *r) {
- r[1] = a[0] * b[0];
- double err1 = r[1] - a_sp[1] * b_sp[1];
- double err2 = err1 - a_sp[0] * b_sp[1];
- double err3 = err2 - a_sp[1] * b_sp[0];
- r[0] = a_sp[0] * b_sp[0] - err3;
- }
-
- static inline void prod_2_1(const double *a, const double *b, double *r) {
- double b_sp[2]; split(b[0], b_sp);
- double t1[2]; prod_1_1s(a+0, b, b_sp, t1);
- r[0] = t1[0];
- double t2[2]; prod_1_1s(a+1, b, b_sp, t2);
- double t3[2]; op<1,1>::add(t1+1, t2, t3);
- r[1] = t3[0];
- double t4[2]; op<1,1>::add_fast(t2+1, t3+1, r + 2);
- }
-
- static inline void prod_1_2(const double *a, const double *b, double *r) {
- prod_2_1(b, a, r);
- }
-
- static inline void prod_4_1(const double *a, const double *b, double *r) {
- double b_sp[2]; split(b[0], b_sp);
- double t1[2]; prod_1_1s(a+0, b, b_sp, t1);
- r[0] = t1[0];
- double t2[2]; prod_1_1s(a+1, b, b_sp, t2);
- double t3[2]; op<1,1>::add(t1+1, t2, t3);
- r[1] = t3[0];
- double t4[2]; op<1,1>::add_fast(t2+1, t3+1, t4);
- r[2] = t4[0];
- double t5[2]; prod_1_1s(a+2, b, b_sp, t5);
- double t6[2]; op<1,1>::add(t4+1, t5, t6);
- r[3] = t6[0];
- double t7[2]; op<1,1>::add_fast(t5+1, t6+1, t7);
- r[4] = t7[0];
- double t8[2]; prod_1_1s(a+3, b, b_sp, t8);
- double t9[2]; op<1,1>::add(t7+1, t8, t9);
- r[5] = t9[0];
- op<1,1>::add_fast(t8+1, t9+1, r + 6);
- }
-
- static inline void prod_1_4(const double *a, const double *b, double *r) {
- prod_4_1(b, a, r);
- }
-
- static inline void prod_2_2(const double *a, const double *b, double *r) {
- double a1_sp[2]; split(a[1], a1_sp);
- double a0_sp[2]; split(a[0], a0_sp);
- double b1_sp[2]; split(b[1], b1_sp);
- double b0_sp[2]; split(b[0], b0_sp);
-
- double t1[2]; prod_1s_1s(a+0, a0_sp, b+0, b0_sp, t1);
- r[0] = t1[0];
- double t2[2]; prod_1s_1s(a+1, a1_sp, b+0, b0_sp, t2);
-
- double t3[2]; op<1,1>::add(t1+1, t2, t3);
- double t4[2]; op<1,1>::add_fast(t2+1, t3+1, t4);
-
- double t5[2]; prod_1s_1s(a+0, a0_sp, b+1, b1_sp, t5);
-
- double t6[2]; op<1,1>::add(t3, t5, t6);
- r[1] = t6[0];
- double t7[2]; op<1,1>::add(t4, t6+1, t7);
- double t8[2]; op<1,1>::add(t4+1, t7+1, t8);
-
- double t9[2]; prod_1s_1s(a+1, a1_sp, b+1, b1_sp, t9);
-
- double t10[2]; op<1,1>::add(t5+1, t9, t10);
- double t11[2]; op<1,1>::add(t7, t10, t11);
- r[2] = t11[0];
- double t12[2]; op<1,1>::add(t8, t11+1, t12);
- double t13[2]; op<1,1>::add(t8+1, t12+1, t13);
- double t14[2]; op<1,1>::add(t9+1, t10+1, t14);
- double t15[2]; op<1,1>::add(t12, t14, t15);
- r[3] = t15[0];
- double t16[2]; op<1,1>::add(t13, t15+1, t16);
- double t17[2]; op<1,1>::add(t13+1, t16+1, t17);
- double t18[2]; op<1,1>::add(t16, t14+1, t18);
- r[4] = t18[0];
- double t19[2]; op<1,1>::add(t17, t18+1, t19);
- r[5] = t19[0];
- double t20[2]; op<1,1>::add(t17+1, t19+1, t20);
- r[6] = t20[0];
- r[7] = t20[1];
- }
-
-
-
- static inline void square(const double a, double *r) {
- r[1] = a * a;
- double a_sp[2]; split(a, a_sp);
- double err1 = r[1] - (a_sp[1] * a_sp[1]);
- double err3 = err1 - ((a_sp[1] + a_sp[1]) * a_sp[0]);
- r[0] = a_sp[0] * a_sp[0] - err3;
- }
-
- static inline void square_2(const double *a, double *r) {
- double t1[2]; square(a[0], t1);
- r[0] = t1[0];
- double t2 = a[0] + a[0];
- double t3[2]; prod_1_1(a+1, &t2, t3);
- double t4[3]; op<2,1>::add(t3, t1 + 1, t4);
- r[1] = t4[0];
- double t5[2]; square(a[1], t5);
- double t6[4]; op<2,2>::add(t5, t4 + 1, r + 2);
- }
- }
-
-
-
- void exact_t::compress() {
- double sum[2];
-
- int j = size() - 1;
- double Q = (*this)[j];
- for (int i = (int)size()-2; i >= 0; --i) {
- detail::op<1,1>::add_fast(&Q, &(*this)[i], sum);
- if (sum[0] != 0) {
- (*this)[j--] = sum[1];
- Q = sum[0];
- } else {
- Q = sum[1];
- }
- }
- int j2 = 0;
- for (int i = j + 1; i < (int)size(); ++i) {
- detail::op<1,1>::add_fast(&(*this)[i], &Q, sum);
- if (sum[0] != 0) {
- (*this)[j2++] = sum[0];
- }
- Q = sum[1];
- }
- (*this)[j2++] = Q;
-
- erase(begin() + j2, end());
- }
-
- template<typename iter_t>
- void negate(iter_t begin, iter_t end) {
- while (begin != end) { *begin = -*begin; ++begin; }
- }
-
- void negate(exact_t &e) {
- negate(&e[0], &e[e.size()]);
- }
-
- template<typename iter_t>
- void scale_zeroelim(iter_t ebegin,
- iter_t eend,
- double b,
- exact_t &h) {
- double Q;
-
- h.clear();
- double b_sp[2]; detail::split(b, b_sp);
-
- double prod[2], sum[2];
-
- detail::prod_1_1s((double *)ebegin++, &b, b_sp, prod);
- Q = prod[1];
- if (prod[0] != 0.0) {
- h.push_back(prod[0]);
- }
- while (ebegin != eend) {
- double enow = *ebegin++;
- detail::prod_1_1s(&enow, &b, b_sp, prod);
- detail::op<1,1>::add(&Q, prod, sum);
- if (sum[0] != 0) {
- h.push_back(sum[0]);
- }
- detail::op<1,1>::add_fast(prod+1, sum+1, sum);
- Q = sum[1];
- if (sum[0] != 0) {
- h.push_back(sum[0]);
- }
- }
- if ((Q != 0.0) || (h.size() == 0)) {
- h.push_back(Q);
- }
- }
-
- void scale_zeroelim(const exact_t &e,
- double b,
- exact_t &h) {
- scale_zeroelim(&e[0], &e[e.size()], b, h);
- }
-
- template<typename iter_t>
- void sum_zeroelim(iter_t ebegin,
- iter_t eend,
- iter_t fbegin,
- iter_t fend,
- exact_t &h) {
- double Q;
- double enow, fnow;
-
- double sum[2];
-
- enow = *ebegin;
- fnow = *fbegin;
-
- h.clear();
-
- if ((fnow > enow) == (fnow > -enow)) {
- Q = enow;
- enow = *++ebegin;
- } else {
- Q = fnow;
- fnow = *++fbegin;
- }
-
- if (ebegin != eend && fbegin != fend) {
- if ((fnow > enow) == (fnow > -enow)) {
- detail::op<1,1>::add_fast(&enow, &Q, sum);
- enow = *++ebegin;
- } else {
- detail::op<1,1>::add_fast(&fnow, &Q, sum);
- fnow = *++fbegin;
- }
- Q = sum[1];
- if (sum[0] != 0.0) {
- h.push_back(sum[0]);
- }
- while (ebegin != eend && fbegin != fend) {
- if ((fnow > enow) == (fnow > -enow)) {
- detail::op<1,1>::add(&Q, &enow, sum);
- enow = *++ebegin;
- } else {
- detail::op<1,1>::add(&Q, &fnow, sum);
- fnow = *++fbegin;
- }
- Q = sum[1];
- if (sum[0] != 0.0) {
- h.push_back(sum[0]);
- }
- }
- }
-
- while (ebegin != eend) {
- detail::op<1,1>::add(&Q, &enow, sum);
- enow = *++ebegin;
- Q = sum[1];
- if (sum[0] != 0.0) {
- h.push_back(sum[0]);
- }
- }
- while (fbegin != fend) {
- detail::op<1,1>::add(&Q, &fnow, sum);
- fnow = *++fbegin;
- Q = sum[1];
- if (sum[0] != 0.0) {
- h.push_back(sum[0]);
- }
- }
-
- if (Q != 0.0 || !h.size()) {
- h.push_back(Q);
- }
- }
-
- void sum_zeroelim(const exact_t &e,
- const exact_t &f,
- exact_t &h) {
- sum_zeroelim(&e[0], &e[e.size()], &f[0], &f[f.size()], h);
- }
-
- void sum_zeroelim(const double *ebegin,
- const double *eend,
- const exact_t &f,
- exact_t &h) {
- sum_zeroelim(ebegin, eend, &f[0], &f[f.size()], h);
- }
-
- void sum_zeroelim(const exact_t &e,
- const double *fbegin,
- const double *fend,
- exact_t &h) {
- sum_zeroelim(&e[0], &e[e.size()], fbegin, fend, h);
- }
-
-
- exact_t operator+(const exact_t &a, const exact_t &b) {
- exact_t r;
- sum_zeroelim(a, b, r);
- return r;
- }
-
-
-
- void diffprod(const double a, const double b, const double c, const double d, double *r) {
- // return ab - cd;
- double ab[2], cd[2];
- detail::prod_1_1(&a, &b, ab);
- detail::prod_1_1(&c, &d, cd);
- detail::op<2,2>::sub(ab, cd, r);
- }
-
- double orient3dexact(const double *pa,
- const double *pb,
- const double *pc,
- const double *pd) {
- using namespace detail;
-
- double ab[4]; diffprod(pa[0], pb[1], pb[0], pa[1], ab);
- double bc[4]; diffprod(pb[0], pc[1], pc[0], pb[1], bc);
- double cd[4]; diffprod(pc[0], pd[1], pd[0], pc[1], cd);
- double da[4]; diffprod(pd[0], pa[1], pa[0], pd[1], da);
- double ac[4]; diffprod(pa[0], pc[1], pc[0], pa[1], ac);
- double bd[4]; diffprod(pb[0], pd[1], pd[0], pb[1], bd);
-
- exact_t temp;
- exact_t cda, dab, abc, bcd;
- exact_t adet, bdet, cdet, ddet, abdet, cddet, det;
-
- sum_zeroelim(cd, cd + 4, da, da + 4, temp);
- sum_zeroelim(temp, ac, ac + 4, cda);
-
- sum_zeroelim(da, da + 4, ab, ab + 4, temp);
- sum_zeroelim(temp, bd, bd + 4, dab);
-
- negate(bd, bd + 4);
- negate(ac, bd + 4);
-
- sum_zeroelim(ab, ab + 4, bc, bc + 4, temp);
- sum_zeroelim(temp, ac, ac + 4, abc);
-
- sum_zeroelim(bc, bc + 4, cd, cd + 4, temp);
- sum_zeroelim(temp, bd, bd + 4, bcd);
-
- scale_zeroelim(bcd, +pa[2], adet);
- scale_zeroelim(cda, -pb[2], bdet);
- scale_zeroelim(dab, +pc[2], cdet);
- scale_zeroelim(abc, -pd[2], ddet);
-
- sum_zeroelim(adet, bdet, abdet);
- sum_zeroelim(cdet, ddet, cddet);
-
- sum_zeroelim(abdet, cddet, det);
-
- return det[det.size() - 1];
- }
-
- }
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-08 19:11:55 +0300