Skip to content
GitLab
Unverified Commit f125ab9a authored by Andrew Goetz's avatar Andrew Goetz
Browse files

Added dedicated class in org.orekit.utils for TwoSum algorithm. 

Refactored existing instances of TwoSum algorithm to use class.
parent 46b0e5c5
Pipeline #996 failed with stage
Hide whitespace changes
Inline Side-by-side
src/main/java/org/orekit/data/PoissonSeries.java
View file @ f125ab9a
......@@ -21,6 +21,8 @@ import java.util.Map;
import org.hipparchus.RealFieldElement;
import org.hipparchus.util.MathArrays;
import org.orekit.utils.TwoSum;
import org.orekit.utils.TwoSum.SumAndResidual;
/**
* Class representing a Poisson series for nutation or ephemeris computations.
......@@ -77,20 +79,14 @@ public class PoissonSeries {
final double p = polynomial.value(elements.getTC());
// non-polynomial part
// compute sum accurately, using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
double npHigh = 0;
double npLow = 0;
for (final Map.Entry<Long, SeriesTerm> entry : series.entrySet()) {
final double v = entry.getValue().value(elements)[0];
final double sum = npHigh + v;
final double sPrime = sum - v;
final double tPrime = sum - sPrime;
final double deltaS = npHigh - sPrime;
final double deltaT = v - tPrime;
npLow += deltaS + deltaT;
npHigh = sum;
// Use TwoSum algorithm for high precision.
final SumAndResidual sumAndResidual = TwoSum.twoSum(npHigh, v);
npHigh = sumAndResidual.getSum();
npLow += sumAndResidual.getResidual();
}
// add the polynomial and the non-polynomial parts
......@@ -220,22 +216,15 @@ public class PoissonSeries {
public double[] value(final BodiesElements elements) {
// non-polynomial part
// compute sum accurately, using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
final double[] npHigh = new double[polynomials.length];
final double[] npLow = new double[polynomials.length];
for (final SeriesTerm term : joinedTerms) {
final double[] termValue = term.value(elements);
for (int i = 0; i < termValue.length; ++i) {
final double v = termValue[i];
final double sum = npHigh[i] + v;
final double sPrime = sum - v;
final double tPrime = sum - sPrime;
final double deltaS = npHigh[i] - sPrime;
final double deltaT = v - tPrime;
npLow[i] += deltaS + deltaT;
npHigh[i] = sum;
// Use TwoSum algorithm for high precision.
final SumAndResidual sumAndResidual = TwoSum.twoSum(npHigh[i], termValue[i]);
npHigh[i] = sumAndResidual.getSum();
npLow[i] += sumAndResidual.getResidual();
}
}
......
src/main/java/org/orekit/time/AbsoluteDate.java
View file @ f125ab9a
......@@ -27,6 +27,8 @@ import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitIllegalArgumentException;
import org.orekit.errors.OrekitMessages;
import org.orekit.utils.Constants;
import org.orekit.utils.TwoSum;
import org.orekit.utils.TwoSum.SumAndResidual;
/** This class represents a specific instant in time.
......@@ -309,21 +311,11 @@ public class AbsoluteDate
final double seconds = time.getSecond();
final double tsOffset = timeScale.offsetToTAI(date, time);
// compute sum exactly, using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
// at the end, the EXACT result of addition seconds + tsOffset
// is sum + residual, where sum is the closest representable number to the exact
// result and residual is the missing part that does not fit in the first number
final double sum = seconds + tsOffset;
final double sPrime = sum - tsOffset;
final double tPrime = sum - sPrime;
final double deltaS = seconds - sPrime;
final double deltaT = tsOffset - tPrime;
final double residual = deltaS + deltaT;
final long dl = (long) FastMath.floor(sum);
offset = (sum - dl) + residual;
// Use TwoSum for high precision.
final SumAndResidual sumAndResidual = TwoSum.twoSum(seconds, tsOffset);
final long dl = (long) FastMath.floor(sumAndResidual.getSum());
offset = (sumAndResidual.getSum() - dl) + sumAndResidual.getResidual();
epoch = 60l * ((date.getJ2000Day() * 24l + time.getHour()) * 60l +
time.getMinute() - time.getMinutesFromUTC() - 720l) + dl;
......@@ -430,26 +422,15 @@ public class AbsoluteDate
* @see #durationFrom(AbsoluteDate)
*/
public AbsoluteDate(final AbsoluteDate since, final double elapsedDuration) {
final double sum = since.offset + elapsedDuration;
if (Double.isInfinite(sum)) {
offset = sum;
epoch = (sum < 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
// Use TwoSum for high precision.
final SumAndResidual sumAndResidual = TwoSum.twoSum(since.offset, elapsedDuration);
if (Double.isInfinite(sumAndResidual.getSum())) {
offset = sumAndResidual.getSum();
epoch = (sumAndResidual.getSum() < 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
} else {
// compute sum exactly, using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
// at the end, the EXACT result of addition since.offset + elapsedDuration
// is sum + residual, where sum is the closest representable number to the exact
// result and residual is the missing part that does not fit in the first number
final double oPrime = sum - elapsedDuration;
final double dPrime = sum - oPrime;
final double deltaO = since.offset - oPrime;
final double deltaD = elapsedDuration - dPrime;
final double residual = deltaO + deltaD;
final long dl = (long) FastMath.floor(sum);
offset = (sum - dl) + residual;
epoch = since.epoch + dl;
final long dl = (long) FastMath.floor(sumAndResidual.getSum());
offset = (sumAndResidual.getSum() - dl) + sumAndResidual.getResidual();
epoch = since.epoch + dl;
}
}
......@@ -1057,24 +1038,14 @@ public class AbsoluteDate
}
}
// compute offset from 2000-01-01T00:00:00 in specified time scale exactly,
// using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
// at the end, the EXACT result of addition offset + timeScale.offsetFromTAI(this)
// is sum + residual, where sum is the closest representable number to the exact
// result and residual is the missing part that does not fit in the first number
// Compute offset from 2000-01-01T00:00:00 in specified time scale.
// Use TwoSum for high precision.
final double taiOffset = timeScale.offsetFromTAI(this);
final double sum = offset + taiOffset;
final double oPrime = sum - taiOffset;
final double dPrime = sum - oPrime;
final double deltaO = offset - oPrime;
final double deltaD = taiOffset - dPrime;
final double residual = deltaO + deltaD;
final SumAndResidual sumAndResidual = TwoSum.twoSum(offset, taiOffset);
// split date and time
final long carry = (long) FastMath.floor(sum);
double offset2000B = (sum - carry) + residual;
final long carry = (long) FastMath.floor(sumAndResidual.getSum());
double offset2000B = (sumAndResidual.getSum() - carry) + sumAndResidual.getResidual();
long offset2000A = epoch + carry + 43200l;
if (offset2000B < 0) {
offset2000A -= 1;
......
src/main/java/org/orekit/time/FieldAbsoluteDate.java
View file @ f125ab9a
......@@ -27,6 +27,9 @@ import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.utils.Constants;
import org.orekit.utils.TwoSum;
import org.orekit.utils.TwoSum.FieldSumAndResidual;
import org.orekit.utils.TwoSum.SumAndResidual;
/** This class represents a specific instant in time.
......@@ -151,25 +154,15 @@ public class FieldAbsoluteDate<T extends RealFieldElement<T>>
*/
public FieldAbsoluteDate(final FieldAbsoluteDate<T> since, final T elapsedDuration) {
this.field = since.field;
final T sum = since.offset.add(elapsedDuration);
if (Double.isInfinite(sum.getReal())) {
offset = sum;
epoch = (sum.getReal() < 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
// Use TwoSum for high precision.
final FieldSumAndResidual<T> sumAndResidual = TwoSum.twoSum(since.offset, elapsedDuration);
if (Double.isInfinite(sumAndResidual.getSum().getReal())) {
offset = sumAndResidual.getSum();
epoch = (sumAndResidual.getSum().getReal() < 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
} else {
// compute sum exactly, using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
// at the end, the EXACT result of addition since.offset + elapsedDuration
// is sum + residual, where sum is the closest representable number to the exact
// result and residual is the missing part that does not fit in the first number
final double oPrime = sum.getReal() - elapsedDuration.getReal();
final double dPrime = sum.getReal() - oPrime;
final double deltaO = since.offset.getReal() - oPrime;
final double deltaD = elapsedDuration.getReal() - dPrime;
final double residual = deltaO + deltaD;
final long dl = (long) FastMath.floor(sum.getReal());
offset = sum.subtract(dl).add(residual);
epoch = since.epoch + dl;
final long dl = (long) FastMath.floor(sumAndResidual.getSum().getReal());
offset = sumAndResidual.getSum().subtract(dl).add(sumAndResidual.getResidual());
epoch = since.epoch + dl;
}
}
......@@ -212,21 +205,11 @@ public class FieldAbsoluteDate<T extends RealFieldElement<T>>
final double seconds = time.getSecond();
final double tsOffset = timeScale.offsetToTAI(date, time);
// compute sum exactly, using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
// at the end, the EXACT result of addition seconds + tsOffset
// is sum + residual, where sum is the closest representable number to the exact
// result and residual is the missing part that does not fit in the first number
final double sum = seconds + tsOffset;
final double sPrime = sum - tsOffset;
final double tPrime = sum - sPrime;
final double deltaS = seconds - sPrime;
final double deltaT = tsOffset - tPrime;
final double residual = deltaS + deltaT;
final long dl = (long) FastMath.floor(sum);
offset = field.getZero().add((sum - dl) + residual);
// Use TwoSum for high precision.
final SumAndResidual sumAndResidual = TwoSum.twoSum(seconds, tsOffset);
final long dl = (long) FastMath.floor(sumAndResidual.getSum());
offset = field.getZero().add((sumAndResidual.getSum() - dl) + sumAndResidual.getResidual());
epoch = 60l * ((date.getJ2000Day() * 24l + time.getHour()) * 60l +
time.getMinute() - time.getMinutesFromUTC() - 720l) + dl;
......@@ -379,25 +362,15 @@ public class FieldAbsoluteDate<T extends RealFieldElement<T>>
*/
private FieldAbsoluteDate(final long epoch, final double tA, final T tB) {
this.field = tB.getField();
final T sum = tB.add(tA);
if (Double.isInfinite(sum.getReal())) {
this.offset = sum;
this.epoch = (sum.getReal() < 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
// Use TwoSum for high precision.
final FieldSumAndResidual<T> sumAndResidual = TwoSum.twoSum(field.getZero().add(tA), tB);
if (Double.isInfinite(sumAndResidual.getSum().getReal())) {
this.offset = sumAndResidual.getSum();
this.epoch = (sumAndResidual.getSum().getReal() < 0) ? Long.MIN_VALUE : Long.MAX_VALUE;
} else {
// compute sum exactly, using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
// at the end, the EXACT result of addition tA + tB
// is sum + residual, where sum is the closest representable number to the exact
// result and residual is the missing part that does not fit in the first number
final double oPrime = sum.getReal() - tA;
final double dPrime = sum.getReal() - oPrime;
final double deltaO = tB.getReal() - oPrime;
final double deltaD = tA - dPrime;
final double residual = deltaO + deltaD;
final long dl = (long) FastMath.floor(sum.getReal());
this.offset = sum.subtract(dl).add(residual);
this.epoch = epoch + dl;
final long dl = (long) FastMath.floor(sumAndResidual.getSum().getReal());
this.offset = sumAndResidual.getSum().subtract(dl).add(sumAndResidual.getResidual());
this.epoch = epoch + dl;
}
}
......@@ -1239,24 +1212,14 @@ public class FieldAbsoluteDate<T extends RealFieldElement<T>>
}
}
// compute offset from 2000-01-01T00:00:00 in specified time scale exactly,
// using Møller-Knuth TwoSum algorithm without branching
// the following statements must NOT be simplified, they rely on floating point
// arithmetic properties (rounding and representable numbers)
// at the end, the EXACT result of addition offset + timeScale.offsetFromTAI(this)
// is sum + residual, where sum is the closest representable number to the exact
// result and residual is the missing part that does not fit in the first number
// Compute offset from 2000-01-01T00:00:00 in specified time scale.
// Use TwoSum for high accuracy.
final double taiOffset = timeScale.offsetFromTAI(this).getReal();
final double sum = offset.getReal() + taiOffset;
final double oPrime = sum - taiOffset;
final double dPrime = sum - oPrime;
final double deltaO = offset.getReal() - oPrime;
final double deltaD = taiOffset - dPrime;
final double residual = deltaO + deltaD;
final SumAndResidual sumAndResidual = TwoSum.twoSum(offset.getReal(), taiOffset);
// split date and time
final long carry = (long) FastMath.floor(sum);
double offset2000B = (sum - carry) + residual;
final long carry = (long) FastMath.floor(sumAndResidual.getSum());
double offset2000B = (sumAndResidual.getSum() - carry) + sumAndResidual.getResidual();
long offset2000A = epoch + carry + 43200l;
if (offset2000B < 0) {
offset2000A -= 1;
......
src/main/java/org/orekit/utils/TwoSum.java 0 → 100644
View file @ f125ab9a
/* Copyright 2002-2021 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.utils;
import org.hipparchus.RealFieldElement;
/**
* Møller-Knuth TwoSum algorithm for adding two floating-point numbers while tracking the residual
* error.
*/
public final class TwoSum {
private TwoSum() {
}
/**
* Sums {@code a} and {@code b} using the TwoSum algorithm.
* @param a first summand
* @param b second summand
* @return sum and residual error
*/
public static SumAndResidual twoSum(final double a, final double b) {
final double s = a + b;
final double aPrime = s - b;
final double bPrime = s - aPrime;
final double deltaA = a - aPrime;
final double deltaB = b - bPrime;
final double t = deltaA + deltaB;
return new SumAndResidual(s, t);
}
/**
* Sums {@code a} and {@code b} using the TwoSum algorithm.
* @param <T> field element type
* @param a first summand
* @param b second summand
* @return sum and residual error
*/
public static <T extends RealFieldElement<T>> FieldSumAndResidual<T> twoSum(final T a, final T b) {
final T s = a.add(b);
final T aPrime = s.subtract(b);
final T bPrime = s.subtract(aPrime);
final T deltaA = a.subtract(aPrime);
final T deltaB = b.subtract(bPrime);
final T t = deltaA.add(deltaB);
return new FieldSumAndResidual<>(s, t);
}
/**
* Result class containing the sum and the residual error in the sum.
*/
public static final class SumAndResidual {
/** Sum. */
private final double sum;
/** Residual error. */
private final double residual;
/**
* Constructs a {@link SumAndResidual} instance.
* @param sum sum
* @param residual residual
*/
private SumAndResidual(final double sum, final double residual) {
this.sum = sum;
this.residual = residual;
}
/**
* Returns the sum.
* @return sum
*/
public double getSum() {
return sum;
}
/**
* Returns the residual error.
* @return residual error
*/
public double getResidual() {
return residual;
}
}
/**
* Result class containing the sum and the residual error in the sum.
* @param <T> field element type
*/
public static final class FieldSumAndResidual<T extends RealFieldElement<T>> {
/** Sum. */
private final T sum;
/** Residual error. */
private final T residual;
/**
* Constructs a {@link SumAndResidual} instance.
* @param sum sum
* @param residual residual
*/
private FieldSumAndResidual(final T sum, final T residual) {
this.sum = sum;
this.residual = residual;
}
/**
* Returns the sum.
* @return sum
*/
public T getSum() {
return sum;
}
/**
* Returns the residual error.
* @return residual error
*/
public T getResidual() {
return residual;
}
}
}
src/test/java/org/orekit/utils/TwoSumTest.java 0 → 100644
View file @ f125ab9a
package org.orekit.utils;
import org.hipparchus.RealFieldElement;
import org.hipparchus.util.Tuple;
import org.junit.Assert;
import org.junit.Test;
import org.orekit.utils.TwoSum.FieldSumAndResidual;
import org.orekit.utils.TwoSum.SumAndResidual;
/**
* Unit tests for {@link TwoSum}.
*/
public final class TwoSumTest {
/**
* Tests {@link TwoSum#twoSum(double, double)}.
*/
@Test
public void testTwoSum() {
// sum = 0.30000000000000004
// residual = -2.7755575615628914E-17
final double a1 = 0.1;
final double b1 = 0.2;
final SumAndResidual result1 = TwoSum.twoSum(a1, b1);
Assert.assertEquals(a1 + b1, result1.getSum(), 0.);
Assert.assertEquals(a1 + b1, result1.getSum() + result1.getResidual(), 0.);
// sum = -1580.3205849419005
// residual = -1.1368683772161603E-13
final double a2 = -615.7212034581913;
final double b2 = -964.5993814837093;
final SumAndResidual result2 = TwoSum.twoSum(a2, b2);
Assert.assertEquals(a2 + b2, result2.getSum(), 0.);
Assert.assertEquals(a2 + b2, result2.getSum() + result2.getResidual(), 0.);
// sum = 251.8625825973395
// residual = 1.4210854715202004E-14
final double a3 = 60.348375484313706;
final double b3 = 191.5142071130258;
final SumAndResidual result3 = TwoSum.twoSum(a3, b3);
Assert.assertEquals(a3 + b3, result3.getSum(), 0.);
Assert.assertEquals(a3 + b3, result3.getSum() + result3.getResidual(), 0.);
// sum = 622.8319023175123
// residual = -4.3315557304163255E-14
final double a4 = 622.8314146170453;
final double b4 = 0.0004877004669900762;
final SumAndResidual result4 = TwoSum.twoSum(a4, b4);
Assert.assertEquals(a4 + b4, result4.getSum(), 0.);
Assert.assertEquals(a4 + b4, result4.getSum() + result4.getResidual(), 0.);
}
/**
* Tests {@link TwoSum#twoSum(RealFieldElement, RealFieldElement)}.
*/
@Test
public void testTwoSumField() {
final Tuple a = new Tuple(0.1, -615.7212034581913, 60.348375484313706, 622.8314146170453);
final Tuple b = new Tuple(0.2, -964.5993814837093, 191.5142071130258, 0.0004877004669900762);
final FieldSumAndResidual<Tuple> result = TwoSum.twoSum(a, b);
for (int i = 0; i < a.getDimension(); ++i) {
Assert.assertEquals(a.getComponent(i) + b.getComponent(i), result.getSum().getComponent(i), 0.);
Assert.assertEquals(a.getComponent(i) + b.getComponent(i),
result.getSum().getComponent(i) + result.getResidual().getComponent(i), 0.);
}
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment