Skip to content

Commit

Permalink
Add typing and better identifiers to FukushimaEllipticBDJ.
Browse files Browse the repository at this point in the history
pleroy committed Jun 16, 2019
1 parent ceb5994 commit f9b07aa
Showing 4 changed files with 78 additions and 47 deletions.
11 changes: 8 additions & 3 deletions benchmarks/elliptic_integrals_benchmark.cpp
Original file line number Diff line number Diff line change
@@ -7,8 +7,13 @@
#include "benchmark/benchmark.h"
#include "numerics/elliptic_integrals.hpp"
#include "quantities/numbers.hpp"
#include "quantities/si.hpp"

namespace principia {

using quantities::Angle;
using quantities::si::Radian;

namespace numerics {

void BM_FukushimaEllipticBDJ(benchmark::State& state) {
@@ -18,11 +23,11 @@ void BM_FukushimaEllipticBDJ(benchmark::State& state) {
std::uniform_real_distribution<> distribution_φ(0.0, π / 2);
std::uniform_real_distribution<> distribution_n(0.0, 1.0);
std::uniform_real_distribution<> distribution_mc(0.0, 1.0);
std::vector<double> φs;
std::vector<Angle> φs;
std::vector<double> ns;
std::vector<double> mcs;
for (int i = 0; i < size; ++i) {
φs.push_back(distribution_φ(random));
φs.push_back(distribution_φ(random) * Radian);
ns.push_back(distribution_n(random));
mcs.push_back(distribution_mc(random));
}
@@ -31,7 +36,7 @@ void BM_FukushimaEllipticBDJ(benchmark::State& state) {
double b;
double d;
double j;
for (double const φ : φs) {
for (Angle const φ : φs) {
for (double const n : ns) {
for (double const mc : mcs) {
FukushimaEllipticBDJ(φ, n, mc, b, d, j);
72 changes: 46 additions & 26 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -2,7 +2,10 @@
#include "glog/logging.h"

#include "numerics/elliptic_integrals.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/numbers.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"

// TODO(phl):
// 1. Use arrays for the coefficients.
@@ -24,6 +27,13 @@
// Mathematical Functions.

namespace principia {

using quantities::Angle;
using quantities::Cos;
using quantities::Sin;
using quantities::Sqrt;
using quantities::si::Radian;

namespace numerics {

namespace {
@@ -1351,48 +1361,58 @@ double FukushimaT(double const t, double const h) {
//
// Outputs: b, d, j
//
void FukushimaEllipticBDJ(double const φ,
void FukushimaEllipticBDJ(Angle const& φ,
double const n,
double const mc,
double& b,
double& d,
double& j) {
double m, nc, h, c, x, d2, z, bc, dc, jc, sz, t, v, t2;

// NOTE(phl): The original Fortran code had 1.345, which, according to the
// above-mentioned paper, is suitable for single precision. However, this is
// double precision. Importantly, this number should be roughly
// ArcSin[Sqrt[0.9]] where 0.9 is the factor appearing in the next if
// statement. The discrepancy has a 5-10% impact on performance. I am not
// NOTE(phl): The original Fortran code had φs = 1.345 * Radian, which,
// according to the above-mentioned paper, is suitable for single precision.
// However, this is double precision. Importantly, this doesn't match the
// value of ys. The discrepancy has a 5-10% impact on performance. I am not
// sure if it has an impact on correctness.
if (φ < 1.249) {
FukushimaEllipticBsDsJs(sin(φ), n, mc, b, d, j);

// Sin(φs)^2 must be approximately ys.
constexpr Angle φs = 1.249 * Radian;
constexpr double ys = 0.9;

// The selection rule in [Fuku11b] section 2.1, equations (7-11) and [Fuku11c]
// section 3.2, equations (22) and (23). The identifiers follow Fukushima's
// notation.
// NOTE(phl): The computation of 1 - c² loses accuracy with respect to the
// evaluation of Sin(φ).
if (φ < φs) {
FukushimaEllipticBsDsJs(Sin(φ), n, mc, b, d, j);
} else {
m = 1.0 - mc;
nc = 1.0 - n;
h = n * nc * (n - m);
c = cos(φ);
x = c * c;
d2 = mc + m * x;
if (x < 0.9 * d2) {
z = c / sqrt(d2);
double const m = 1.0 - mc;
double const nc = 1.0 - n;
double const h = n * nc * (n - m);
double const c = Cos(φ);
double const c² = c * c;
double const z²_denominator = mc + m * c²;
if (c² < ys * z²_denominator) {
double bc, dc, jc;
double const z = c / Sqrt(z²_denominator);
FukushimaEllipticBsDsJs(z, n, mc, b, d, j);
FukushimaEllipticBDJ(nc, mc, bc, dc, jc);
sz = z * sqrt(1.0 - x);
t = sz / nc;
double const sz = z * Sqrt(1.0 - );
double const t = sz / nc;
b = bc - (b - sz);
d = dc - (d + sz);
j = jc - (j + FukushimaT(t, h));
} else {
v = mc * (1.0 - x);
if (v < x * d2) {
double const w²_numerator = mc * (1.0 - );
if (w²_numerator < * z²_denominator) {
FukushimaEllipticBcDcJc(c, n, mc, b, d, j);
} else {
t2 = (1.0 - x) / d2;
FukushimaEllipticBcDcJc(sqrt(mc * t2), n, mc, b, d, j);
double bc, dc, jc;
double const w²_denominator = z²_denominator;
double const w²_over_mc = (1.0 - c²) / w²_denominator;
FukushimaEllipticBcDcJc(Sqrt(mc * w²_over_mc), n, mc, b, d, j);
FukushimaEllipticBDJ(nc, mc, bc, dc, jc);
sz = c * sqrt(t2);
t = sz / nc;
double const sz = c * Sqrt(w²_over_mc);
double const t = sz / nc;
b = bc - (b - sz);
d = dc - (d + sz);
j = jc - (j + FukushimaT(t, h));
4 changes: 3 additions & 1 deletion numerics/elliptic_integrals.hpp
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#pragma once

#include "quantities/quantities.hpp"

// This code is a straightforward translation in C++ of:
// Fukushima, Toshio. (2018). xelbdj.txt: Fortran test driver for
// "elbdj"/"relbdj", subroutines to compute the double/single precision general
@@ -10,7 +12,7 @@
namespace principia {
namespace numerics {

void FukushimaEllipticBDJ(double φ,
void FukushimaEllipticBDJ(quantities::Angle const& φ,
double n,
double mc,
double& b,
38 changes: 21 additions & 17 deletions numerics/elliptic_integrals_test.cpp
Original file line number Diff line number Diff line change
@@ -6,6 +6,7 @@
#include "numerics/elliptic_integrals.hpp"
#include "gtest/gtest.h"
#include "quantities/numbers.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/serialization.hpp"
@@ -14,6 +15,8 @@

namespace principia {

using quantities::Angle;
using quantities::si::Radian;
using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::ReadFromTabulatedData;
@@ -28,42 +31,42 @@ class EllipticIntegralsTest : public ::testing::Test {};
TEST_F(EllipticIntegralsTest, Xelbdj) {
auto const xeldbj_expected =
ReadFromTabulatedData(SOLUTION_DIR / "numerics" / "xelbdj.proto.txt");
double Δnc, Δmc, Δφ, nc, nn, mc, mm, φ, b, d, j;
int lend, kend, iend;

lend = 5;
kend = 5;
iend = 4;
Δnc = 1.0 / static_cast<double>(lend - 1);
Δmc = 1.0 / static_cast<double>(kend - 1);
Δφ = (π / 2) / static_cast<double>(iend);

constexpr int lend = 5;
constexpr int kend = 5;
constexpr int iend = 4;
constexpr double Δnc = 1.0 / static_cast<double>(lend - 1);
constexpr double Δmc = 1.0 / static_cast<double>(kend - 1);
constexpr Angle Δφ = (π / 2) * Radian / static_cast<double>(iend);

This comment has been minimized.

Copy link
@eggrobin

eggrobin Jun 16, 2019

Member

I don't think the explicit casts are needed.

This comment has been minimized.

Copy link
@pleroy

pleroy Jun 16, 2019

Author Member

Good point. They were in the Fortran code, of course.

(No LGTM because of this?)

std::printf(
"%10s%10s%10s%25s%25s%25s\n", "n", "m", "φ / π", "elb", "eld", "elj");

int expected_index = 0;
for (int l = 1; l <= lend; ++l) {
nc = static_cast<double>(l - 1) * Δnc;
double nc = static_cast<double>(l - 1) * Δnc;
if (nc <= 1.05e-8) {
nc = 1.05e-8;
}
float rnc = static_cast<float>(nc);
if (rnc <= 2.44e-4) {
nc = 2.44e-4;
}
nn = 1.0 - nc;
double const nn = 1.0 - nc;
for (int k = 1; k <= kend; ++k) {
std::printf("\n");
mc = static_cast<double>(k - 1) * Δmc;
double mc = static_cast<double>(k - 1) * Δmc;
if (mc <= 0.0) {
mc = 1.21e-32;
}
mm = 1.0 - mc;
double const mm = 1.0 - mc;
for (int i = 0; i <= iend; ++i) {
φ = Δφ * static_cast<double>(i);
double b, d, j;
Angle const φ = Δφ * static_cast<double>(i);
FukushimaEllipticBDJ(φ, nn, mc, b, d, j);
std::printf("%10.5f%10.5f%10.5f%25.15f%25.15f%25.15f\n",
nn,
mm,
φ / π,
φ / (π * Radian),
b,
d,
j);
@@ -77,7 +80,8 @@ TEST_F(EllipticIntegralsTest, Xelbdj) {
auto const expected_value_j = expected_entry.value(2);
EXPECT_THAT(nn, IsNear(expected_argument_n, 1.001));
EXPECT_THAT(mm, IsNear(expected_argument_m, 1.001));
EXPECT_THAT(φ / π, IsNear(expected_argument_φ_over_π, 1.001));
EXPECT_THAT(φ / (π * Radian),
IsNear(expected_argument_φ_over_π, 1.001));
EXPECT_THAT(b, AlmostEquals(expected_value_b, 0, 8));

// TODO(phl): xelbdj_all.txt enshrines values that are incorrect for
@@ -104,7 +108,7 @@ TEST_F(EllipticIntegralsTest, MathematicaMNear1) {
for (auto const& entry : elliptic_integrals_expected.entry()) {
double const argument_n = entry.argument(0);
double const argument_m = entry.argument(1);
double const argument_φ = entry.argument(2);
Angle const argument_φ = entry.argument(2) * Radian;
double const expected_value_b = entry.value(0);
double const expected_value_d = entry.value(1);
double const expected_value_j = entry.value(2);

0 comments on commit f9b07aa

Please sign in to comment.