Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 4a1eeff37e08
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: e4dad90d9531
Choose a head ref
  • 4 commits
  • 1 file changed
  • 1 contributor

Commits on Jun 18, 2019

  1. Copy the full SHA
    8e241e6 View commit details
  2. Rewrite FukushimaT.

    pleroy committed Jun 18, 2019
    Copy the full SHA
    996b489 View commit details
  3. After egg's review.

    pleroy committed Jun 18, 2019
    Copy the full SHA
    ae2175b View commit details
  4. Merge pull request #2212 from pleroy/FunctionRewrite

     Rewrite FukushimaT
    pleroy authored Jun 18, 2019

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    e4dad90 View commit details
Showing with 83 additions and 110 deletions.
  1. +83 −110 numerics/elliptic_integrals.cpp
193 changes: 83 additions & 110 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -2,8 +2,12 @@
#include "glog/logging.h"

#include <limits>
#include <tuple>
#include <utility>

#include "numerics/elliptic_integrals.hpp"
#include "numerics/polynomial.hpp"
#include "numerics/polynomial_evaluators.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/numbers.hpp"
#include "quantities/quantities.hpp"
@@ -88,6 +92,47 @@ double FukushimaEllipticJsMaclaurinSeries(double y, double n, double m);
// Fukushima's T function [Fuku11c].
double FukushimaT(double t, double h);

// A generator for the Maclaurin series for Fukushima's T function.
template<int n>
class FukushimaTMaclaurin {
template<typename>
struct Generator;

template<int... k>
struct Generator<std::index_sequence<k...>> {
static auto constexpr series = std::make_tuple(1.0 / (2.0 * k + 1.0)...);
};

public:
static constexpr auto series =
Generator<std::make_index_sequence<n + 1>>::series;
};

PolynomialInMonomialBasis<double, double, 1, HornerEvaluator>
fukushima_t_maclaurin_1(FukushimaTMaclaurin<1>::series);
PolynomialInMonomialBasis<double, double, 2, HornerEvaluator>
fukushima_t_maclaurin_2(FukushimaTMaclaurin<2>::series);
PolynomialInMonomialBasis<double, double, 3, HornerEvaluator>
fukushima_t_maclaurin_3(FukushimaTMaclaurin<3>::series);
PolynomialInMonomialBasis<double, double, 4, EstrinEvaluator>
fukushima_t_maclaurin_4(FukushimaTMaclaurin<4>::series);
PolynomialInMonomialBasis<double, double, 5, EstrinEvaluator>
fukushima_t_maclaurin_5(FukushimaTMaclaurin<5>::series);
PolynomialInMonomialBasis<double, double, 6, EstrinEvaluator>
fukushima_t_maclaurin_6(FukushimaTMaclaurin<6>::series);
PolynomialInMonomialBasis<double, double, 7, EstrinEvaluator>
fukushima_t_maclaurin_7(FukushimaTMaclaurin<7>::series);
PolynomialInMonomialBasis<double, double, 8, EstrinEvaluator>
fukushima_t_maclaurin_8(FukushimaTMaclaurin<8>::series);
PolynomialInMonomialBasis<double, double, 9, EstrinEvaluator>
fukushima_t_maclaurin_9(FukushimaTMaclaurin<9>::series);
PolynomialInMonomialBasis<double, double, 10, EstrinEvaluator>
fukushima_t_maclaurin_10(FukushimaTMaclaurin<10>::series);
PolynomialInMonomialBasis<double, double, 11, EstrinEvaluator>
fukushima_t_maclaurin_11(FukushimaTMaclaurin<11>::series);
PolynomialInMonomialBasis<double, double, 12, EstrinEvaluator>
fukushima_t_maclaurin_12(FukushimaTMaclaurin<12>::series);

// Double precision general complete elliptic integral "cel"
//
// created by Burlisch
@@ -1244,119 +1289,47 @@ double FukushimaEllipticJsMaclaurinSeries(double const y,
}

double FukushimaT(double const t, double const h) {
double z, y, x, a, r, ri;
constexpr double A3 = 1.0 / 3.0;
constexpr double A5 = 1.0 / 5.0;
constexpr double A7 = 1.0 / 7.0;
constexpr double A9 = 1.0 / 9.0;
constexpr double A11 = 1.0 / 11.0;
constexpr double A13 = 1.0 / 13.0;
constexpr double A15 = 1.0 / 15.0;
constexpr double A17 = 1.0 / 17.0;
constexpr double A19 = 1.0 / 19.0;
constexpr double A21 = 1.0 / 21.0;
constexpr double A23 = 1.0 / 23.0;
constexpr double A25 = 1.0 / 25.0;

z = -h * t * t;
a = abs(z);

if (a < 3.3306691e-16) {
double const z = -h * t * t;
double const abs_z = abs(z);

// NOTE(phl): One might be tempted to rewrite this statement using a binary
// split of the interval [0, 1], but according to Table 1 of [Fuku11c] the
// distribution of z is very biased towards the small values, so this is
// simpler and probably better. (It also explains the position of z < 0 in
// the list.)
if (abs_z < 3.3306691e-16) {
return t;
} else if (a < 2.3560805e-08) {
return t * (1.0 + z * A3);
} else if (a < 9.1939631e-06) {
return t * (1.0 + z * (A3 + z * A5));
} else if (a < 1.7779240e-04) {
return t * (1.0 + z * (A3 + z * (A5 + z * A7)));
} else if (a < 1.0407839e-03) {
return t * (1.0 + z * (A3 + z * (A5 + z * (A7 + z * A9))));
} else if (a < 3.3616998e-03) {
return t * (1.0 + z * (A3 + z * (A5 + z * (A7 + z * (A9 + z * A11)))));
} else if (a < 7.7408014e-03) {
return
t * (1.0 +
z * (A3 + z * (A5 + z * (A7 + z * (A9 + z * (A11 + z * A13))))));
} else if (a < 1.4437181e-02) {
return
t * (1.0 +
z * (A3 +
z * (A5 +
z * (A7 + z * (A9 + z * (A11 + z * (A13 + z * A15)))))));
} else if (a < 2.3407312e-02) {
return
t *
(1.0 +
z * (A3 +
z * (A5 +
z * (A7 +
z * (A9 +
z * (A11 + z * (A13 + z * (A15 + z * A17))))))));
} else if (a < 3.4416203e-02) {
return
t *
(1.0 +
z * (A3 +
z * (A5 +
z * (A7 +
z * (A9 +
z * (A11 +
z * (A13 +
z * (A15 + z * (A17 + z * A19)))))))));
} else if (abs_z < 2.3560805e-08) {
return t * fukushima_t_maclaurin_1.Evaluate(z);
} else if (abs_z < 9.1939631e-06) {
return t * fukushima_t_maclaurin_2.Evaluate(z);
} else if (abs_z < 1.7779240e-04) {
return t * fukushima_t_maclaurin_3.Evaluate(z);
} else if (abs_z < 1.0407839e-03) {
return t * fukushima_t_maclaurin_4.Evaluate(z);
} else if (abs_z < 3.3616998e-03) {
return t * fukushima_t_maclaurin_5.Evaluate(z);
} else if (abs_z < 7.7408014e-03) {
return t * fukushima_t_maclaurin_6.Evaluate(z);
} else if (abs_z < 1.4437181e-02) {
return t * fukushima_t_maclaurin_7.Evaluate(z);
} else if (abs_z < 2.3407312e-02) {
return t * fukushima_t_maclaurin_8.Evaluate(z);
} else if (abs_z < 3.4416203e-02) {
return t * fukushima_t_maclaurin_9.Evaluate(z);
} else if (z < 0.0) {
r = sqrt(h);
ri = 1.0 / r;
return atan(r * t) * ri;
} else if (a < 4.7138547e-02) {
return
t *
(1.0 +
z * (A3 +
z * (A5 +
z * (A7 +
z * (A9 +
z * (A11 +
z * (A13 +
z * (A15 +
z * (A17 +
z * (A19 + z * A21))))))))));
} else if (a < 6.1227405e-02) {
return
t *
(1.0 +
z * (A3 +
z * (A5 +
z * (A7 +
z * (A9 +
z * (A11 +
z * (A13 +
z * (A15 +
z * (A17 +
z * (A19 +
z * (A21 +
z * A23)))))))))));
} else if (a < 7.6353468e-02) {
return
t *
(1.0 +
z * (A3 +
z * (A5 +
z * (A7 +
z * (A9 +
z * (A11 +
z * (A13 +
z * (A15 +
z * (A17 +
z * (A19 +
z * (A21 +
z * (A23 +
z * A25))))))))))));
double const r = Sqrt(h);
double const ri = 1.0 / r;
return std::atan(r * t) / r;
} else if (abs_z < 4.7138547e-02) {
return t * fukushima_t_maclaurin_10.Evaluate(z);
} else if (abs_z < 6.1227405e-02) {
return t * fukushima_t_maclaurin_11.Evaluate(z);
} else if (abs_z < 7.6353468e-02) {
return t * fukushima_t_maclaurin_12.Evaluate(z);
} else {
r = sqrt(-h);
ri = 1.0 / r;
y = r * t;
x = log((1.0 + y) / (1.0 - y)) * 0.5;
return x * ri;
double const r = Sqrt(-h);
return std::atanh(r * t) / r;
}
}