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

Commits on Jul 7, 2019

  1. Copy the full SHA
    246142a View commit details
  2. Unverified

    This commit is not signed, but one or more authors requires that any commit attributed to them is signed.
    Copy the full SHA
    6f2c364 View commit details
  3. Tests passing.

    pleroy committed Jul 7, 2019
    Copy the full SHA
    aeb4ae8 View commit details
  4. Comment.

    pleroy committed Jul 7, 2019
    Copy the full SHA
    e6caf72 View commit details
  5. Faster reduction.

    pleroy committed Jul 7, 2019
    Copy the full SHA
    bd7d4b2 View commit details
  6. Merge pull request #2238 from pleroy/ArgumentReduction

    Simple argument reduction for elliptic integrals
    pleroy authored Jul 7, 2019
    Copy the full SHA
    df57b9d View commit details
91 changes: 89 additions & 2 deletions mathematica/elliptic_integrals.wl
Original file line number Diff line number Diff line change
@@ -85,14 +85,58 @@
(*Join[randomargs2,randomvals2,2]];*)


(* ::Subsubsection:: *)
(*Argument reduction*)


(* ::Input:: *)
(*redm={1/10,1/2,9/10}*)


(* ::Input:: *)
(*red\[CurlyPhi]=Table[i/2,{i,-20,20}]*)


(* ::Input:: *)
(*redargs2=Flatten[Outer[List,red\[CurlyPhi],redm],1];*)


(* ::Input:: *)
(*redvals2=Map[*)
(*N[{*)
(*If[#[[2]]==0,*)
(*Limit[fukushimaB[#[[1]],m],m->0],*)
(*fukushimaB[#[[1]],#[[2]]]],*)
(*If[#[[2]]==0,*)
(*Limit[fukushimaD[#[[1]],m],m->0],*)
(*fukushimaD[#[[1]],#[[2]]]],*)
(*EllipticE[#[[1]],#[[2]]],*)
(*EllipticF[#[[1]],#[[2]]]*)
(*},20]&,*)
(*redargs2];*)


(* ::Input:: *)
(*redstrs2=Map[*)
(*"entry { argument: "<>decimalFloatLiteral[#[[1]],2]<>*)
(*" argument: "<>decimalFloatLiteral[#[[2]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[3]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[4]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[5]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[6]],2]<>*)
(*"}"*)
(*&,*)
(*Join[N[redargs2,20],redvals2,2]];*)


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


(* ::Input:: *)
(*Export[*)
(*"..\\numerics\\bivariate_elliptic_integrals.proto.txt",*)
(*StringRiffle[randomstrs2,"\n"],*)
(*StringRiffle[Flatten[{randomstrs2,redstrs2}],"\n"],*)
(*"text"]*)


@@ -188,14 +232,57 @@
(*Join[random1args3,random1vals3,2]];*)


(* ::Subsubsection:: *)
(*Argument reduction*)


(* ::Input:: *)
(*redm={1/10,1/2,9/10}*)


(* ::Input:: *)
(*red\[CurlyPhi]=Table[i/2,{i,-20,20}]*)


(* ::Input:: *)
(*redn={1/10,1/2,9/10}*)


(* ::Input:: *)
(*redargs3=Flatten[Outer[List,red\[CurlyPhi],redn,redm],2];*)


(* ::Input:: *)
(*redvals3=Map[*)
(*N[{*)
(*If[#[[2]]==0,*)
(*Limit[fukushimaJ[#[[1]],n,#[[3]]],n->0],*)
(*fukushimaJ[#[[1]],#[[2]],#[[3]]]],*)
(*EllipticPi[#[[2]],#[[1]],#[[3]]]*)
(*},20]&,*)
(*redargs3];*)


(* ::Input:: *)
(*redstrs3=Map[*)
(*"entry { argument: "<>decimalFloatLiteral[#[[1]],2]<>*)
(*" argument: "<>decimalFloatLiteral[#[[2]],2]<>*)
(*" argument: "<>decimalFloatLiteral[#[[3]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[4]],2]<>*)
(*" value: "<>decimalFloatLiteral[#[[5]],2]<>*)
(*"}"*)
(*&,*)
(*Join[N[redargs3,20],redvals3,2]];*)


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


(* ::Input:: *)
(*Export[*)
(*"..\\numerics\\trivariate_elliptic_integrals.proto.txt",*)
(*StringRiffle[Flatten[{randomstrs3,random1strs3}],"\n"],*)
(*StringRiffle[Flatten[{randomstrs3,random1strs3,redstrs3}],"\n"],*)
(*"text"]*)


125 changes: 124 additions & 1 deletion numerics/bivariate_elliptic_integrals.proto.txt

Large diffs are not rendered by default.

70 changes: 61 additions & 9 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -91,6 +91,12 @@ double FukushimaEllipticJsMaclaurinSeries(double y, double n, double m);
// Fukushima's T function [Fuku11c].
double FukushimaT(double t, double h);

// Argument reduction: angle = fractional_part + integer_part * π where
// fractional_part is in [-π/2, π/2].
void Reduce(Angle const& angle,
Angle& fractional_part,
std::int64_t& integer_part);

// A generator for the Maclaurin series for q(m) / m where q is Jacobi's nome
// function.
template<int n, template<typename, typename, int> class Evaluator>
@@ -1618,6 +1624,27 @@ double FukushimaT(double const t, double const h) {
}
}

// TODO(phl): This is extremely imprecise near large multiples of π. Use a
// better algorithm (Payne-Hanek?).
void Reduce(Angle const& angle,
Angle& fractional_part,
std::int64_t& integer_part) {
double const angle_in_half_cycles = angle / (π * Radian);
double reduced_in_half_cycles;
#if PRINCIPIA_USE_SSE3_INTRINSICS
auto const& x = angle_in_half_cycles;
__m128d const x_128d = _mm_set_sd(x);
integer_part = _mm_cvtsd_si64(x_128d);
reduced_in_half_cycles = _mm_cvtsd_f64(
_mm_sub_sd(x_128d,
_mm_cvtsi64_sd(__m128d{}, integer_part)));
#else
integer_part = std::nearbyint(angle_in_half_cycles);
reduced_in_half_cycles = angle_in_half_cycles - integer_part;
#endif
fractional_part = reduced_in_half_cycles * π * Radian;
}

} // namespace

// Double precision general incomplete elliptic integrals of all three kinds
@@ -1640,13 +1667,20 @@ void FukushimaEllipticBDJ(Angle const& φ,
double& b,
double& d,
double& j) {
DCHECK_LE(0 * Radian, φ);
DCHECK_GE(π / 2 * Radian, φ);
DCHECK_LE(0, n);
DCHECK_GE(1, n);
DCHECK_LE(0, mc);
DCHECK_GE(1, mc);

// See Appendix B of [Fuku11b] and Appendix A.1 of [Fuku11c] for argument
// reduction.
// TODO(phl): This is extremely imprecise near large multiples of π. Use a
// better algorithm (Payne-Hanek?).
Angle φ_reduced;
std::int64_t count;
Reduce(φ, φ_reduced, count);
Angle const abs_φ_reduced = Abs(φ_reduced);

// 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
@@ -1657,22 +1691,26 @@ void FukushimaEllipticBDJ(Angle const& φ,
constexpr Angle φs = 1.249 * Radian;
constexpr double ys = 0.9;

bool has_computed_complete_integrals = false;
double bc;
double dc;
double jc;

// 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);
if (abs_φ_reduced < φs) {
FukushimaEllipticBsDsJs(Sin(abs_φ_reduced), n, mc, b, d, j);
} else {
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 = Cos(abs_φ_reduced);
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);
@@ -1681,12 +1719,12 @@ void FukushimaEllipticBDJ(Angle const& φ,
b = bc - (b - sz);
d = dc - (d + sz);
j = jc - (j + FukushimaT(t, h));
has_computed_complete_integrals = true;
} else {
double const w²_numerator = mc * (1.0 - c²);
if (w²_numerator < c² * z²_denominator) {
FukushimaEllipticBcDcJc(c, n, mc, b, d, j);
} else {
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);
@@ -1696,9 +1734,25 @@ void FukushimaEllipticBDJ(Angle const& φ,
b = bc - (b - sz);
d = dc - (d + sz);
j = jc - (j + FukushimaT(t, h));
has_computed_complete_integrals = true;
}
}
}

if (φ_reduced < 0.0 * Radian) {
b = -b;
d = -d;
j = -j;
}
if (count != 0) {
double const nc = 1.0 - n;
if (!has_computed_complete_integrals) {
FukushimaEllipticBDJ(nc, mc, bc, dc, jc);
}
b = 2 * count * bc + b;
d = 2 * count * dc + d;
j = 2 * count * jc + j;
}
}

void EllipticEFΠ(quantities::Angle const& φ,
@@ -1707,8 +1761,6 @@ void EllipticEFΠ(quantities::Angle const& φ,
double& e,
double& f,
double& ᴨ) {
DCHECK_LE(0 * Radian, φ);
DCHECK_GE(π / 2 * Radian, φ);
DCHECK_LE(0, n);
DCHECK_GE(1, n);
DCHECK_LE(0, mc);
2 changes: 1 addition & 1 deletion numerics/elliptic_integrals_test.cpp
Original file line number Diff line number Diff line change
@@ -140,7 +140,7 @@ TEST_F(EllipticIntegralsTest, MathematicaBivariate) {
actual_value_f,
actual_value_ᴨ);

EXPECT_THAT(actual_value_b, AlmostEquals(expected_value_b, 0, 5))
EXPECT_THAT(actual_value_b, AlmostEquals(expected_value_b, 0, 6))
<< argument_φ << " " << argument_m;
EXPECT_THAT(actual_value_d, AlmostEquals(expected_value_d, 0, 44))
<< argument_φ << " " << argument_m;
371 changes: 370 additions & 1 deletion numerics/trivariate_elliptic_integrals.proto.txt

Large diffs are not rendered by default.