Unverified Commit 5b240e06 authored by Andrew Goetz's avatar Andrew Goetz
Browse files

Updated uses of 2Sum algorithm to use Hipparchus MathUtils.twoSum.

parent 00d1b70b
......@@ -21,6 +21,8 @@ import java.util.Map;
import org.hipparchus.RealFieldElement;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.MathUtils.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 2Sum for high precision.
final SumAndResidual sumAndResidual = MathUtils.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 2Sum for high precision.
final SumAndResidual sumAndResidual = MathUtils.twoSum(npHigh[i], termValue[i]);
npHigh[i] = sumAndResidual.getSum();
npLow[i] += sumAndResidual.getResidual();
}
}
......
......@@ -21,6 +21,8 @@ import java.util.Date;
import java.util.TimeZone;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.MathUtils.SumAndResidual;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
......@@ -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 2Sum for high precision.
final SumAndResidual sumAndResidual = MathUtils.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 2Sum for high precision.
final SumAndResidual sumAndResidual = MathUtils.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 2Sum 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 = MathUtils.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;
......
......@@ -22,6 +22,9 @@ import java.util.TimeZone;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.MathUtils.FieldSumAndResidual;
import org.hipparchus.util.MathUtils.SumAndResidual;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
......@@ -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 2Sum for high precision.
final FieldSumAndResidual<T> sumAndResidual = MathUtils.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 2Sum for high precision.
final SumAndResidual sumAndResidual = MathUtils.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 2Sum for high precision.
final FieldSumAndResidual<T> sumAndResidual = MathUtils.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 2Sum 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 = MathUtils.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;
......
Markdown is supported
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