1// Copyright ©2017 The Gonum Authors. All rights reserved. 2// Use of this source code is governed by a BSD-style 3// license that can be found in the LICENSE file. 4 5package mathext 6 7import ( 8 "math" 9) 10 11// CompleteK computes the complete elliptic integral of the 1st kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1]. 12// 13// K(m) = \int_{0}^{π/2} 1/{\sqrt{1-m{\sin^2θ}}} dθ 14func CompleteK(m float64) float64 { 15 // Reference: 16 // Toshio Fukushima, Precise and fast computation of complete elliptic integrals 17 // by piecewise minimax rational function approximation, 18 // Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76. 19 // https://doi.org/10.1016/j.cam.2014.12.038 20 // Original Fortran code available at: 21 // https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation 22 if m < 0 || 1 < m || math.IsNaN(m) { 23 return math.NaN() 24 } 25 26 mc := 1 - m 27 28 if mc > 0.592990 { 29 t := 2.45694208987494165*mc - 1.45694208987494165 30 t2 := t * t 31 p := ((3703.75266375099019 + t2*(2744.82029097576810+t2*36.2381612593459565)) + t*(5462.47093231923466+t2*(543.839017382099411+t2*0.393188651542789784))) 32 q := ((2077.94377067058435 + t2*(1959.05960044399275+t2*43.5464368440078942)) + t*(3398.00069767755460+t2*(472.794455487539279+t2))) 33 return p / q 34 } 35 if mc > 0.350756 { 36 t := 4.12823963605439369*mc - 1.44800482178389491 37 t2 := t * t 38 p := ((4264.28203103974630 + t2*(3214.59187442783167+t2*43.2589626155454993)) + t*(6341.90978213264024+t2*(642.790566685354573+t2*0.475223892294445943))) 39 q := ((2125.06914237062279 + t2*(2006.03187933518870+t2*44.1848041560412224)) + t*(3479.95663350926514+t2*(482.900172581418890+t2))) 40 return p / q 41 } 42 if mc > 0.206924 { 43 t := 6.95255575949719117*mc - 1.43865064797819679 44 t2 := t * t 45 p := ((4870.25402224986382 + t2*(3738.29369283392307+t2*51.3609902253065926)) + t*(7307.18826377416591+t2*(754.928587580583704+t2*0.571948962277566451))) 46 q := ((2172.51745704102287 + t2*(2056.13612019430497+t2*44.9026847057686146)) + t*(3565.04737778032566+t2*(493.962405117599400+t2))) 47 return p / q 48 } 49 if mc > 0.121734 { 50 t := 11.7384669562155183*mc - 1.42897053644793990 51 t2 := t * t 52 p := ((5514.8512729127464 + t2*(4313.60788246750934+t2*60.598720224393536)) + t*(8350.4595896779631+t2*(880.27903031894216+t2*0.68504458747933773))) 53 q := ((2218.41682813309737 + t2*(2107.97379949034285+t2*45.6911096775045314)) + t*(3650.41829123846319+t2*(505.74295207655096+t2))) 54 return p / q 55 } 56 if mc > 0.071412 { 57 t := 19.8720241643813839*mc - 1.41910098962680339 58 t2 := t * t 59 p := ((6188.8743957372448 + t2*(4935.41351498551527+t2*70.981049144472361)) + t*(9459.3331440432847+t2*(1018.21910476032105+t2*0.81599895108245948))) 60 q := ((2260.73112539748448 + t2*(2159.68721749761492+t2*46.5298955058476510)) + t*(3732.66955095581621+t2*(517.86964191812384+t2))) 61 return p / q 62 } 63 if mc > 0.041770 { 64 t := 33.7359152553808785*mc - 1.40914918021725929 65 t2 := t * t 66 p := ((6879.5170681289562 + t2*(5594.8381504799829+t2*82.452856129147838)) + t*(10615.0836403687221+t2*(1167.26108955935542+t2*0.96592719058503951))) 67 q := ((2296.88303450660439 + t2*(2208.74949754945558+t2*47.3844470709989137)) + t*(3807.37745652028212+t2*(529.79651353072921+t2))) 68 return p / q 69 } 70 if mc > 0.024360 { 71 t := 57.4382538770821367*mc - 1.39919586444572085 72 t2 := t * t 73 p := ((7570.6827538712100 + t2*(6279.2661370014890+t2*94.886883830605940)) + t*(11792.9392624454532+t2*(1325.01058966228180+t2*1.13537029594409690))) 74 q := ((2324.04824540459984 + t2*(2252.22250562615338+t2*48.2089280211559345)) + t*(3869.56755306385732+t2*(540.85752251676412+t2))) 75 return p / q 76 } 77 if mc > 0.014165 { 78 t := 98.0872976949485042*mc - 1.38940657184894556 79 t2 := t * t 80 p := ((8247.2601660137746 + t2*(6974.7495213178613+t2*108.098282908839979)) + t*(12967.7060124572914+t2*(1488.54008220335966+t2*1.32411616748380686))) 81 q := ((2340.47337508405427 + t2*(2287.70677154700516+t2*48.9575432570382154)) + t*(3915.63324533769906+t2*(550.45072377717361+t2))) 82 return p / q 83 } 84 if mc > 0.008213 { 85 t := 168.010752688172043*mc - 1.37987231182795699 86 t2 := t * t 87 p := ((8894.2961573611293 + t2*(7666.5611739483371+t2*121.863474964652041)) + t*(14113.7038749808951+t2*(1654.60731579994159+t2*1.53112170837206117))) 88 q := ((2344.88618943372377 + t2*(2313.28396270968662+t2*49.5906602613891184)) + t*(3942.81065054556536+t2*(558.07615380622169+t2))) 89 return p / q 90 } 91 if mc > 0 { 92 t := 1.0 - 121.758188238159016*mc 93 p := -math.Log(mc*0.0625) * (34813.4518336350547 + t*(235.767716637974271+t*0.199792723884069485)) / (69483.5736412906324 + t*(614.265044703187382+t)) 94 q := -mc * (9382.53386835986099 + t*(51.6478985993381223+t*0.00410754154682816898)) / (37327.7262507318317 + t*(408.017247271148538+t)) 95 return p + q 96 } 97 98 return math.Inf(1) 99} 100 101// CompleteE computes the complete elliptic integral of the 2nd kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1]. 102// 103// E(m) = \int_{0}^{π/2} {\sqrt{1-m{\sin^2θ}}} dθ 104func CompleteE(m float64) float64 { 105 // Reference: 106 // Toshio Fukushima, Precise and fast computation of complete elliptic integrals 107 // by piecewise minimax rational function approximation, 108 // Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76. 109 // https://doi.org/10.1016/j.cam.2014.12.038 110 // Original Fortran code available at: 111 // https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation 112 if m < 0 || 1 < m || math.IsNaN(m) { 113 return math.NaN() 114 } 115 116 mc := 1 - m 117 118 if mc > 0.566638 { 119 t := 2.30753965506897236*mc - 1.30753965506897236 120 t2 := t * t 121 p := ((19702.2363352671642 + t2*(18177.1879313824040+t2*409.975559128654710)) + t*(31904.1559574281609+t2*(4362.94760768571862+t2*10.3244775335024885))) 122 q := ((14241.2135819448616 + t2*(10266.4884503526076+t2*117.162100771599098)) + t*(20909.9899599927367+t2*(1934.86289070792954+t2))) 123 return p / q 124 } 125 if mc > 0.315153 { 126 t := 3.97638030101198879*mc - 1.25316818100483130 127 t2 := t * t 128 p := ((16317.0721393008221 + t2*(15129.4009798463159+t2*326.113727011739428)) + t*(26627.8852140835023+t2*(3574.15857605556033+t2*7.93163724081373477))) 129 q := ((13047.1505096551210 + t2*(9964.25173735060361+t2*117.670514069579649)) + t*(19753.5762165922376+t2*(1918.72232033637537+t2))) 130 return p / q 131 } 132 if mc > 0.171355 { 133 t := 6.95419964116329852*mc - 1.19163687951153702 134 t2 := t * t 135 p := ((13577.3850240991520 + t2*(12871.9137872656293+t2*263.964361648520708)) + t*(22545.4744699553993+t2*(3000.74575264868572+t2*6.08522443139677663))) 136 q := ((11717.3306408059832 + t2*(9619.40382323874064+t2*118.690522739531267)) + t*(18431.1264424290258+t2*(1904.06010727307491+t2))) 137 return p / q 138 } 139 if mc > 0.090670 { 140 t := 12.3938774245522712*mc - 1.12375286608415443 141 t2 := t * t 142 p := ((11307.9485341543712 + t2*(11208.6068472959372+t2*219.253495956962613)) + t*(19328.6173704569489+t2*(2596.54874477084334+t2*4.66931143174036616))) 143 q := ((10307.6837501971393 + t2*(9241.7604666150102+t2*120.498555754227847)) + t*(16982.2450249024383+t2*(1893.41905403040679+t2))) 144 return p / q 145 } 146 if mc > 0.046453 { 147 t := 22.6157360291290680*mc - 1.05056878576113260 148 t2 := t * t 149 p := ((9383.1490856819874 + t2*(9977.2498973537718+t2*188.618148076418837)) + t*(16718.9730458676860+t2*(2323.49987246555537+t2*3.59313532204509922))) 150 q := ((8877.1964704758383 + t2*(8840.2771293410661+t2*123.422125687316355)) + t*(15450.0537230364062+t2*(1889.13672102820913+t2))) 151 return p / q 152 } 153 if mc > 0.022912 { 154 t := 42.4790790535661187*mc - 0.973280659275306911 155 t2 := t * t 156 p := ((7719.1171817802054 + t2*(9045.3996063894006+t2*169.386557799782496)) + t*(14521.7363804934985+t2*(2149.92068078627829+t2*2.78515570453129137))) 157 q := ((7479.7539074698012 + t2*(8420.3848818926324+t2*127.802109608726363)) + t*(13874.4978011497847+t2*(1892.69753150329759+t2))) 158 return p / q 159 } 160 if mc > 0.010809 { 161 t := 82.6241427745187144*mc - 0.893084359249772784 162 t2 := t * t 163 p := ((6261.6095608987273 + t2*(8304.3265605809870+t2*159.371262600702237)) + t*(12593.0874916293982+t2*(2048.68391263416822+t2*2.18867046462858104))) 164 q := ((6156.4532048239501 + t2*(7979.7435857665227+t2*133.911640385965187)) + t*(12283.8373999680518+t2*(1903.60556312663537+t2))) 165 return p / q 166 } 167 if mc > 0.004841 { 168 t := 167.560321715817694*mc - 0.811159517426273458 169 t2 := t * t 170 p := ((4978.06146583586728 + t2*(7664.6703673290453+t2*156.689647694892782)) + t*(10831.7178150656694+t2*(1995.66437151562090+t2*1.75859085945198570))) 171 q := ((4935.56743322938333 + t2*(7506.8028283118051+t2*141.854303920116856)) + t*(10694.5510113880077+t2*(1918.38517009740321+t2))) 172 return p / q 173 } 174 if mc > 0 { 175 t := 1.0 - 206.568890725056806*mc 176 p := -mc * math.Log(mc*0.0625) * (41566.6612602868736 + t*(154.034981522913482+t*0.0618072471798575991)) / (165964.442527585615 + t*(917.589668642251803+t)) 177 q := (132232.803956682877 + t*(353.375480007017643-t*1.40105837312528026)) / (132393.665743088043 + t*(192.112635228732532-t)) 178 return p + q 179 } 180 181 return 1 182} 183 184// 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]. 185// 186// B(m) = \int_{0}^{π/2} {\cos^2θ} / {\sqrt{1-m{\sin^2θ}}} dθ 187func CompleteB(m float64) float64 { 188 // Reference: 189 // Toshio Fukushima, Precise and fast computation of complete elliptic integrals 190 // by piecewise minimax rational function approximation, 191 // Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76. 192 // https://doi.org/10.1016/j.cam.2014.12.038 193 // Original Fortran code available at: 194 // https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation 195 if m < 0 || 1 < m || math.IsNaN(m) { 196 return math.NaN() 197 } 198 199 mc := 1 - m 200 201 if mc > 0.555073 { 202 t := 2.24755971204264969*mc - 1.24755971204264969 203 t2 := t * t 204 p := ((2030.25011505956379 + t2*(1727.60635612511943+t2*25.0715510300422010)) + t*(3223.16236100954529+t2*(361.164121995173076+t2*0.280355207707726826))) 205 q := ((2420.64907902774675 + t2*(2327.48464880306840+t2*47.9870997057202318)) + t*(4034.28168313496638+t2*(549.234220839203960+t2))) 206 return p / q 207 } 208 if mc > 0.302367 { 209 t := 3.95716761770595079*mc - 1.19651690106289522 210 t2 := t * t 211 p := ((2209.26925068374373 + t2*(1981.37862223307242+t2*29.7612810087709299)) + t*(3606.58475322372526+t2*(422.693774742063054+t2*0.334623999861181980))) 212 q := ((2499.57898767250755 + t2*(2467.63998386656941+t2*50.0198090806651216)) + t*(4236.30953048456334+t2*(581.879599221457589+t2))) 213 return p / q 214 } 215 if mc > 0.161052 { 216 t := 7.07638962601280827*mc - 1.13966670204861480 217 t2 := t * t 218 p := ((2359.14823394150129 + t2*(2254.30785457761760+t2*35.2259786264917876)) + t*(3983.28520266051676+t2*(492.601686517364701+t2*0.396605124984359783))) 219 q := ((2563.95563932625156 + t2*(2633.23323959119935+t2*52.6711647124832948)) + t*(4450.19076667898892+t2*(622.983787815718489+t2))) 220 return p / q 221 } 222 if mc > 0.083522 { 223 t := 12.8982329420869341*mc - 1.07728621178898491 224 t2 := t * t 225 p := ((2464.65334987833736 + t2*(2541.68516994216007+t2*41.5832527504007778)) + t*(4333.38639187691528+t2*(571.53606797524881+t2*0.465975784547025267))) 226 q := ((2600.66956117247726 + t2*(2823.69445052534842+t2*56.136001230010910)) + t*(4661.64381841490914+t2*(674.25435972414302+t2))) 227 return p / q 228 } 229 if mc > 0.041966 { 230 t := 24.0639137549331023*mc - 1.00986620463952257 231 t2 := t * t 232 p := ((2509.86724450741259 + t2*(2835.27071287535469+t2*48.9701196718008345)) + t*(4631.12336462339975+t2*(659.86172161727281+t2*0.54158304771955794))) 233 q := ((2594.15983397593723 + t2*(3034.20118545214106+t2*60.652838995496991)) + t*(4848.17491604384532+t2*(737.15143838356850+t2))) 234 return p / q 235 } 236 if mc > 0.020313 { 237 t := 46.1829769546944996*mc - 0.938114810880709371 238 t2 := t * t 239 p := ((2480.58307884128017 + t2*(3122.00900554841322+t2*57.541132641218839)) + t*(4845.57861173250699+t2*(757.31633816400643+t2*0.62119950515996627))) 240 q := ((2528.85218300581396 + t2*(3253.86151324157460+t2*66.496093157522450)) + t*(4979.31783250484768+t2*(812.40556572486862+t2))) 241 return p / q 242 } 243 if mc > 0.009408 { 244 t := 91.7010545621274645*mc - 0.862723521320495186 245 t2 := t * t 246 p := ((2365.25385348859592 + t2*(3381.09304915246175+t2*67.442026950538221)) + t*(4939.53925884558687+t2*(862.16657576129841+t2*0.70143698925710129))) 247 q := ((2390.48737882063755 + t2*(3462.34808443022907+t2*73.934680452209164)) + t*(5015.4675579215077+t2*(898.99542983710459+t2))) 248 return p / q 249 } 250 if mc > 0.004136 { 251 t := 189.681335356600910*mc - 0.784522003034901366 252 t2 := t * t 253 p := ((2160.82916040868119 + t2*(3584.53058926175721+t2*78.769178005879162)) + t*(4877.14832623847052+t2*(970.53716686804832+t2*0.77797110431753920))) 254 q := ((2172.70451405048305 + t2*(3630.52345460629336+t2*83.173163222639080)) + t*(4916.35263668839769+t2*(993.36676027886685+t2))) 255 return p / q 256 } 257 if mc > 0 { 258 t := 1 - 106.292517006802721*mc 259 p := mc * math.Log(mc*0.0625) * (6607.46457640413908 + t*(19.0287633783211078-t*0.00625368946932704460)) / (26150.3443630974309 + t*(354.603981274536040+t)) 260 q := (26251.5678902584870 + t*(168.788023807915689+t*0.352150236262724288)) / (26065.7912239203873 + t*(353.916840382280456+t)) 261 return p + q 262 } 263 264 return 1 265} 266 267// 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]. 268// 269// D(m) = \int_{0}^{π/2} {\sin^2θ} / {\sqrt{1-m{\sin^2θ}}} dθ 270func CompleteD(m float64) float64 { 271 // Reference: 272 // Toshio Fukushima, Precise and fast computation of complete elliptic integrals 273 // by piecewise minimax rational function approximation, 274 // Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76. 275 // https://doi.org/10.1016/j.cam.2014.12.038 276 // Original Fortran code available at: 277 // https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation 278 if m < 0 || 1 < m || math.IsNaN(m) { 279 return math.NaN() 280 } 281 282 mc := 1 - m 283 284 if mc > 0.599909 { 285 t := 2.49943137936119533*mc - 1.49943137936119533 286 t2 := t * t 287 p := ((1593.39813781813498 + t2*(1058.56241259843217+t2*11.7584241242587571)) + t*(2233.25576544961714+t2*(195.247394601357872+t2*0.101486443490307517))) 288 q := ((1685.47865546030468 + t2*(1604.88100543517015+t2*38.6743012128666717)) + t*(2756.20968383181114+t2*(397.504162950935944+t2))) 289 return p / q 290 } 291 if mc > 0.359180 { 292 t := 4.15404874360795750*mc - 1.49205122772910617 293 t2 := t * t 294 p := ((1967.01442513777287 + t2*(1329.30058268219177+t2*15.0447805948342760)) + t*(2779.87604145516343+t2*(247.475085945854673+t2*0.130547566005491628))) 295 q := ((1749.70634057327467 + t2*(1654.40804288486242+t2*39.1895256017535337)) + t*(2853.92630369567765+t2*(406.925098588378587+t2))) 296 return p / q 297 } 298 if mc > 0.214574 { 299 t := 6.91534237860116454*mc - 1.48385267554596628 300 t2 := t * t 301 p := ((2409.64196912091452 + t2*(1659.30176823041376+t2*19.1942111405094383)) + t*(3436.40744503228691+t2*(312.186468430688790+t2*0.167847673021897479))) 302 q := ((1824.89205701262525 + t2*(1715.38574780156913+t2*39.8798253173462218)) + t*(2971.02216287936566+t2*(418.929791715319490+t2))) 303 return p / q 304 } 305 if mc > 0.127875 { 306 t := 11.5341584101316047*mc - 1.47493050669557896 307 t2 := t * t 308 p := ((2926.81143179637839 + t2*(2056.45624281065334+t2*24.3811986813439843)) + t*(4214.52119721241319+t2*(391.420514384925370+t2*0.215574280659075512))) 309 q := ((1910.33091918583314 + t2*(1787.99942542734799+t2*40.7663012893484449)) + t*(3107.04531802441481+t2*(433.673494280825971+t2))) 310 return p / q 311 } 312 if mc > 0.076007 { 313 t := 19.2797100331611013*mc - 1.46539292049047582 314 t2 := t * t 315 p := ((3520.63614251102960 + t2*(2526.67111759550923+t2*30.7739877519417978)) + t*(5121.2842239226937+t2*(486.926821696342529+t2*0.276315678908126399))) 316 q := ((2003.81997889501324 + t2*(1871.05914195570669+t2*41.8489850490387023)) + t*(3259.09205279874214+t2*(451.007555352632053+t2))) 317 return p / q 318 } 319 if mc > 0.045052 { 320 t := 32.3049588111775157*mc - 1.45540300436116944 321 t2 := t * t 322 p := ((4188.00087087025347 + t2*(3072.05695847158556+t2*38.5070211470790031)) + t*(6156.0080960857764+t2*(599.76666155374012+t2*0.352955925261363680))) 323 q := ((2101.60113938424690 + t2*(1961.76794074710108+t2*43.0997999502743622)) + t*(3421.55151253792527+t2*(470.407158843118117+t2))) 324 return p / q 325 } 326 if mc > 0.026626 { 327 t := 54.2711386084880061*mc - 1.44502333658960165 328 t2 := t * t 329 p := ((4916.74442376570733 + t2*(3688.12811638360551+t2*47.6447145147811350)) + t*(7304.6632479558695+t2*(729.75841970840314+t2*0.448422756936257635))) 330 q := ((2197.49982676612397 + t2*(2055.19657857622715+t2*44.4576261146308645)) + t*(3584.94502590860852+t2*(490.880160668822953+t2))) 331 return p / q 332 } 333 if mc > 0.015689 { 334 t := 91.4327512114839536*mc - 1.43448843375697175 335 t2 := t * t 336 p := ((5688.7542903989517 + t2*(4364.21513060078954+t2*58.159468141567195)) + t*(8542.6096475195826+t2*(875.35992968472914+t2*0.56528145509695951))) 337 q := ((2285.44062680812883 + t2*(2145.80779422696555+t2*45.8427480379028781)) + t*(3739.30422133833258+t2*(511.23253971875808+t2))) 338 return p / q 339 } 340 if mc > 0.009216 { 341 t := 154.487872701992894*mc - 1.42376023482156651 342 t2 := t * t 343 p := ((6475.3392225234969 + t2*(5081.2997108708577+t2*69.910123337464043)) + t*(9829.1138694605662+t2*(1033.32687775311981+t2*0.70526087421186325))) 344 q := ((2357.74885505777295 + t2*(2226.89527217032394+t2*47.1609071069631012)) + t*(3872.32565152553360+t2*(530.03943432061149+t2))) 345 return p / q 346 } 347 if mc > 0 { 348 t := 1 - 108.506944444444444*mc 349 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))) 350 q := -(27356.1090344387530 + t*(107.767403612304371-t*0.0827769227048233593)) / (27104.0854889805978 + t*(358.708172147752755+t)) 351 return p + q 352 } 353 354 return math.Inf(1) 355} 356