LSST Applications g1653933729+34a971ddd9,g28da252d5a+6ecf39ab43,g2bbee38e9b+7117419c07,g2bc492864f+7117419c07,g2cdde0e794+704103fe75,g3156d2b45e+6e87dc994a,g347aa1857d+7117419c07,g35bb328faa+34a971ddd9,g3a166c0a6a+7117419c07,g3e281a1b8c+8ec26ec694,g4005a62e65+ba0306790b,g414038480c+40f8eb5d96,g569e0e2b34+d8e684fe9f,g717e5f8c0f+b9f74a9651,g76fe02a360+28e3f83f63,g80478fca09+fdb1f299e2,g82479be7b0+2fffa722b6,g858d7b2824+b9f74a9651,g8b782ad322+b9f74a9651,g8cd86fa7b1+ed49507c07,g90e5cb8612+3a7ace1691,g9125e01d80+34a971ddd9,g979bb04a14+232ddb636d,g9a700aaafa+34a971ddd9,ga5288a1d22+3504ff4bc9,gae0086650b+34a971ddd9,gae74b0b5c6+761ee125be,gb3d89ef091+14ce01cc3a,gb58c049af0+ace264a4f2,gc28159a63d+7117419c07,gcf0d15dbbd+34f38ce56c,gda6a2b7d83+34f38ce56c,gdaeeff99f8+7774323b41,ge33fd446bb+b9f74a9651,ge5cf12406b+d328917b0a,ge79ae78c31+7117419c07,gf0baf85859+890af219f9,gf5289d68f6+fe3ed47e0d,w.2024.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
orientation.cc
Go to the documentation of this file.
1/*
2 * This file is part of sphgeom.
3 *
4 * Developed for the LSST Data Management System.
5 * This product includes software developed by the LSST Project
6 * (http://www.lsst.org).
7 * See the COPYRIGHT file at the top-level directory of this distribution
8 * for details of code ownership.
9 *
10 * This software is dual licensed under the GNU General Public License and also
11 * under a 3-clause BSD license. Recipients may choose which of these licenses
12 * to use; please see the files gpl-3.0.txt and/or bsd_license.txt,
13 * respectively. If you choose the GPL option then the following text applies
14 * (but note that there is still no warranty even if you opt for BSD instead):
15 *
16 * This program is free software: you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation, either version 3 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
29
32
34
35#include <algorithm>
36
38
39
40namespace lsst {
41namespace sphgeom {
42
43namespace {
44
45// `BigFloat` is an exact floating point type.
46struct BigFloat {
47 BigInteger * mantissa;
49
50 BigFloat() : mantissa(0), exponent(0) {}
51 BigFloat(BigInteger * m) : mantissa(m), exponent(0) {}
52};
53
54// `computeProduct` computes the product of 3 doubles exactly and stores the
55// result in p.
56void computeProduct(BigFloat & p, double d0, double d1, double d2) {
57 // This constant (2^53) is used to scale the fractions returned by
58 // std::frexp and turn them into integer mantissas.
59 static double const SCALE = 9007199254740992.0;
60 std::uint32_t buf[2];
61 BigInteger i(buf, sizeof(buf) / sizeof(std::uint32_t));
62 // Unpack the 3 input doubles into integer mantissas and exponents.
63 int e0 = 0;
64 int e1 = 0;
65 int e2 = 0;
66 double m0 = std::frexp(d0, &e0) * SCALE;
67 double m1 = std::frexp(d1, &e1) * SCALE;
68 double m2 = std::frexp(d2, &e2) * SCALE;
69 // Compute the product of the 3 input doubles using exact arithmetic.
70 p.mantissa->setTo(static_cast<std::int64_t>(m0));
71 i.setTo(static_cast<std::int64_t>(m1));
72 p.mantissa->multiply(i);
73 i.setTo(static_cast<std::int64_t>(m2));
74 p.mantissa->multiply(i);
75 // Finally, adjust the exponent of the result to compensate for the 3
76 // multiplications by 2^53 performed above.
77 p.exponent = e0 + e1 + e2 - 3 * 53;
78}
79
80} // unnamed namespace
81
82
84 Vector3d const & b,
85 Vector3d const & c)
86{
87 // Product mantissa storage buffers.
88 std::uint32_t mantissaBuffers[6][6];
89 // Product mantissas.
90 BigInteger mantissas[6] = {
91 BigInteger(mantissaBuffers[0], 6),
92 BigInteger(mantissaBuffers[1], 6),
93 BigInteger(mantissaBuffers[2], 6),
94 BigInteger(mantissaBuffers[3], 6),
95 BigInteger(mantissaBuffers[4], 6),
96 BigInteger(mantissaBuffers[5], 6)
97 };
98 BigFloat products[6] = {
99 BigFloat(&mantissas[0]),
100 BigFloat(&mantissas[1]),
101 BigFloat(&mantissas[2]),
102 BigFloat(&mantissas[3]),
103 BigFloat(&mantissas[4]),
104 BigFloat(&mantissas[5])
105 };
106 // An accumulator and its storage.
107 std::uint32_t accumulatorBuffer[512];
108 BigInteger accumulator(accumulatorBuffer,
109 sizeof(accumulatorBuffer) / sizeof(std::uint32_t));
110 // Compute the products in the determinant. Performing all multiplication
111 // up front means that each product mantissa occupies at most 3*53 bits.
112 computeProduct(products[0], a.x(), b.y(), c.z());
113 computeProduct(products[1], a.x(), b.z(), c.y());
114 computeProduct(products[2], a.y(), b.z(), c.x());
115 computeProduct(products[3], a.y(), b.x(), c.z());
116 computeProduct(products[4], a.z(), b.x(), c.y());
117 computeProduct(products[5], a.z(), b.y(), c.x());
118 mantissas[1].negate();
119 mantissas[3].negate();
120 mantissas[5].negate();
121 // Sort the array of products in descending exponent order.
122 std::sort(products, products + 6, [](BigFloat const & a, BigFloat const & b) {
123 return a.exponent > b.exponent;
124 });
125 // First, initialize the accumulator to the product with the highest
126 // exponent, then add the remaining products. Prior to each addition, we
127 // must shift the accumulated value so that its radix point lines up with
128 // the the radix point of the product to add.
129 //
130 // More precisely, at each step we have an accumulated value A·2ʲ and a
131 // product P·2ᵏ, and we update the accumulator to equal (A·2ʲ⁻ᵏ + P)·2ᵏ.
132 // Because the products were sorted beforehand, j ≥ k and 2ʲ⁻ᵏ is an
133 // integer.
134 accumulator = *products[0].mantissa;
135 for (int i = 1; i < 6; ++i) {
136 accumulator.multiplyPow2(products[i - 1].exponent - products[i].exponent);
137 accumulator.add(*products[i].mantissa);
138 }
139 return accumulator.getSign();
140}
141
143 UnitVector3d const & b,
144 UnitVector3d const & c)
145{
146 // This constant is a little more than 5ε, where ε = 2^-53. When multiplied
147 // by the permanent of |M|, it gives an error bound on the determinant of
148 // M. Here, M is a 3x3 matrix and |M| denotes the matrix obtained by
149 // taking the absolute value of each of its components. The derivation of
150 // this proceeds in the same manner as the derivation of the error bounds
151 // in section 4.3 of:
152 //
153 // Adaptive Precision Floating-Point Arithmetic
154 // and Fast Robust Geometric Predicates,
155 // Jonathan Richard Shewchuk,
156 // Discrete & Computational Geometry 18(3):305–363, October 1997.
157 //
158 // available online at http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf
159 static double const relativeError = 5.6e-16;
160 // Because all 3 unit vectors are normalized, the maximum absolute value of
161 // any vector component, cross product component or dot product term in
162 // the calculation is very close to 1. The permanent of |M| must therefore
163 // be below 3 + c, where c is some small multiple of ε. This constant, a
164 // little larger than 3 * 5ε, is an upper bound on the absolute error in
165 // the determinant calculation.
166 static double const maxAbsoluteError = 1.7e-15;
167 // This constant accounts for floating point underflow (assuming hardware
168 // without gradual underflow, just to be conservative) in the computation
169 // of det(M). It is a little more than 14 * 2^-1022.
170 static double const minAbsoluteError = 4.0e-307;
171
172 double bycz = b.y() * c.z();
173 double bzcy = b.z() * c.y();
174 double bzcx = b.z() * c.x();
175 double bxcz = b.x() * c.z();
176 double bxcy = b.x() * c.y();
177 double bycx = b.y() * c.x();
178 double determinant = a.x() * (bycz - bzcy) +
179 a.y() * (bzcx - bxcz) +
180 a.z() * (bxcy - bycx);
181 if (determinant > maxAbsoluteError) {
182 return 1;
183 } else if (determinant < -maxAbsoluteError) {
184 return -1;
185 }
186 // Expend some more effort on what is hopefully a tighter error bound
187 // before falling back on arbitrary precision arithmetic.
188 double permanent = std::fabs(a.x()) * (std::fabs(bycz) + std::fabs(bzcy)) +
189 std::fabs(a.y()) * (std::fabs(bzcx) + std::fabs(bxcz)) +
190 std::fabs(a.z()) * (std::fabs(bxcy) + std::fabs(bycx));
191 double maxError = relativeError * permanent + minAbsoluteError;
192 if (determinant > maxError) {
193 return 1;
194 } else if (determinant < -maxError) {
195 return -1;
196 }
197 // Avoid the slow path when any two inputs are identical or antipodal.
198 if (a == b || b == c || a == c || a == -b || b == -c || a == -c) {
199 return 0;
200 }
201 return orientationExact(a, b, c);
202}
203
204
205namespace {
206
207 inline int _orientationXYZ(double ab, double ba) {
208 // Calling orientation() with a first argument of (1,0,0), (0,1,0) or
209 // (0,0,1) corresponds to computing the sign of a 2x2 determinant
210 // rather than a 3x3 determinant. The corresponding error bounds
211 // are also tighter.
212 static double const relativeError = 1.12e-16; // > 2^-53
213 static double const maxAbsoluteError = 1.12e-16; // > 2^-53
214 static double const minAbsoluteError = 1.0e-307; // > 3 * 2^-1022
215
216 double determinant = ab - ba;
217 if (determinant > maxAbsoluteError) {
218 return 1;
219 } else if (determinant < -maxAbsoluteError) {
220 return -1;
221 }
222 double permanent = std::fabs(ab) + std::fabs(ba);
223 double maxError = relativeError * permanent + minAbsoluteError;
224 if (determinant > maxError) {
225 return 1;
226 } else if (determinant < -maxError) {
227 return -1;
228 }
229 return 0;
230 }
231
232}
233
234int orientationX(UnitVector3d const & b, UnitVector3d const & c) {
235 int o = _orientationXYZ(b.y() * c.z(), b.z() * c.y());
236 return (o != 0) ? o : orientationExact(UnitVector3d::X(), b, c);
237}
238
239int orientationY(UnitVector3d const & b, UnitVector3d const & c) {
240 int o = _orientationXYZ(b.z() * c.x(), b.x() * c.z());
241 return (o != 0) ? o : orientationExact(UnitVector3d::Y(), b, c);
242}
243
244int orientationZ(UnitVector3d const & b, UnitVector3d const & c) {
245 int o = _orientationXYZ(b.x() * c.y(), b.y() * c.x());
246 return (o != 0) ? o : orientationExact(UnitVector3d::Z(), b, c);
247}
248
249}} // namespace lsst::sphgeom
This file declares a class for arbitrary precision integers.
int m
Definition SpanSet.cc:48
table::Key< int > b
BigInteger is an arbitrary precision signed integer class.
Definition BigInteger.h:51
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.
Definition BigInteger.h:74
void negate()
negate multiplies this integer by -1.
Definition BigInteger.h:109
BigInteger & add(BigInteger const &b)
add adds b to this integer.
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
static UnitVector3d Z()
static UnitVector3d Y()
static UnitVector3d X()
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition Vector3d.h:51
double x() const
Definition Vector3d.h:73
double y() const
Definition Vector3d.h:75
double z() const
Definition Vector3d.h:77
T fabs(T... args)
T frexp(T... args)
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).
T sort(T... args)
This file declares functions for orienting points on the sphere.
BigInteger * mantissa
int exponent