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

Commits on Jun 20, 2019

  1. Rewrite EllipticK.

    pleroy committed Jun 20, 2019
    Copy the full SHA
    05f7611 View commit details

Commits on Jun 21, 2019

  1. Merge pull request #2214 from pleroy/Elk

    Rewrite EllipticK
    pleroy authored Jun 21, 2019
    Copy the full SHA
    292a2a2 View commit details
Showing with 209 additions and 189 deletions.
  1. +207 −189 numerics/elliptic_integrals.cpp
  2. +2 −0 quantities/numbers.hpp
396 changes: 207 additions & 189 deletions numerics/elliptic_integrals.cpp
Original file line number Diff line number Diff line change
@@ -160,6 +160,194 @@ using FukushimaTMaclaurin10 = FukushimaTMaclaurin<10, EstrinEvaluator>;
using FukushimaTMaclaurin11 = FukushimaTMaclaurin<11, EstrinEvaluator>;
using FukushimaTMaclaurin12 = FukushimaTMaclaurin<12, EstrinEvaluator>;

// Polynomials for EllipticK. The last part of the name indicates the value of
// m around which the approximation is valid.
// TODO(phl): Some of these polynomials use the Horner scheme because the Estrin
// scheme introduces inaccuracies. Investigate why.
PolynomialInMonomialBasis<double, double, 10, HornerEvaluator> const
k_maclaurin_0_05(std::make_tuple(1.591003453790792180,
0.416000743991786912,
0.245791514264103415,
0.179481482914906162,
0.144556057087555150,
0.123200993312427711,
0.108938811574293531,
0.098853409871592910,
0.091439629201749751,
0.085842591595413900,
0.081541118718303215));
PolynomialInMonomialBasis<double, double, 12, EstrinEvaluator> const
k_maclaurin_0_15(std::make_tuple(1.635256732264579992,
0.471190626148732291,
0.309728410831499587,
0.252208311773135699,
0.226725623219684650,

0.215774446729585976,
0.213108771877348910,
0.216029124605188282,
0.223255831633057896,
0.234180501294209925,

0.248557682972264071,
0.266363809892617521,
0.287728452156114668));
PolynomialInMonomialBasis<double, double, 11, EstrinEvaluator> const
k_maclaurin_0_25(std::make_tuple(1.685750354812596043,
0.541731848613280329,
0.401524438390690257,
0.369642473420889090,
0.376060715354583645,

0.405235887085125919,
0.453294381753999079,
0.520518947651184205,
0.609426039204995055,
0.724263522282908870,

0.871013847709812357,
1.057652872753547036));
PolynomialInMonomialBasis<double, double, 12, EstrinEvaluator> const
k_maclaurin_0_35(std::make_tuple(1.744350597225613243,
0.634864275371935304,
0.539842564164445538,
0.571892705193787391,
0.670295136265406100,

0.832586590010977199,
1.073857448247933265,
1.422091460675497751,
1.920387183402304829,
2.632552548331654201,

3.652109747319039160,
5.115867135558865806,
7.224080007363877411));
PolynomialInMonomialBasis<double, double, 13, EstrinEvaluator> const
k_maclaurin_0_45(std::make_tuple(1.813883936816982644,
0.763163245700557246,
0.761928605321595831,
0.951074653668427927,
1.315180671703161215,

1.928560693477410941,
2.937509342531378755,
4.594894405442878062,
7.330071221881720772,
11.87151259742530180,

19.45851374822937738,
32.20638657246426863,
53.73749198700554656,
90.27388602940998849));
PolynomialInMonomialBasis<double, double, 14, HornerEvaluator> const
k_maclaurin_0_55(std::make_tuple(1.898924910271553526,
0.950521794618244435,
1.151077589959015808,
1.750239106986300540,
2.952676812636875180,

5.285800396121450889,
9.832485716659979747,
18.78714868327559562,
36.61468615273698145,
72.45292395127771801,

145.1079577347069102,
293.4786396308497026,
598.3851815055010179,
1228.420013075863451,
2536.529755382764488));
PolynomialInMonomialBasis<double, double, 16, EstrinEvaluator> const
k_maclaurin_0_65(std::make_tuple(2.007598398424376302,
1.248457231212347337,
1.926234657076479729,
3.751289640087587680,
8.119944554932045802,

18.66572130873555361,
44.60392484291437063,
109.5092054309498377,
274.2779548232413480,
697.5598008606326163,

1795.716014500247129,
4668.381716790389910,
12235.76246813664335,
32290.17809718320818,
85713.07608195964685,

228672.1890493117096,
612757.2711915852774));
PolynomialInMonomialBasis<double, double, 19, EstrinEvaluator> const
k_maclaurin_0_75(std::make_tuple(2.156515647499643235,
1.791805641849463243,
3.826751287465713147,
10.38672468363797208,
31.40331405468070290,

100.9237039498695416,
337.3268282632272897,
1158.707930567827917,
4060.990742193632092,
14454.00184034344795,

52076.66107599404803,
189493.6591462156887,
695184.5762413896145,
2.567994048255284686e6,
9.541921966748386322e6,

3.563492744218076174e7,
1.336692984612040871e8,
5.033521866866284541e8,
1.901975729538660119e9,
7.208915015330103756e9));
PolynomialInMonomialBasis<double, double, 15, EstrinEvaluator> const
k_maclaurin_0_825(std::make_tuple(2.318122621712510589,
2.616920150291232841,
7.897935075731355823,
30.50239715446672327,
131.4869365523528456,

602.9847637356491617,
2877.024617809972641,
14110.51991915180325,
70621.44088156540229,
358977.2665825309926,

1.847238263723971684e6,
9.600515416049214109e6,
5.030767708502366879e7,
2.654441886527127967e8,
1.408862325028702687e9,

7.515687935373774627e9));
PolynomialInMonomialBasis<double, double, 19, EstrinEvaluator> const
k_maclaurin_0_875(std::make_tuple(2.473596173751343912,
3.727624244118099310,
15.60739303554930496,
84.12850842805887747,
506.9818197040613935,

3252.277058145123644,
21713.24241957434256,
149037.0451890932766,
1.043999331089990839e6,
7.427974817042038995e6,

5.350383967558661151e7,
3.892498869948708474e8,
2.855288351100810619e9,
2.109007703876684053e10,
1.566998339477902014e11,

1.170222242422439893e12,
8.777948323668937971e12,
6.610124275248495041e13,
4.994880537133887989e14,
3.785974339724029920e15));

// Double precision general complete elliptic integral "cel"
//
// created by Burlisch
@@ -1435,207 +1623,37 @@ void FukushimaEllipticBDJ(Angle const& φ,
// Inputs: mc = complementary parameter 0 <= mc <= 1
//
double EllipticK(double const mc) {
double m, mx;
double kkc, nome;

constexpr double TINY = 1.0e-99;

m = 1.0 - mc;
if (abs(m) < 1.0e-16) {
// TODO(phl): Use a binary split of [0, 1] to reduce the number of
// comparisons.
double const m = 1.0 - mc;
if (m == 0.0) {
return π / 2;
} else if (mc < TINY) {
return 1.3862943611198906 - 0.5 * log(TINY);
} else if (mc < 1.11e-16) {
return 1.3862943611198906 - 0.5 * log(mc);
} else if (mc < std::numeric_limits<double>::epsilon() / 2.0) {
return 2.0 * log_2 - 0.5 * std::log(mc);
} else if (mc < 0.1) {
nome = EllipticNomeQ<14>(mc);
mx = mc - 0.05;
//
// K'
//
kkc = 1.591003453790792180 + mx * (
0.416000743991786912 + mx * (
0.245791514264103415 + mx * (
0.179481482914906162 + mx * (
0.144556057087555150 + mx * (
0.123200993312427711 + mx * (
0.108938811574293531 + mx * (
0.098853409871592910 + mx * (
0.091439629201749751 + mx * (
0.085842591595413900 + mx * (
0.081541118718303215))))))))));
return -kkc * (1 / π) * log(nome);
double const nome = EllipticNomeQ<14>(mc);
// Evaluate K'.
return -k_maclaurin_0_05.Evaluate(mc - 0.05) * (1 / π) * std::log(nome);
} else if (m <= 0.1) {
mx = m - 0.05;
return 1.591003453790792180 + mx * (
0.416000743991786912 + mx * (
0.245791514264103415 + mx * (
0.179481482914906162 + mx * (
0.144556057087555150 + mx * (
0.123200993312427711 + mx * (
0.108938811574293531 + mx * (
0.098853409871592910 + mx * (
0.091439629201749751 + mx * (
0.085842591595413900 + mx * (
0.081541118718303215))))))))));
return k_maclaurin_0_05.Evaluate(m - 0.05);
} else if (m <= 0.2) {
mx = m - 0.15;
return 1.635256732264579992 + mx * (
0.471190626148732291 + mx * (
0.309728410831499587 + mx * (
0.252208311773135699 + mx * (
0.226725623219684650 + mx * (
0.215774446729585976 + mx * (
0.213108771877348910 + mx * (
0.216029124605188282 + mx * (
0.223255831633057896 + mx * (
0.234180501294209925 + mx * (
0.248557682972264071 + mx * (
0.266363809892617521 + mx * (
0.287728452156114668))))))))))));
return k_maclaurin_0_15.Evaluate(m - 0.15);
} else if (m <= 0.3) {
mx = m - 0.25;
return 1.685750354812596043 + mx * (
0.541731848613280329 + mx * (
0.401524438390690257 + mx * (
0.369642473420889090 + mx * (
0.376060715354583645 + mx * (
0.405235887085125919 + mx * (
0.453294381753999079 + mx * (
0.520518947651184205 + mx * (
0.609426039204995055 + mx * (
0.724263522282908870 + mx * (
0.871013847709812357 + mx * (
1.057652872753547036)))))))))));
return k_maclaurin_0_25.Evaluate(m - 0.25);
} else if (m <= 0.4) {
mx = m - 0.35;
return 1.744350597225613243 + mx * (
0.634864275371935304 + mx * (
0.539842564164445538 + mx * (
0.571892705193787391 + mx * (
0.670295136265406100 + mx * (
0.832586590010977199 + mx * (
1.073857448247933265 + mx * (
1.422091460675497751 + mx * (
1.920387183402304829 + mx * (
2.632552548331654201 + mx * (
3.652109747319039160 + mx * (
5.115867135558865806 + mx * (
7.224080007363877411))))))))))));
return k_maclaurin_0_35.Evaluate(m - 0.35);
} else if (m <= 0.5) {
mx = m - 0.45;
return 1.813883936816982644 + mx * (
0.763163245700557246 + mx * (
0.761928605321595831 + mx * (
0.951074653668427927 + mx * (
1.315180671703161215 + mx * (
1.928560693477410941 + mx * (
2.937509342531378755 + mx * (
4.594894405442878062 + mx * (
7.330071221881720772 + mx * (
11.87151259742530180 + mx * (
19.45851374822937738 + mx * (
32.20638657246426863 + mx * (
53.73749198700554656 + mx * (
90.27388602940998849)))))))))))));
return k_maclaurin_0_45.Evaluate(m - 0.45);
} else if (m <= 0.6) {
mx = m - 0.55;
return 1.898924910271553526 + mx * (
0.950521794618244435 + mx * (
1.151077589959015808 + mx * (
1.750239106986300540 + mx * (
2.952676812636875180 + mx * (
5.285800396121450889 + mx * (
9.832485716659979747 + mx * (
18.78714868327559562 + mx * (
36.61468615273698145 + mx * (
72.45292395127771801 + mx * (
145.1079577347069102 + mx * (
293.4786396308497026 + mx * (
598.3851815055010179 + mx * (
1228.420013075863451 + mx * (
2536.529755382764488))))))))))))));
return k_maclaurin_0_55.Evaluate(m - 0.55);
} else if (m <= 0.7) {
mx = m - 0.65;
return 2.007598398424376302 + mx * (
1.248457231212347337 + mx * (
1.926234657076479729 + mx * (
3.751289640087587680 + mx * (
8.119944554932045802 + mx * (
18.66572130873555361 + mx * (
44.60392484291437063 + mx * (
109.5092054309498377 + mx * (
274.2779548232413480 + mx * (
697.5598008606326163 + mx * (
1795.716014500247129 + mx * (
4668.381716790389910 + mx * (
12235.76246813664335 + mx * (
32290.17809718320818 + mx * (
85713.07608195964685 + mx * (
228672.1890493117096 + mx * (
612757.2711915852774))))))))))))))));
return k_maclaurin_0_65.Evaluate(m - 0.65);
} else if (m <= 0.8) {
mx = m - 0.75;
return 2.156515647499643235 + mx * (
1.791805641849463243 + mx * (
3.826751287465713147 + mx * (
10.38672468363797208 + mx * (
31.40331405468070290 + mx * (
100.9237039498695416 + mx * (
337.3268282632272897 + mx * (
1158.707930567827917 + mx * (
4060.990742193632092 + mx * (
14454.00184034344795 + mx * (
52076.66107599404803 + mx * (
189493.6591462156887 + mx * (
695184.5762413896145 + mx * (
2.567994048255284686e6 + mx * (
9.541921966748386322e6 + mx * (
3.563492744218076174e7 + mx * (
1.336692984612040871e8 + mx * (
5.033521866866284541e8 + mx * (
1.901975729538660119e9 + mx * (
7.208915015330103756e9)))))))))))))))))));
return k_maclaurin_0_75.Evaluate(m - 0.75);
} else if (m <= 0.85) {
mx = m - 0.825;
return 2.318122621712510589 + mx * (
2.616920150291232841 + mx * (
7.897935075731355823 + mx * (
30.50239715446672327 + mx * (
131.4869365523528456 + mx * (
602.9847637356491617 + mx * (
2877.024617809972641 + mx * (
14110.51991915180325 + mx * (
70621.44088156540229 + mx * (
358977.2665825309926 + mx * (
1.847238263723971684e6 + mx * (
9.600515416049214109e6 + mx * (
5.030767708502366879e7 + mx * (
2.654441886527127967e8 + mx * (
1.408862325028702687e9 + mx * (
7.515687935373774627e9)))))))))))))));
return k_maclaurin_0_825.Evaluate(m - 0.825);
} else {
mx = m - 0.875;
return 2.473596173751343912 + mx * (
3.727624244118099310 + mx * (
15.60739303554930496 + mx * (
84.12850842805887747 + mx * (
506.9818197040613935 + mx * (
3252.277058145123644 + mx * (
21713.24241957434256 + mx * (
149037.0451890932766 + mx * (
1.043999331089990839e6 + mx * (
7.427974817042038995e6 + mx * (
5.350383967558661151e7 + mx * (
3.892498869948708474e8 + mx * (
2.855288351100810619e9 + mx * (
2.109007703876684053e10 + mx * (
1.566998339477902014e11 + mx * (
1.170222242422439893e12 + mx * (
8.777948323668937971e12 + mx * (
6.610124275248495041e13 + mx * (
4.994880537133887989e14 + mx * (
3.785974339724029920e15)))))))))))))))))));
return k_maclaurin_0_875.Evaluate(m - 0.875);
}
}

2 changes: 2 additions & 0 deletions quantities/numbers.hpp
Original file line number Diff line number Diff line change
@@ -6,5 +6,7 @@ namespace principia {
constexpr double φ = 1.61803398874989484820458683436563811772030917980576;
constexpr double π = 3.14159265358979323846264338327950288419716939937511;
constexpr double e = 2.71828182845904523536028747135266249775724709369996;
constexpr double log_2 =
6.93147180559945309417232121458176568075500134360255e-1;

} // namespace principia