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

Commits on Jun 14, 2019

  1. Copy the full SHA
    f2781e6 View commit details

Commits on Jun 15, 2019

  1. Unverified

    This commit is not signed, but one or more authors requires that any commit attributed to them is signed.
    Copy the full SHA
    a09a66b View commit details
  2. Cleanup.

    pleroy committed Jun 15, 2019
    Copy the full SHA
    9ecc65f View commit details
  3. Copy the full SHA
    8a1cb3d View commit details
  4. After egg's review.

    pleroy committed Jun 15, 2019
    Copy the full SHA
    3b046eb View commit details
  5. Merge pull request #2202 from pleroy/BetterTests

    Move the elliptic integrals test to a proto, add a new test, fix bugs
    pleroy authored Jun 15, 2019
    Copy the full SHA
    ae45711 View commit details
83 changes: 83 additions & 0 deletions mathematica/elliptic_integrals.wl
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
(* ::Package:: *)

(* ::Input:: *)
(*fukushimaB[\[CurlyPhi]_,m_]:=(EllipticE[\[CurlyPhi],m]-(1-m)EllipticF[\[CurlyPhi],m])/m*)


(* ::Input:: *)
(*fukushimaD[\[CurlyPhi]_,m_]:=(EllipticF[\[CurlyPhi],m]-EllipticE[\[CurlyPhi],m])/m*)


(* ::Input:: *)
(*fukushimaJ[\[CurlyPhi]_,n_,m_]:=(EllipticPi[n,\[CurlyPhi],m]-EllipticF[\[CurlyPhi],m])/n*)


(* ::Input:: *)
(*decimalFloatLiteral[x_Real,exponentWidth_Integer]:=*)
(* With[*)
(* {m=MantissaExponent[x][[1]]*10,*)
(* e=MantissaExponent[x][[2]]-1},*)
(* StringJoin[*)
(* StringRiffle[#,""]&/@*)
(* {{#[[1]]},*)
(* If[Length[#]>1,{"."},Nothing],*)
(* If[Length[#]>1,StringPartition[#[[2]],UpTo[5]],Nothing],*)
(* If[e!=0,"e"<>If[e>0,"+","-"]<>IntegerString[e,10,exponentWidth],Nothing]}&[*)
(* StringSplit[ToString[m],"."]]]]*)


(* ::Input:: *)
(*SeedRandom[666]*)


(* ::Input:: *)
(*random\[CurlyPhi]=Sort[RandomReal[{0,\[Pi]/2},10,WorkingPrecision->20]]*)


(* ::Input:: *)
(*randomn=Sort[RandomReal[{0,1},10,WorkingPrecision->20]]*)


(* ::Input:: *)
(*randomm=Sort[RandomReal[{999/1000,1},10,WorkingPrecision->20]]*)


(* ::Input:: *)
(*args=Flatten[Outer[List,randomn,randomm,random\[CurlyPhi]],2];*)


(* ::Input:: *)
(*vals=Map[*)
(*{*)
(*#[[1]],#[[2]],#[[3]],*)
(*fukushimaB[#[[3]],#[[2]]],*)
(*fukushimaD[#[[3]],#[[2]]],*)
(*fukushimaJ[#[[3]],#[[1]],#[[2]]]*)
(*}&,*)
(*args,*)
(*{1}];*)


(* ::Input:: *)
(*strs=Map[*)
(*"entry { argument: "<>decimalFloatLiteral[#[[1]],2]<>*)
(*" argument: "<>decimalFloatLiteral[#[[2]],2]<>*)
(*" argument: "<>decimalFloatLiteral[#[[3]],2]<>*)
(*" value: "<>decimalFloatLiteral[fukushimaB[#[[3]],#[[2]]],2]<>*)
(*" value: "<>decimalFloatLiteral[fukushimaD[#[[3]],#[[2]]],2]<>*)
(*" value: "<>decimalFloatLiteral[fukushimaJ[#[[3]],#[[1]],#[[2]]],2] <>*)
(*"}"*)
(*&,*)
(*args,*)
(*{1}];*)


(* ::Input:: *)
(*SetDirectory[NotebookDirectory[]]*)


(* ::Input:: *)
(*Export[*)
(*"..\\numerics\\elliptic_integrals.proto.txt",*)
(*StringRiffle[strs,"\n"],*)
(*"text"]*)
37 changes: 10 additions & 27 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -15,39 +15,21 @@
namespace principia {
namespace numerics {

double Cel(double const kc0,
double const nc,
double const aa,
double const bb,
int& err);
double Cel(double kc0, double nc, double aa, double bb, int& err);

void Celbd(double const mc, double& elb, double& eld);
void Celbd(double mc, double& elb, double& eld);

void Celbdj(double const nc,
double const mc,
double& bc,
double& dc,
double& jc);
void Celbdj(double nc, double mc, double& bc, double& dc, double& jc);

void Elcbdj(double const c0,
double const n,
double const mc,
double& b,
double& d,
double& j);
void Elcbdj(double c0, double n, double mc, double& b, double& d, double& j);

void Elsbdj(double const s0,
double const n,
double const mc,
double& b,
double& d,
double& j);
void Elsbdj(double s0, double n, double mc, double& b, double& d, double& j);

void Serbd(double const y, double const m, double& b, double& d);
void Serbd(double y, double m, double& b, double& d);

double Serj(double const y, double const n, double const m);
double Serj(double y, double n, double m);

double Uatan(double const t, double const h);
double Uatan(double t, double h);

// Double precision general incomplete elliptic integrals of all three kinds
//
@@ -697,7 +679,8 @@ void Elcbdj(double const c0,
return;
}
m = 1.0 - mc;
// TODO(phl): Speculative: h = n * (1.0 - n) * (n - m);
h = n * (1.0 - n) * (n - m);
yy[1] = y;
ss[1] = s;
int i;
for (i = 1; i <= 10; ++i) {
8 changes: 4 additions & 4 deletions numerics/elliptic_integrals.hpp
Original file line number Diff line number Diff line change
@@ -10,10 +10,10 @@
namespace principia {
namespace numerics {

void Elbdj(double const phi,
double const phic,
double const n,
double const mc,
void Elbdj(double phi,
double phic,
double n,
double mc,
double& b,
double& d,
double& j);
1,000 changes: 1,000 additions & 0 deletions numerics/elliptic_integrals.proto.txt

Large diffs are not rendered by default.

372 changes: 0 additions & 372 deletions numerics/elliptic_integrals_test.cc

This file was deleted.

125 changes: 125 additions & 0 deletions numerics/elliptic_integrals_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@

#include <filesystem>
#include <vector>
#include <utility>

#include "numerics/elliptic_integrals.hpp"
#include "gtest/gtest.h"
#include "quantities/numbers.hpp"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/serialization.hpp"

namespace principia {

using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::ReadFromTabulatedData;

namespace numerics {

class EllipticIntegralsTest : public ::testing::Test {};

// TODO(phl): This is far from covering all the code paths. In particular, it
// doesn't seem to go through Elcbdj.
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;
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);
std::printf(
"%10s%10s%10s%25s%25s%25s\n", "n", "m", "phi/PI", "elb", "eld", "elj");
int expected_index = 0;
for (int l = 1; l <= lend; ++l) {
nc = static_cast<double>(l - 1) * dnc;
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;
for (int k = 1; k <= kend; ++k) {
std::printf("\n");
mc = static_cast<double>(k - 1) * dmc;
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);
Elbdj(phi, phic, 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_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(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));
++expected_index;
}
}
}
}

// Tabulated values produced with Mathematica for m close to 1 (where there were
// bugs).
TEST_F(EllipticIntegralsTest, MathematicaMNear1) {
auto const elliptic_integrals_expected = ReadFromTabulatedData(
SOLUTION_DIR / "numerics" / "elliptic_integrals.proto.txt");

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);
double const expected_value_b = entry.value(0);
double const expected_value_d = entry.value(1);
double const expected_value_j = entry.value(2);

double actual_value_b;
double actual_value_d;
double actual_value_j;
Elbdj(argument_φ,
π / 2 - argument_φ,
argument_n,
1.0 - argument_m,
actual_value_b,
actual_value_d,
actual_value_j);

EXPECT_THAT(actual_value_b, AlmostEquals(expected_value_b, 0, 7))
<< argument_n << " " << argument_m << " " << argument_φ;
EXPECT_THAT(actual_value_d, AlmostEquals(expected_value_d, 0, 27))
<< argument_n << " " << argument_m << " " << argument_φ;
EXPECT_THAT(actual_value_j, AlmostEquals(expected_value_j, 0, 18835))
<< argument_n << " " << argument_m << " " << argument_φ;
}
}

} // namespace numerics
} // namespace principia
3 changes: 2 additions & 1 deletion numerics/numerics.vcxproj
Original file line number Diff line number Diff line change
@@ -47,7 +47,7 @@
<ClCompile Include="combinatorics_test.cpp" />
<ClCompile Include="double_precision_test.cpp" />
<ClCompile Include="elliptic_integrals.cpp" />
<ClCompile Include="elliptic_integrals_test.cc" />
<ClCompile Include="elliptic_integrals_test.cpp" />
<ClCompile Include="elliptic_functions.cpp" />
<ClCompile Include="elliptic_functions_test.cpp" />
<ClCompile Include="fast_sin_cos_2π.cpp" />
@@ -65,6 +65,7 @@
<ClCompile Include="чебышёв_series_test.cpp" />
</ItemGroup>
<ItemGroup>
<Text Include="xelbdj.proto.txt" />
<Text Include="xgscd.proto.txt" />
</ItemGroup>
</Project>
9 changes: 6 additions & 3 deletions numerics/numerics.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -175,19 +175,22 @@
<ClCompile Include="elliptic_integrals.cpp">
<Filter>Source Files</Filter>
</ClCompile>
<ClCompile Include="elliptic_integrals_test.cc">
<Filter>Test Files</Filter>
</ClCompile>
<ClCompile Include="elliptic_functions.cpp">
<Filter>Source Files</Filter>
</ClCompile>
<ClCompile Include="elliptic_functions_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
<ClCompile Include="elliptic_integrals_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<Text Include="xgscd.proto.txt">
<Filter>Resource Files</Filter>
</Text>
<Text Include="xelbdj.proto.txt">
<Filter>Resource Files</Filter>
</Text>
</ItemGroup>
</Project>
153 changes: 153 additions & 0 deletions numerics/xelbdj.proto.txt

Large diffs are not rendered by default.