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: e9d68968a5c4
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 343571bb91d0
Choose a head ref
  • 5 commits
  • 2 files changed
  • 1 contributor

Commits on Jun 24, 2019

  1. Rewrite JacobiSNCNDN.

    pleroy committed Jun 24, 2019
    Copy the full SHA
    ddbeba3 View commit details
  2. Remove TODO.

    pleroy committed Jun 24, 2019
    Copy the full SHA
    c689760 View commit details
  3. Some simplification.

    pleroy committed Jun 24, 2019
    Copy the full SHA
    e02eb3a View commit details

Commits on Jun 25, 2019

  1. After egg's review.

    pleroy committed Jun 25, 2019
    Copy the full SHA
    3c18495 View commit details
  2. Merge pull request #2221 from pleroy/SnCnDn

     Rewrite JacobiSNCNDN
    pleroy authored Jun 25, 2019
    Copy the full SHA
    343571b View commit details
Showing with 44 additions and 57 deletions.
  1. +44 −51 numerics/elliptic_functions.cpp
  2. +0 −6 numerics/elliptic_integrals.cpp
95 changes: 44 additions & 51 deletions numerics/elliptic_functions.cpp
Original file line number Diff line number Diff line change
@@ -13,6 +13,7 @@

namespace principia {

using quantities::Abs;
using quantities::Sqrt;

namespace numerics {
@@ -127,8 +128,6 @@ void JacobiSNCNDNReduced(double const u,
//
// Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
//
// Used subprograms: scd2, elk
//
// Inputs: u = argument, mc = 1-m, 0 < mc <= 1
//
// Output: s = sn(u|m), c=cn(u|m), d=dn(u|m)
@@ -138,65 +137,59 @@ void JacobiSNCNDN(double const u,
double& s,
double& c,
double& d) {
double m, kc, ux, k, kh, kh3, kh5, kh7, k2, k3, k4, sx;
constexpr double k_over_2_lower_bound = π / 4.0;

m = 1.0 - mc;
kc = sqrt(mc);
ux = abs(u);
if (ux < 0.785) {
JacobiSNCNDNReduced(ux, mc, s, c, d);
// The argument reduction follows Fukushima (2009), Fast computation of
// Jacobian elliptic function and incomplete elliptic integrals for constant
// values of elliptic parameter and elliptic characteristic, sections 2.4 and
// 3.5.2.
double const m = 1.0 - mc;
double const kʹ = Sqrt(mc);
double abs_u = Abs(u);
if (abs_u < k_over_2_lower_bound) {
JacobiSNCNDNReduced(abs_u, mc, s, c, d);
} else {
k = EllipticK(mc);
kh = k * 0.5;
kh3 = k * 1.5;
kh5 = k * 2.5;
kh7 = k * 3.5;
k2 = k * 2.0;
k3 = k * 3.0;
k4 = k * 4.0;
ux = ux - k4 * static_cast<double>(static_cast<int>(ux / k4));
if (ux < kh) {
JacobiSNCNDNReduced(ux, mc, s, c, d);
} else if (ux < k) {
ux = k - ux;
JacobiSNCNDNReduced(ux, mc, s, c, d);
sx = c / d;
c = kc * s / d;
double const k = EllipticK(mc);
double const two_k = 2.0 * k;
double const three_k = 3.0 * k;
double const four_k = 4.0 * k;
abs_u =
abs_u - four_k * static_cast<double>(static_cast<int>(abs_u / four_k));
if (abs_u < 0.5 * k) {
JacobiSNCNDNReduced(abs_u, mc, s, c, d);
} else if (abs_u < k) {
JacobiSNCNDNReduced(k - abs_u, mc, s, c, d);
double const sx = c / d;
c = kʹ * s / d;
s = sx;
d = kc / d;
} else if (ux < kh3) {
ux = ux - k;
JacobiSNCNDNReduced(ux, mc, s, c, d);
sx = c / d;
c = -kc * s / d;
d = kʹ / d;
} else if (abs_u < 1.5 * k) {
JacobiSNCNDNReduced(abs_u - k, mc, s, c, d);
double const sx = c / d;
c = -kʹ * s / d;
s = sx;
d = kc / d;
} else if (ux < k2) {
ux = k2 - ux;
JacobiSNCNDNReduced(ux, mc, s, c, d);
d = kʹ / d;
} else if (abs_u < two_k) {
JacobiSNCNDNReduced(two_k - abs_u, mc, s, c, d);
c = -c;
} else if (ux < kh5) {
ux = ux - k2;
JacobiSNCNDNReduced(ux, mc, s, c, d);
} else if (abs_u < 2.5 * k) {
JacobiSNCNDNReduced(abs_u - two_k, mc, s, c, d);
s = -s;
c = -c;
} else if (ux < k3) {
ux = k3 - ux;
JacobiSNCNDNReduced(ux, mc, s, c, d);
sx = -c / d;
c = -kc * s / d;
} else if (abs_u < three_k) {
JacobiSNCNDNReduced(three_k - abs_u, mc, s, c, d);
double const sx = -c / d;
c = -kʹ * s / d;
s = sx;
d = kc / d;
} else if (ux < kh7) {
ux = ux - k3;
JacobiSNCNDNReduced(ux, mc, s, c, d);
sx = -c / d;
c = kc * s / d;
d = kʹ / d;
} else if (abs_u < 3.5 * k) {
JacobiSNCNDNReduced(abs_u - three_k, mc, s, c, d);
double const sx = -c / d;
c = kʹ * s / d;
s = sx;
d = kc / d;
d = / d;
} else {
ux = k4 - ux;
JacobiSNCNDNReduced(ux, mc, s, c, d);
JacobiSNCNDNReduced(four_k - abs_u, mc, s, c, d);
s = -s;
}
}
6 changes: 0 additions & 6 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -14,12 +14,6 @@
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"

// TODO(phl):
// 1. Use arrays for the coefficients.
// 2. Use Estrin evaluation for polynomials of high degree (possibly adding
// support for polynomials of two and three variables).
// 3. Figure something for the uninitialized variables.

// Bibliography:
// [Buli69] Bulirsch (1969), Numerical Calculation of Elliptic Integrals and
// Elliptic Fuctions. III.