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

Commits on Jun 16, 2019

  1. Change the identifiers that contain Greek letters to use the proper U…

    …nicode character.
    
    Remove φc as it introduces errors.
    
    Delta.
    pleroy committed Jun 16, 2019
    Copy the full SHA
    ceb5994 View commit details
  2. Merge pull request #2207 from pleroy/Greek

    Use Greek letters instead of pi, phi, etc. and remove phic
    pleroy authored Jun 16, 2019
    Copy the full SHA
    6cde2f4 View commit details
11 changes: 5 additions & 6 deletions benchmarks/elliptic_integrals_benchmark.cpp
Original file line number Diff line number Diff line change
@@ -15,14 +15,14 @@ void BM_FukushimaEllipticBDJ(benchmark::State& state) {
constexpr int size = 20;

std::mt19937_64 random(42);
std::uniform_real_distribution<> distribution_phi(0.0, π / 2);
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> phis;
std::vector<double> φs;
std::vector<double> ns;
std::vector<double> mcs;
for (int i = 0; i < size; ++i) {
phis.push_back(distribution_phi(random));
φs.push_back(distribution_φ(random));
ns.push_back(distribution_n(random));
mcs.push_back(distribution_mc(random));
}
@@ -31,11 +31,10 @@ void BM_FukushimaEllipticBDJ(benchmark::State& state) {
double b;
double d;
double j;
for (double const phi : phis) {
double const phic = π / 2 - phi;
for (double const φ : φs) {
for (double const n : ns) {
for (double const mc : mcs) {
FukushimaEllipticBDJ(phi, phic, n, mc, b, d, j);
FukushimaEllipticBDJ(φ, n, mc, b, d, j);
}
}
}
10 changes: 4 additions & 6 deletions numerics/elliptic_functions_test.cpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@


#include "numerics/elliptic_functions.hpp"

#include "glog/logging.h"
@@ -23,17 +23,15 @@ class EllipticFunctionsTest : public ::testing::Test {};
TEST_F(EllipticFunctionsTest, Xgscd) {
auto const xgscd_expected =
ReadFromTabulatedData(SOLUTION_DIR / "numerics" / "xgscd.proto.txt");
constexpr double PI = 3.1415926535897932384626433;
constexpr double PIHALF = PI * 0.5;
double dmc, mc, m, du, u, s, c, d;
double Δmc, mc, m, du, u, s, c, d;
int jend, iend;
jend = 10;
iend = 8;
dmc = 1.0 / static_cast<double>(jend);
Δmc = 1.0 / static_cast<double>(jend);
std::printf("%10s%10s%25s%25s%25s\n", "m", "u", "s", "c", "d");
int expected_index = 0;
for (int j = 1; j <= jend; ++j) {
mc = static_cast<double>(j) * dmc;
mc = static_cast<double>(j) * Δmc;
m = 1.0 - mc;
du = EllipticK(mc) / static_cast<double>(iend);
for (int i = 0; i <= iend * 8; ++i) {
48 changes: 17 additions & 31 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -161,11 +161,6 @@ double BulirschCel(double const kc0,
void FukushimaEllipticBD(double const mc, double& elb, double& eld) {
double elk, m, mx, kkc, nome;

constexpr double PIQ = 0.78539816339744830961566084581988;
constexpr double PIHALF = 1.5707963267948966192313216916398;
constexpr double PI = 3.1415926535897932384626433832795;
constexpr double PIINV = 0.31830988618379067153776752674503;

constexpr double Q1 = 1.0 / 16.0;
constexpr double Q2 = 1.0 / 32.0;
constexpr double Q3 = 21.0 / 1024.0;
@@ -213,8 +208,8 @@ void FukushimaEllipticBD(double const mc, double& elb, double& eld) {

m = 1.0 - mc;
if (m < 1.11e-16) {
elb = PIQ;
eld = PIQ;
elb = π / 4;
eld = π / 4;
} else if (mc < 1.11e-16) {
elb = 1.0;
eld = 0.3862943611198906188344642429164 - 0.5 * log(mc);
@@ -230,7 +225,7 @@ void FukushimaEllipticBD(double const mc, double& elb, double& eld) {
} else {
mx = mc - 0.05;

// (K'-1)/(pi/2)
// (K'-1)/(π / 2)
dkkc = 0.01286425658832983978282698630501405107893
+ mx * (0.26483429894479586582278131697637750604652
+ mx * (0.15647573786069663900214275050014481397750
@@ -245,7 +240,7 @@ void FukushimaEllipticBD(double const mc, double& elb, double& eld) {
+ mx * (0.04978344771840508342564702588639140680363
+ mx * (0.04812375496807025605361215168677991360500))))))))))));

// (K'-E')/(pi/2)
// (K'-E')/(π / 2)
dddc = 0.02548395442966088473597712420249483947953
+ mx * (0.51967384324140471318255255900132590084179
+ mx * (0.20644951110163173131719312525729037023377
@@ -269,10 +264,10 @@ void FukushimaEllipticBD(double const mc, double& elb, double& eld) {
elb = 1.0 + delb;
eld = elk1 - delb;
} else if (m <= 0.01) {
elb = PIHALF*(B1 + m * (B2 + m * (B3 + m * (B4 + m * (B5 + m * (B6 + m *
(B7 + m * B8)))))));
eld = PIHALF*(D1 + m * (D2 + m * (D3 + m * (D4 + m * (D5 + m * (D6 + m *
(D7 + m * D8)))))));
elb = (π / 2) * (B1 + m * (B2 + m * (B3 + m * (B4 + m * (B5 + m * (B6 + m *
(B7 + m * B8)))))));
eld = (π / 2) * (D1 + m * (D2 + m * (D3 + m * (D4 + m * (D5 + m * (D6 + m *
(D7 + m * D8)))))));
} else if (m <= 0.1) {
mx = 0.95 - mc;
elb = 0.790401413584395132310045630540381158921005
@@ -1350,25 +1345,18 @@ double FukushimaT(double const t, double const h) {
//
// Author: T. Fukushima Toshio.Fukushima@nao.ac.jp
//
// Used subprograms: cel,celbd,celbdj,elcbdj,serbd,serj,uatan
//
// Inputs: phi = argument 0 <= phi <= PI/2
// phic = complementar argument 0 <= phic <= PI/2
// n = characteristic 0 <= n <= 1
// mc = complementary parameter 0 <= mc <= 1
// Inputs: φ = argument 0 <= φ <= π / 2
// n = characteristic 0 <= n <= 1
// mc = complementary parameter 0 <= mc <= 1
//
// Outputs: b, d, j
//
// CAUTION: phi and phic must satisfy condition, phi + phic = PI/2
//
void FukushimaEllipticBDJ(double const phi,
double const phic,
void FukushimaEllipticBDJ(double const φ,
double const n,
double const mc,
double& b,
double& d,
double& j) {
// TODO(phl): CHECK_EQ(π / 2, phi + phic);
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
@@ -1377,13 +1365,13 @@ void FukushimaEllipticBDJ(double const phi,
// 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
// sure if it has an impact on correctness.
if (phi < 1.249) {
FukushimaEllipticBsDsJs(sin(phi), n, mc, b, d, j);
if (φ < 1.249) {
FukushimaEllipticBsDsJs(sin(φ), n, mc, b, d, j);
} else {
m = 1.0 - mc;
nc = 1.0 - n;
h = n * nc * (n - m);
c = sin(phic);
c = cos);
x = c * c;
d2 = mc + m * x;
if (x < 0.9 * d2) {
@@ -1441,13 +1429,11 @@ double EllipticK(double const mc) {
constexpr double D12 = 40784671953.0 / 8796093022208.0;
constexpr double D13 = 9569130097211.0 / 2251799813685248.0;
constexpr double D14 = 17652604545791.0 / 4503599627370496.0;
constexpr double PIHALF = π / 2;
constexpr double PIINV = 0.5 / PIHALF;
constexpr double TINY = 1.0e-99;

m = 1.0 - mc;
if (abs(m) < 1.0e-16) {
return PIHALF;
return π / 2;
} else if (mc < TINY) {
return 1.3862943611198906 - 0.5 * log(TINY);
} else if (mc < 1.11e-16) {
@@ -1471,7 +1457,7 @@ double EllipticK(double const mc) {
0.091439629201749751 + mx * (
0.085842591595413900 + mx * (
0.081541118718303215))))))))));
return -kkc * PIINV * log(nome);
return -kkc * (1 / π) * log(nome);
} else if (m <= 0.1) {
mx = m - 0.05;
return 1.591003453790792180 + mx * (
5 changes: 2 additions & 3 deletions numerics/elliptic_integrals.hpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#pragma once
#pragma once

// This code is a straightforward translation in C++ of:
// Fukushima, Toshio. (2018). xelbdj.txt: Fortran test driver for
@@ -10,8 +10,7 @@
namespace principia {
namespace numerics {

void FukushimaEllipticBDJ(double phi,
double phic,
void FukushimaEllipticBDJ(double φ,
double n,
double mc,
double& b,
2 changes: 2 additions & 0 deletions numerics/elliptic_integrals.proto.txt
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
# n m φ / π B D J
#--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
entry { argument: 2.1010855918191289907e-01 argument: 9.9903585338494881831e-01 argument: 3.5068934540537104588e-01 value: 3.4353825891999272e-01 value: 1.456087381239668e-02 value: 1.478463239022e-02}
entry { argument: 2.1010855918191289907e-01 argument: 9.9903585338494881831e-01 argument: 4.0310686403184424756e-01 value: 3.9226736277424449e-01 value: 2.221037436174530e-02 value: 2.265998746066e-02}
entry { argument: 2.1010855918191289907e-01 argument: 9.9903585338494881831e-01 argument: 6.5294929932643112307e-01 value: 6.0748468394713446e-01 value: 9.74527946532982e-02 value: 1.0251627307820e-01}
42 changes: 24 additions & 18 deletions numerics/elliptic_integrals_test.cpp
Original file line number Diff line number Diff line change
@@ -10,6 +10,8 @@
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/serialization.hpp"

#include "quantities/quantities.hpp"

namespace principia {

using testing_utilities::AlmostEquals;
@@ -26,22 +28,20 @@ class EllipticIntegralsTest : public ::testing::Test {};
TEST_F(EllipticIntegralsTest, Xelbdj) {
auto const xeldbj_expected =
ReadFromTabulatedData(SOLUTION_DIR / "numerics" / "xelbdj.proto.txt");
constexpr double PI = 3.1415926535897932384626433;
constexpr double PIHALF = 1.5707963267948966192313216916398;
double dnc, dmc, dphi, nc, nn, mc, mm, phi, phic, b, d, j;
double Δnc, Δmc, Δφ, nc, nn, mc, mm, φ, b, d, j;
int lend, kend, iend;

lend = 5;
kend = 5;
iend = 4;
dnc = 1.0 / static_cast<double>(lend - 1);
dmc = 1.0 / static_cast<double>(kend - 1);
dphi = PIHALF / static_cast<double>(iend);
Δnc = 1.0 / static_cast<double>(lend - 1);
Δmc = 1.0 / static_cast<double>(kend - 1);
Δφ = (π / 2) / static_cast<double>(iend);
std::printf(
"%10s%10s%10s%25s%25s%25s\n", "n", "m", "phi/PI", "elb", "eld", "elj");
"%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) * dnc;
nc = static_cast<double>(l - 1) * Δnc;
if (nc <= 1.05e-8) {
nc = 1.05e-8;
}
@@ -52,36 +52,43 @@ TEST_F(EllipticIntegralsTest, Xelbdj) {
nn = 1.0 - nc;
for (int k = 1; k <= kend; ++k) {
std::printf("\n");
mc = static_cast<double>(k - 1) * dmc;
mc = static_cast<double>(k - 1) * Δmc;
if (mc <= 0.0) {
mc = 1.21e-32;
}
mm = 1.0 - mc;
for (int i = 0; i <= iend; ++i) {
phi = dphi * static_cast<double>(i);
phic = dphi * static_cast<double>(iend - i);
FukushimaEllipticBDJ(phi, phic, nn, mc, b, d, j);
φ = Δφ * 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,
phi / PI,
φ / π,
b,
d,
j);

auto const& expected_entry = xeldbj_expected.entry(expected_index);
auto const expected_argument_n = expected_entry.argument(0);
auto const expected_argument_m = expected_entry.argument(1);
auto const expected_argument_phi_over_pi = expected_entry.argument(2);
auto const expected_argument_φ_over_π = expected_entry.argument(2);
auto const expected_value_b = expected_entry.value(0);
auto const expected_value_d = expected_entry.value(1);
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(phi / PI, IsNear(expected_argument_phi_over_pi, 1.001));
EXPECT_THAT(φ / π, IsNear(expected_argument_φ_over_π, 1.001));
EXPECT_THAT(b, AlmostEquals(expected_value_b, 0, 8));
EXPECT_THAT(d, AlmostEquals(expected_value_d, 0, 97));
EXPECT_THAT(j, AlmostEquals(expected_value_j, 0, 135));

// TODO(phl): xelbdj_all.txt enshrines values that are incorrect for
// m close to 1 and φ close to (the machine representation of) π / 2.
// That seems to stem from the computation of the complementary angle
// φc = π / 2 - φ, which introduces a 1 ULP inconsistency between φc and
// φ.
if (mm != 1.0 || i != iend) {
EXPECT_THAT(d, AlmostEquals(expected_value_d, 0, 97));
EXPECT_THAT(j, AlmostEquals(expected_value_j, 0, 135));
}
++expected_index;
}
}
@@ -106,7 +113,6 @@ TEST_F(EllipticIntegralsTest, MathematicaMNear1) {
double actual_value_d;
double actual_value_j;
FukushimaEllipticBDJ(argument_φ,
π / 2 - argument_φ,
argument_n,
1.0 - argument_m,
actual_value_b,
1 change: 1 addition & 0 deletions numerics/numerics.vcxproj
Original file line number Diff line number Diff line change
@@ -65,6 +65,7 @@
<ClCompile Include="чебышёв_series_test.cpp" />
</ItemGroup>
<ItemGroup>
<Text Include="elliptic_integrals.proto.txt" />
<Text Include="xelbdj.proto.txt" />
<Text Include="xgscd.proto.txt" />
</ItemGroup>
3 changes: 3 additions & 0 deletions numerics/numerics.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -192,5 +192,8 @@
<Text Include="xelbdj.proto.txt">
<Filter>Resource Files</Filter>
</Text>
<Text Include="elliptic_integrals.proto.txt">
<Filter>Resource Files</Filter>
</Text>
</ItemGroup>
</Project>
2 changes: 1 addition & 1 deletion numerics/xelbdj.proto.txt
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# These tabulated values come from the output given by Fukushima at the end of
# his xelbdj_all.txt.
# n m phi/PI B D J
# n m φ / π B D J
#----------------------------------------------------------------------------------------------------------------------------------------------------------------
entry { argument: 0.99976, argument: 1.00000, argument: 0.00000, value: 0.000000000000000, value: 0.000000000000000, value: 0.000000000000000 }
entry { argument: 0.99976, argument: 1.00000, argument: 0.12500, value: 0.382683432365090, value: 0.020516286796422, value: 0.022570343177599 }