github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/ell_complete.go

// Copyright ©2017 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package mathext

import (
	"math"
)

// CompleteK computes the complete elliptic integral of the 1st kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
//
//  K(m) = \int_{0}^{π/2} 1/{\sqrt{1-m{\sin^2θ}}} dθ
func CompleteK(m float64) float64 {
	// Reference:
	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
	// by piecewise minimax rational function approximation,
	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
	// https://doi.org/10.1016/j.cam.2014.12.038
	// Original Fortran code available at:
	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
	if m < 0 || 1 < m || math.IsNaN(m) {
		return math.NaN()
	}

	mc := 1 - m

	if mc > 0.592990 {
		t := 2.45694208987494165*mc - 1.45694208987494165
		t2 := t * t
		p := ((3703.75266375099019 + t2*(2744.82029097576810+t2*36.2381612593459565)) + t*(5462.47093231923466+t2*(543.839017382099411+t2*0.393188651542789784)))
		q := ((2077.94377067058435 + t2*(1959.05960044399275+t2*43.5464368440078942)) + t*(3398.00069767755460+t2*(472.794455487539279+t2)))
		return p / q
	}
	if mc > 0.350756 {
		t := 4.12823963605439369*mc - 1.44800482178389491
		t2 := t * t
		p := ((4264.28203103974630 + t2*(3214.59187442783167+t2*43.2589626155454993)) + t*(6341.90978213264024+t2*(642.790566685354573+t2*0.475223892294445943)))
		q := ((2125.06914237062279 + t2*(2006.03187933518870+t2*44.1848041560412224)) + t*(3479.95663350926514+t2*(482.900172581418890+t2)))
		return p / q
	}
	if mc > 0.206924 {
		t := 6.95255575949719117*mc - 1.43865064797819679
		t2 := t * t
		p := ((4870.25402224986382 + t2*(3738.29369283392307+t2*51.3609902253065926)) + t*(7307.18826377416591+t2*(754.928587580583704+t2*0.571948962277566451)))
		q := ((2172.51745704102287 + t2*(2056.13612019430497+t2*44.9026847057686146)) + t*(3565.04737778032566+t2*(493.962405117599400+t2)))
		return p / q
	}
	if mc > 0.121734 {
		t := 11.7384669562155183*mc - 1.42897053644793990
		t2 := t * t
		p := ((5514.8512729127464 + t2*(4313.60788246750934+t2*60.598720224393536)) + t*(8350.4595896779631+t2*(880.27903031894216+t2*0.68504458747933773)))
		q := ((2218.41682813309737 + t2*(2107.97379949034285+t2*45.6911096775045314)) + t*(3650.41829123846319+t2*(505.74295207655096+t2)))
		return p / q
	}
	if mc > 0.071412 {
		t := 19.8720241643813839*mc - 1.41910098962680339
		t2 := t * t
		p := ((6188.8743957372448 + t2*(4935.41351498551527+t2*70.981049144472361)) + t*(9459.3331440432847+t2*(1018.21910476032105+t2*0.81599895108245948)))
		q := ((2260.73112539748448 + t2*(2159.68721749761492+t2*46.5298955058476510)) + t*(3732.66955095581621+t2*(517.86964191812384+t2)))
		return p / q
	}
	if mc > 0.041770 {
		t := 33.7359152553808785*mc - 1.40914918021725929
		t2 := t * t
		p := ((6879.5170681289562 + t2*(5594.8381504799829+t2*82.452856129147838)) + t*(10615.0836403687221+t2*(1167.26108955935542+t2*0.96592719058503951)))
		q := ((2296.88303450660439 + t2*(2208.74949754945558+t2*47.3844470709989137)) + t*(3807.37745652028212+t2*(529.79651353072921+t2)))
		return p / q
	}
	if mc > 0.024360 {
		t := 57.4382538770821367*mc - 1.39919586444572085
		t2 := t * t
		p := ((7570.6827538712100 + t2*(6279.2661370014890+t2*94.886883830605940)) + t*(11792.9392624454532+t2*(1325.01058966228180+t2*1.13537029594409690)))
		q := ((2324.04824540459984 + t2*(2252.22250562615338+t2*48.2089280211559345)) + t*(3869.56755306385732+t2*(540.85752251676412+t2)))
		return p / q
	}
	if mc > 0.014165 {
		t := 98.0872976949485042*mc - 1.38940657184894556
		t2 := t * t
		p := ((8247.2601660137746 + t2*(6974.7495213178613+t2*108.098282908839979)) + t*(12967.7060124572914+t2*(1488.54008220335966+t2*1.32411616748380686)))
		q := ((2340.47337508405427 + t2*(2287.70677154700516+t2*48.9575432570382154)) + t*(3915.63324533769906+t2*(550.45072377717361+t2)))
		return p / q
	}
	if mc > 0.008213 {
		t := 168.010752688172043*mc - 1.37987231182795699
		t2 := t * t
		p := ((8894.2961573611293 + t2*(7666.5611739483371+t2*121.863474964652041)) + t*(14113.7038749808951+t2*(1654.60731579994159+t2*1.53112170837206117)))
		q := ((2344.88618943372377 + t2*(2313.28396270968662+t2*49.5906602613891184)) + t*(3942.81065054556536+t2*(558.07615380622169+t2)))
		return p / q
	}
	if mc > 0 {
		t := 1.0 - 121.758188238159016*mc
		p := -math.Log(mc*0.0625) * (34813.4518336350547 + t*(235.767716637974271+t*0.199792723884069485)) / (69483.5736412906324 + t*(614.265044703187382+t))
		q := -mc * (9382.53386835986099 + t*(51.6478985993381223+t*0.00410754154682816898)) / (37327.7262507318317 + t*(408.017247271148538+t))
		return p + q
	}

	return math.Inf(1)
}

// CompleteE computes the complete elliptic integral of the 2nd kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
//
//  E(m) = \int_{0}^{π/2} {\sqrt{1-m{\sin^2θ}}} dθ
func CompleteE(m float64) float64 {
	// Reference:
	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
	// by piecewise minimax rational function approximation,
	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
	// https://doi.org/10.1016/j.cam.2014.12.038
	// Original Fortran code available at:
	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
	if m < 0 || 1 < m || math.IsNaN(m) {
		return math.NaN()
	}

	mc := 1 - m

	if mc > 0.566638 {
		t := 2.30753965506897236*mc - 1.30753965506897236
		t2 := t * t
		p := ((19702.2363352671642 + t2*(18177.1879313824040+t2*409.975559128654710)) + t*(31904.1559574281609+t2*(4362.94760768571862+t2*10.3244775335024885)))
		q := ((14241.2135819448616 + t2*(10266.4884503526076+t2*117.162100771599098)) + t*(20909.9899599927367+t2*(1934.86289070792954+t2)))
		return p / q
	}
	if mc > 0.315153 {
		t := 3.97638030101198879*mc - 1.25316818100483130
		t2 := t * t
		p := ((16317.0721393008221 + t2*(15129.4009798463159+t2*326.113727011739428)) + t*(26627.8852140835023+t2*(3574.15857605556033+t2*7.93163724081373477)))
		q := ((13047.1505096551210 + t2*(9964.25173735060361+t2*117.670514069579649)) + t*(19753.5762165922376+t2*(1918.72232033637537+t2)))
		return p / q
	}
	if mc > 0.171355 {
		t := 6.95419964116329852*mc - 1.19163687951153702
		t2 := t * t
		p := ((13577.3850240991520 + t2*(12871.9137872656293+t2*263.964361648520708)) + t*(22545.4744699553993+t2*(3000.74575264868572+t2*6.08522443139677663)))
		q := ((11717.3306408059832 + t2*(9619.40382323874064+t2*118.690522739531267)) + t*(18431.1264424290258+t2*(1904.06010727307491+t2)))
		return p / q
	}
	if mc > 0.090670 {
		t := 12.3938774245522712*mc - 1.12375286608415443
		t2 := t * t
		p := ((11307.9485341543712 + t2*(11208.6068472959372+t2*219.253495956962613)) + t*(19328.6173704569489+t2*(2596.54874477084334+t2*4.66931143174036616)))
		q := ((10307.6837501971393 + t2*(9241.7604666150102+t2*120.498555754227847)) + t*(16982.2450249024383+t2*(1893.41905403040679+t2)))
		return p / q
	}
	if mc > 0.046453 {
		t := 22.6157360291290680*mc - 1.05056878576113260
		t2 := t * t
		p := ((9383.1490856819874 + t2*(9977.2498973537718+t2*188.618148076418837)) + t*(16718.9730458676860+t2*(2323.49987246555537+t2*3.59313532204509922)))
		q := ((8877.1964704758383 + t2*(8840.2771293410661+t2*123.422125687316355)) + t*(15450.0537230364062+t2*(1889.13672102820913+t2)))
		return p / q
	}
	if mc > 0.022912 {
		t := 42.4790790535661187*mc - 0.973280659275306911
		t2 := t * t
		p := ((7719.1171817802054 + t2*(9045.3996063894006+t2*169.386557799782496)) + t*(14521.7363804934985+t2*(2149.92068078627829+t2*2.78515570453129137)))
		q := ((7479.7539074698012 + t2*(8420.3848818926324+t2*127.802109608726363)) + t*(13874.4978011497847+t2*(1892.69753150329759+t2)))
		return p / q
	}
	if mc > 0.010809 {
		t := 82.6241427745187144*mc - 0.893084359249772784
		t2 := t * t
		p := ((6261.6095608987273 + t2*(8304.3265605809870+t2*159.371262600702237)) + t*(12593.0874916293982+t2*(2048.68391263416822+t2*2.18867046462858104)))
		q := ((6156.4532048239501 + t2*(7979.7435857665227+t2*133.911640385965187)) + t*(12283.8373999680518+t2*(1903.60556312663537+t2)))
		return p / q
	}
	if mc > 0.004841 {
		t := 167.560321715817694*mc - 0.811159517426273458
		t2 := t * t
		p := ((4978.06146583586728 + t2*(7664.6703673290453+t2*156.689647694892782)) + t*(10831.7178150656694+t2*(1995.66437151562090+t2*1.75859085945198570)))
		q := ((4935.56743322938333 + t2*(7506.8028283118051+t2*141.854303920116856)) + t*(10694.5510113880077+t2*(1918.38517009740321+t2)))
		return p / q
	}
	if mc > 0 {
		t := 1.0 - 206.568890725056806*mc
		p := -mc * math.Log(mc*0.0625) * (41566.6612602868736 + t*(154.034981522913482+t*0.0618072471798575991)) / (165964.442527585615 + t*(917.589668642251803+t))
		q := (132232.803956682877 + t*(353.375480007017643-t*1.40105837312528026)) / (132393.665743088043 + t*(192.112635228732532-t))
		return p + q
	}

	return 1
}

// CompleteB computes an associate complete elliptic integral of the 2nd kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
//
//  B(m) = \int_{0}^{π/2} {\cos^2θ} / {\sqrt{1-m{\sin^2θ}}} dθ
func CompleteB(m float64) float64 {
	// Reference:
	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
	// by piecewise minimax rational function approximation,
	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
	// https://doi.org/10.1016/j.cam.2014.12.038
	// Original Fortran code available at:
	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
	if m < 0 || 1 < m || math.IsNaN(m) {
		return math.NaN()
	}

	mc := 1 - m

	if mc > 0.555073 {
		t := 2.24755971204264969*mc - 1.24755971204264969
		t2 := t * t
		p := ((2030.25011505956379 + t2*(1727.60635612511943+t2*25.0715510300422010)) + t*(3223.16236100954529+t2*(361.164121995173076+t2*0.280355207707726826)))
		q := ((2420.64907902774675 + t2*(2327.48464880306840+t2*47.9870997057202318)) + t*(4034.28168313496638+t2*(549.234220839203960+t2)))
		return p / q
	}
	if mc > 0.302367 {
		t := 3.95716761770595079*mc - 1.19651690106289522
		t2 := t * t
		p := ((2209.26925068374373 + t2*(1981.37862223307242+t2*29.7612810087709299)) + t*(3606.58475322372526+t2*(422.693774742063054+t2*0.334623999861181980)))
		q := ((2499.57898767250755 + t2*(2467.63998386656941+t2*50.0198090806651216)) + t*(4236.30953048456334+t2*(581.879599221457589+t2)))
		return p / q
	}
	if mc > 0.161052 {
		t := 7.07638962601280827*mc - 1.13966670204861480
		t2 := t * t
		p := ((2359.14823394150129 + t2*(2254.30785457761760+t2*35.2259786264917876)) + t*(3983.28520266051676+t2*(492.601686517364701+t2*0.396605124984359783)))
		q := ((2563.95563932625156 + t2*(2633.23323959119935+t2*52.6711647124832948)) + t*(4450.19076667898892+t2*(622.983787815718489+t2)))
		return p / q
	}
	if mc > 0.083522 {
		t := 12.8982329420869341*mc - 1.07728621178898491
		t2 := t * t
		p := ((2464.65334987833736 + t2*(2541.68516994216007+t2*41.5832527504007778)) + t*(4333.38639187691528+t2*(571.53606797524881+t2*0.465975784547025267)))
		q := ((2600.66956117247726 + t2*(2823.69445052534842+t2*56.136001230010910)) + t*(4661.64381841490914+t2*(674.25435972414302+t2)))
		return p / q
	}
	if mc > 0.041966 {
		t := 24.0639137549331023*mc - 1.00986620463952257
		t2 := t * t
		p := ((2509.86724450741259 + t2*(2835.27071287535469+t2*48.9701196718008345)) + t*(4631.12336462339975+t2*(659.86172161727281+t2*0.54158304771955794)))
		q := ((2594.15983397593723 + t2*(3034.20118545214106+t2*60.652838995496991)) + t*(4848.17491604384532+t2*(737.15143838356850+t2)))
		return p / q
	}
	if mc > 0.020313 {
		t := 46.1829769546944996*mc - 0.938114810880709371
		t2 := t * t
		p := ((2480.58307884128017 + t2*(3122.00900554841322+t2*57.541132641218839)) + t*(4845.57861173250699+t2*(757.31633816400643+t2*0.62119950515996627)))
		q := ((2528.85218300581396 + t2*(3253.86151324157460+t2*66.496093157522450)) + t*(4979.31783250484768+t2*(812.40556572486862+t2)))
		return p / q
	}
	if mc > 0.009408 {
		t := 91.7010545621274645*mc - 0.862723521320495186
		t2 := t * t
		p := ((2365.25385348859592 + t2*(3381.09304915246175+t2*67.442026950538221)) + t*(4939.53925884558687+t2*(862.16657576129841+t2*0.70143698925710129)))
		q := ((2390.48737882063755 + t2*(3462.34808443022907+t2*73.934680452209164)) + t*(5015.4675579215077+t2*(898.99542983710459+t2)))
		return p / q
	}
	if mc > 0.004136 {
		t := 189.681335356600910*mc - 0.784522003034901366
		t2 := t * t
		p := ((2160.82916040868119 + t2*(3584.53058926175721+t2*78.769178005879162)) + t*(4877.14832623847052+t2*(970.53716686804832+t2*0.77797110431753920)))
		q := ((2172.70451405048305 + t2*(3630.52345460629336+t2*83.173163222639080)) + t*(4916.35263668839769+t2*(993.36676027886685+t2)))
		return p / q
	}
	if mc > 0 {
		t := 1 - 106.292517006802721*mc
		p := mc * math.Log(mc*0.0625) * (6607.46457640413908 + t*(19.0287633783211078-t*0.00625368946932704460)) / (26150.3443630974309 + t*(354.603981274536040+t))
		q := (26251.5678902584870 + t*(168.788023807915689+t*0.352150236262724288)) / (26065.7912239203873 + t*(353.916840382280456+t))
		return p + q
	}

	return 1
}

// CompleteD computes an associate complete elliptic integral of the 2nd kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
//
//  D(m) = \int_{0}^{π/2} {\sin^2θ} / {\sqrt{1-m{\sin^2θ}}} dθ
func CompleteD(m float64) float64 {
	// Reference:
	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
	// by piecewise minimax rational function approximation,
	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
	// https://doi.org/10.1016/j.cam.2014.12.038
	// Original Fortran code available at:
	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
	if m < 0 || 1 < m || math.IsNaN(m) {
		return math.NaN()
	}

	mc := 1 - m

	if mc > 0.599909 {
		t := 2.49943137936119533*mc - 1.49943137936119533
		t2 := t * t
		p := ((1593.39813781813498 + t2*(1058.56241259843217+t2*11.7584241242587571)) + t*(2233.25576544961714+t2*(195.247394601357872+t2*0.101486443490307517)))
		q := ((1685.47865546030468 + t2*(1604.88100543517015+t2*38.6743012128666717)) + t*(2756.20968383181114+t2*(397.504162950935944+t2)))
		return p / q
	}
	if mc > 0.359180 {
		t := 4.15404874360795750*mc - 1.49205122772910617
		t2 := t * t
		p := ((1967.01442513777287 + t2*(1329.30058268219177+t2*15.0447805948342760)) + t*(2779.87604145516343+t2*(247.475085945854673+t2*0.130547566005491628)))
		q := ((1749.70634057327467 + t2*(1654.40804288486242+t2*39.1895256017535337)) + t*(2853.92630369567765+t2*(406.925098588378587+t2)))
		return p / q
	}
	if mc > 0.214574 {
		t := 6.91534237860116454*mc - 1.48385267554596628
		t2 := t * t
		p := ((2409.64196912091452 + t2*(1659.30176823041376+t2*19.1942111405094383)) + t*(3436.40744503228691+t2*(312.186468430688790+t2*0.167847673021897479)))
		q := ((1824.89205701262525 + t2*(1715.38574780156913+t2*39.8798253173462218)) + t*(2971.02216287936566+t2*(418.929791715319490+t2)))
		return p / q
	}
	if mc > 0.127875 {
		t := 11.5341584101316047*mc - 1.47493050669557896
		t2 := t * t
		p := ((2926.81143179637839 + t2*(2056.45624281065334+t2*24.3811986813439843)) + t*(4214.52119721241319+t2*(391.420514384925370+t2*0.215574280659075512)))
		q := ((1910.33091918583314 + t2*(1787.99942542734799+t2*40.7663012893484449)) + t*(3107.04531802441481+t2*(433.673494280825971+t2)))
		return p / q
	}
	if mc > 0.076007 {
		t := 19.2797100331611013*mc - 1.46539292049047582
		t2 := t * t
		p := ((3520.63614251102960 + t2*(2526.67111759550923+t2*30.7739877519417978)) + t*(5121.2842239226937+t2*(486.926821696342529+t2*0.276315678908126399)))
		q := ((2003.81997889501324 + t2*(1871.05914195570669+t2*41.8489850490387023)) + t*(3259.09205279874214+t2*(451.007555352632053+t2)))
		return p / q
	}
	if mc > 0.045052 {
		t := 32.3049588111775157*mc - 1.45540300436116944
		t2 := t * t
		p := ((4188.00087087025347 + t2*(3072.05695847158556+t2*38.5070211470790031)) + t*(6156.0080960857764+t2*(599.76666155374012+t2*0.352955925261363680)))
		q := ((2101.60113938424690 + t2*(1961.76794074710108+t2*43.0997999502743622)) + t*(3421.55151253792527+t2*(470.407158843118117+t2)))
		return p / q
	}
	if mc > 0.026626 {
		t := 54.2711386084880061*mc - 1.44502333658960165
		t2 := t * t
		p := ((4916.74442376570733 + t2*(3688.12811638360551+t2*47.6447145147811350)) + t*(7304.6632479558695+t2*(729.75841970840314+t2*0.448422756936257635)))
		q := ((2197.49982676612397 + t2*(2055.19657857622715+t2*44.4576261146308645)) + t*(3584.94502590860852+t2*(490.880160668822953+t2)))
		return p / q
	}
	if mc > 0.015689 {
		t := 91.4327512114839536*mc - 1.43448843375697175
		t2 := t * t
		p := ((5688.7542903989517 + t2*(4364.21513060078954+t2*58.159468141567195)) + t*(8542.6096475195826+t2*(875.35992968472914+t2*0.56528145509695951)))
		q := ((2285.44062680812883 + t2*(2145.80779422696555+t2*45.8427480379028781)) + t*(3739.30422133833258+t2*(511.23253971875808+t2)))
		return p / q
	}
	if mc > 0.009216 {
		t := 154.487872701992894*mc - 1.42376023482156651
		t2 := t * t
		p := ((6475.3392225234969 + t2*(5081.2997108708577+t2*69.910123337464043)) + t*(9829.1138694605662+t2*(1033.32687775311981+t2*0.70526087421186325)))
		q := ((2357.74885505777295 + t2*(2226.89527217032394+t2*47.1609071069631012)) + t*(3872.32565152553360+t2*(530.03943432061149+t2)))
		return p / q
	}
	if mc > 0 {
		t := 1 - 108.506944444444444*mc
		p := -math.Log(mc*0.0625) * (6.2904323649908115e6 + t*(58565.284164780476+t*(131.176674599188545+t*0.0426826410911220304))) / (1.24937550257219890e7 + t*(203580.534005225410+t*(921.17729845011868+t)))
		q := -(27356.1090344387530 + t*(107.767403612304371-t*0.0827769227048233593)) / (27104.0854889805978 + t*(358.708172147752755+t))
		return p + q
	}

	return math.Inf(1)
}