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

Commits on Jun 30, 2019

  1. Implement am(u|m).

    pleroy committed Jun 30, 2019
    Copy the full SHA
    963ce51 View commit details
  2. Test for am(u|m).

    pleroy committed Jun 30, 2019
    Copy the full SHA
    302947b View commit details
  3. Correct implementation.

    pleroy committed Jun 30, 2019
    Copy the full SHA
    a1d5c99 View commit details
  4. A benchmark.

    pleroy committed Jun 30, 2019
    Copy the full SHA
    b863859 View commit details
  5. Copy the full SHA
    f81114f View commit details

Commits on Jul 1, 2019

  1. Copy the full SHA
    a8e40f0 View commit details

Commits on Jul 3, 2019

  1. Copy the full SHA
    aab7039 View commit details
  2. Fix tolerance.

    pleroy committed Jul 3, 2019
    Copy the full SHA
    028d450 View commit details
  3. Test around multiples of K.

    pleroy committed Jul 3, 2019
    Copy the full SHA
    4c5c863 View commit details
  4. Lint.

    pleroy committed Jul 3, 2019
    Copy the full SHA
    77e431f View commit details

Commits on Jul 4, 2019

  1. Merge pull request #2233 from pleroy/Amplitude

    Implement Jacobi amplitude
    pleroy authored Jul 4, 2019
    Copy the full SHA
    9efffef View commit details
30 changes: 29 additions & 1 deletion benchmarks/elliptic_functions_benchmark.cpp
Original file line number Diff line number Diff line change
@@ -1,16 +1,43 @@

// .\Release\x64\benchmarks.exe --benchmark_repetitions=10 --benchmark_min_time=2 --benchmark_filter=JacobiSNCNDN // NOLINT(whitespace/line_length)
// .\Release\x64\benchmarks.exe --benchmark_repetitions=10 --benchmark_min_time=2 --benchmark_filter=Jacobi // NOLINT(whitespace/line_length)

#include <random>
#include <vector>

#include "benchmark/benchmark.h"
#include "numerics/elliptic_functions.hpp"
#include "quantities/numbers.hpp"
#include "quantities/quantities.hpp"

namespace principia {
namespace numerics {

using quantities::Angle;

void BM_JacobiAmplitude(benchmark::State& state) {
constexpr int size = 100;

std::mt19937_64 random(42);
std::uniform_real_distribution<> distribution_u(-10.0, 10.0);
std::uniform_real_distribution<> distribution_mc(0.0, 1.0);
std::vector<double> us;
std::vector<double> mcs;
for (int i = 0; i < size; ++i) {
us.push_back(distribution_u(random));
mcs.push_back(distribution_mc(random));
}

while (state.KeepRunningBatch(size * size)) {
Angle a;
for (double const u : us) {
for (double const mc : mcs) {
a = JacobiAmplitude(u, mc);
}
}
benchmark::DoNotOptimize(a);
}
}

void BM_JacobiSNCNDN(benchmark::State& state) {
constexpr int size = 100;

@@ -39,6 +66,7 @@ void BM_JacobiSNCNDN(benchmark::State& state) {
}
}

BENCHMARK(BM_JacobiAmplitude);
BENCHMARK(BM_JacobiSNCNDN);

} // namespace numerics
3 changes: 2 additions & 1 deletion mathematica/elliptic_functions.wl
Original file line number Diff line number Diff line change
@@ -45,7 +45,7 @@
(* ::Input:: *)
(*randomvals=Map[*)
(*{*)
(*JacobiSN[#[[1]],#[[2]]],JacobiCN[#[[1]],#[[2]]],JacobiDN[#[[1]],#[[2]]]}&,*)
(*JacobiSN[#[[1]],#[[2]]],JacobiCN[#[[1]],#[[2]]],JacobiDN[#[[1]],#[[2]]],JacobiAmplitude[#[[1]],#[[2]]]}&,*)
(*randomargs];*)


@@ -56,6 +56,7 @@
(*" value: "<>decimalFloatLiteral[#[[3]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[4]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[5]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[6]],2]<>*)
(*"}"*)
(*&,*)
(*Join[randomargs,randomvals,2]];*)
92 changes: 81 additions & 11 deletions numerics/elliptic_functions.cpp
Original file line number Diff line number Diff line change
@@ -10,15 +10,21 @@
#include "numerics/polynomial_evaluators.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/numbers.hpp"
#include "quantities/si.hpp"

namespace principia {
namespace numerics {
namespace internal_elliptic_functions {

using quantities::Abs;
using quantities::ArcTan;
using quantities::Sqrt;
using quantities::si::Radian;

namespace numerics {
namespace {

constexpr double k_over_2_lower_bound = π / 4.0;

void JacobiSNCNDNReduced(double u, double mc, double& s, double& c, double& d);

// Maclaurin series for Fukushima b₀. These are polynomials in m that are used
@@ -114,8 +120,6 @@ void JacobiSNCNDNReduced(double const u,
d = Sqrt(1.0 - m * yᵢ);
}

} // namespace

// Double precision subroutine to compute three Jacobian elliptic functions
// simultaneously
//
@@ -132,13 +136,12 @@ void JacobiSNCNDNReduced(double const u,
//
// Output: s = sn(u|m), c=cn(u|m), d=dn(u|m)
//
void JacobiSNCNDN(double const u,
double const mc,
double& s,
double& c,
double& d) {
constexpr double k_over_2_lower_bound = π / 4.0;

void JacobiSNCNDNWithK(double const u,
double const mc,
double const k,
double& s,
double& c,
double& 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
@@ -149,7 +152,6 @@ void JacobiSNCNDN(double const u,
if (abs_u < k_over_2_lower_bound) {
JacobiSNCNDNReduced(abs_u, mc, s, c, d);
} else {
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;
@@ -197,6 +199,74 @@ void JacobiSNCNDN(double const u,
s = -s;
}
}
} // namespace

Angle JacobiAmplitude(double u, double mc) {
constexpr double k_over_2_lower_bound = π / 4.0;
double s;
double c;
double d;
double n;
double abs_u = Abs(u);
if (abs_u < k_over_2_lower_bound) {
JacobiSNCNDNReduced(abs_u, mc, s, c, d);
if (u < 0.0) {
s = -s;
}
n = 0.0;
} else {
// We *don't* follow Fukushima, Fast computation of elliptic functions and
// incomplete integrals for constant values of elliptic parameter and
// elliptic characteristic, formula (20). It calls the ArcTan function with
// negative values of c and values of s close to 0, which corresponds to a
// branch cut: ArcTan can jump from -π or +π (or vice-versa) depending on
// the accuracy of s. Similarly the truncation to integer can jump by 1
// depending on the accuracy of k. These problems may result in jumps of
// 2π for the final value of am(u|m).
// Instead, we explicitly reduce u to the range [-k, k] and thus the ArcTan
// to the range [-π/2, π/2]. We avoid the branch cut, and any inaccuracy in
// the rounding has the innocuous effect of causing the ArcTan to go a bit
// beyond -π/2 or π/2.
double const k = EllipticK(mc);
n = std::nearbyint(u / (2.0 * k));
JacobiSNCNDNWithK(u - 2.0 * n * k, mc, k, s, c, d);
}
return n * π * Radian + ArcTan(s, c);
}

// Double precision subroutine to compute three Jacobian elliptic functions
// simultaneously
//
// For general argument: -infty < u < infty
//
// Reference: T. Fukushima, (2012) Numer. Math.
// DOI 10.1007/s00211-012-0498-0
// "Precise and Fast Computation of Jacobian Elliptic Functions by
// Conditional Duplication"
//
// Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
//
// Inputs: u = argument, mc = 1-m, 0 < mc <= 1
//
// Output: s = sn(u|m), c=cn(u|m), d=dn(u|m)
//
void JacobiSNCNDN(double const u,
double const mc,
double& s,
double& c,
double& d) {
double abs_u = Abs(u);
if (abs_u < k_over_2_lower_bound) {
JacobiSNCNDNReduced(abs_u, mc, s, c, d);
if (u < 0.0) {
s = -s;
}
} else {
double const k = EllipticK(mc);
JacobiSNCNDNWithK(u, mc, k, s, c, d);
}
}

} // namespace internal_elliptic_functions
} // namespace numerics
} // namespace principia
12 changes: 12 additions & 0 deletions numerics/elliptic_functions.hpp
Original file line number Diff line number Diff line change
@@ -1,14 +1,26 @@
#pragma once

#include "quantities/quantities.hpp"

// This code is a straightforward translation in C++ of:
// Fukushima, Toshio. (2012). xgscd.txt (Fortran program package to compute the
// Jacobian elliptic functions, sn(u|m), cn(u|m), dn(u|m)).
// Downloaded from:
// https://www.researchgate.net/publication/233903220_xgscdtxt_Fortran_program_package_to_compute_the_Jacobian_elliptic_functions_snum_cnum_dnum
namespace principia {
namespace numerics {
namespace internal_elliptic_functions {

using quantities::Angle;

Angle JacobiAmplitude(double u, double mc);

void JacobiSNCNDN(double u, double mc, double& s, double& c, double& d);

} // namespace internal_elliptic_functions

using internal_elliptic_functions::JacobiAmplitude;
using internal_elliptic_functions::JacobiSNCNDN;

} // namespace numerics
} // namespace principia
Loading