49 void computeProduct(BigFloat & p,
double d0,
double d1,
double d2) {
52 static double const SCALE = 9007199254740992.0;
54 BigInteger i(buf,
sizeof(buf) /
sizeof(uint32_t));
63 p.mantissa->setTo(
static_cast<int64_t
>(m0));
64 i.setTo(
static_cast<int64_t
>(m1));
65 p.mantissa->multiply(i);
66 i.setTo(
static_cast<int64_t
>(m2));
67 p.mantissa->multiply(i);
70 p.exponent = e0 + e1 + e2 - 3 * 53;
81 uint32_t mantissaBuffers[6][6];
91 BigFloat products[6] = {
92 BigFloat(&mantissas[0]),
93 BigFloat(&mantissas[1]),
94 BigFloat(&mantissas[2]),
95 BigFloat(&mantissas[3]),
96 BigFloat(&mantissas[4]),
97 BigFloat(&mantissas[5])
100 uint32_t accumulatorBuffer[512];
102 sizeof(accumulatorBuffer) /
sizeof(uint32_t));
105 computeProduct(products[0],
a.x(),
b.y(), c.
z());
106 computeProduct(products[1],
a.x(),
b.z(), c.
y());
107 computeProduct(products[2],
a.y(),
b.z(), c.
x());
108 computeProduct(products[3],
a.y(),
b.x(), c.
z());
109 computeProduct(products[4],
a.z(),
b.x(), c.
y());
110 computeProduct(products[5],
a.z(),
b.y(), c.
x());
115 std::sort(products, products + 6, [](BigFloat
const &
a, BigFloat
const &
b) {
116 return a.exponent >
b.exponent;
127 accumulator = *products[0].mantissa;
128 for (
int i = 1; i < 6; ++i) {
152 static double const relativeError = 5.6e-16;
159 static double const maxAbsoluteError = 1.7e-15;
163 static double const minAbsoluteError = 4.0e-307;
165 double bycz =
b.y() * c.
z();
166 double bzcy =
b.z() * c.
y();
167 double bzcx =
b.z() * c.
x();
168 double bxcz =
b.x() * c.
z();
169 double bxcy =
b.x() * c.
y();
170 double bycx =
b.y() * c.
x();
171 double determinant =
a.x() * (bycz - bzcy) +
172 a.y() * (bzcx - bxcz) +
173 a.z() * (bxcy - bycx);
174 if (determinant > maxAbsoluteError) {
176 }
else if (determinant < -maxAbsoluteError) {
184 double maxError = relativeError * permanent + minAbsoluteError;
185 if (determinant > maxError) {
187 }
else if (determinant < -maxError) {
191 if (
a ==
b ||
b == c ||
a == c ||
a == -
b ||
b == -c ||
a == -c) {
200 inline int _orientationXYZ(
double ab,
double ba) {
205 static double const relativeError = 1.12e-16;
206 static double const maxAbsoluteError = 1.12e-16;
207 static double const minAbsoluteError = 1.0e-307;
209 double determinant = ab - ba;
210 if (determinant > maxAbsoluteError) {
212 }
else if (determinant < -maxAbsoluteError) {
216 double maxError = relativeError * permanent + minAbsoluteError;
217 if (determinant > maxError) {
219 }
else if (determinant < -maxError) {
228 int o = _orientationXYZ(
b.y() * c.
z(),
b.z() * c.
y());
233 int o = _orientationXYZ(
b.z() * c.
x(),
b.x() * c.
z());
238 int o = _orientationXYZ(
b.x() * c.
y(),
b.y() * c.
x());
This file declares a class for arbitrary precision integers.
BigInteger is an arbitrary precision signed integer class.
BigInteger & multiplyPow2(unsigned n)
multiplyPow2 multiplies this integer by 2ⁿ.
int getSign() const
getSign returns -1, 0 or 1 if this integer is negative, zero or positive.
void negate()
negate multiplies this integer by -1.
BigInteger & add(BigInteger const &b)
add adds b to this integer.
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Vector3d is a vector in ℝ³ with components stored in double precision.
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
int orientation(UnitVector3d const &a, UnitVector3d const &b, UnitVector3d const &c)
orientation computes and returns the orientations of 3 unit vectors a, b and c.
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
int orientationExact(Vector3d const &a, Vector3d const &b, Vector3d const &c)
orientationExact computes and returns the orientations of 3 vectors a, b and c, which need not be nor...
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
A base class for image defects.
This file declares functions for orienting points on the sphere.