Commit f65ebe0a by Andrew Goetz

### Use Hipparchus MathUtils.twoSum for 2Sum algorithm.

parent 00d1b70b
Pipeline #1002 canceled with stage
 ... ... @@ -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 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; ... ...