Coverage for pygeodesy/elliptic.py: 99%

437 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-28 15:52 -0400

1 

2# -*- coding: utf-8 -*- 

3 

4u'''I{Karney}'s elliptic functions and integrals. 

5 

6Class L{Elliptic} transcoded from I{Charles Karney}'s C++ class U{EllipticFunction 

7<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1EllipticFunction.html>} 

8to pure Python, including symmetric integrals L{Elliptic.fRC}, L{Elliptic.fRD}, 

9L{Elliptic.fRF}, L{Elliptic.fRG} and L{Elliptic.fRJ} as C{static methods}. 

10 

11Python method names follow the C++ member functions, I{except}: 

12 

13 - member functions I{without arguments} are mapped to Python properties 

14 prefixed with C{"c"}, for example C{E()} is property C{cE}, 

15 

16 - member functions with 1 or 3 arguments are renamed to Python methods 

17 starting with an C{"f"}, example C{E(psi)} to C{fE(psi)} and C{E(sn, 

18 cn, dn)} to C{fE(sn, cn, dn)}, 

19 

20 - other Python method names conventionally start with a lower-case 

21 letter or an underscore if private. 

22 

23Following is a copy of I{Karney}'s U{EllipticFunction.hpp 

24<https://GeographicLib.SourceForge.io/C++/doc/EllipticFunction_8hpp_source.html>} 

25file C{Header}. 

26 

27Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023) 

28and licensed under the MIT/X11 License. For more information, see the 

29U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. 

30 

31B{Elliptic integrals and functions.} 

32 

33This provides the elliptic functions and integrals needed for 

34C{Ellipsoid}, C{GeodesicExact}, and C{TransverseMercatorExact}. Two 

35categories of function are provided: 

36 

37 - functions to compute U{symmetric elliptic integrals 

38 <https://DLMF.NIST.gov/19.16.i>} 

39 

40 - methods to compute U{Legrendre's elliptic integrals 

41 <https://DLMF.NIST.gov/19.2.ii>} and U{Jacobi elliptic 

42 functions<https://DLMF.NIST.gov/22.2>}. 

43 

44In the latter case, an object is constructed giving the modulus 

45C{k} (and optionally the parameter C{alpha}). The modulus (and 

46parameter) are always passed as squares which allows C{k} to be 

47pure imaginary. (Confusingly, Abramowitz and Stegun call C{m = k**2} 

48the "parameter" and C{n = alpha**2} the "characteristic".) 

49 

50In geodesic applications, it is convenient to separate the incomplete 

51integrals into secular and periodic components, e.g. 

52 

53I{C{E(phi, k) = (2 E(k) / pi) [ phi + delta E(phi, k) ]}} 

54 

55where I{C{delta E(phi, k)}} is an odd periodic function with 

56period I{C{pi}}. 

57 

58The computation of the elliptic integrals uses the algorithms given 

59in U{B. C. Carlson, Computation of real or complex elliptic integrals 

60<https://DOI.org/10.1007/BF02198293>} (also available U{here 

61<https://ArXiv.org/pdf/math/9409227.pdf>}), Numerical Algorithms 10, 

6213--26 (1995) with the additional optimizations given U{here 

63<https://DLMF.NIST.gov/19.36.i>}. 

64 

65The computation of the Jacobi elliptic functions uses the algorithm 

66given in U{R. Bulirsch, Numerical Calculation of Elliptic Integrals 

67and Elliptic Functions<https://DOI.org/10.1007/BF01397975>}, 

68Numerische Mathematik 7, 78--90 (1965). 

69 

70The notation follows U{NIST Digital Library of Mathematical Functions 

71<https://DLMF.NIST.gov>} chapters U{19<https://DLMF.NIST.gov/19>} and 

72U{22<https://DLMF.NIST.gov/22>}. 

73''' 

74# make sure int/int division yields float quotient, see .basics 

75from __future__ import division as _; del _ # PYCHOK semicolon 

76 

77from pygeodesy.basics import copysign0, map2, neg 

78from pygeodesy.constants import EPS, INF, PI, PI_2, PI_4, \ 

79 _EPStol as _TolJAC, _over, \ 

80 _0_0, _0_125, _0_25, _0_5, _1_64th, \ 

81 _1_0, _2_0, _N_2_0, _3_0, _4_0, _6_0, \ 

82 _8_0, _180_0, _360_0 

83from pygeodesy.errors import _ValueError, _xkwds_pop 

84from pygeodesy.fmath import fdot, hypot1 

85from pygeodesy.fsums import Fsum, _sum 

86from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _SPACE_ 

87from pygeodesy.karney import _ALL_LAZY, _signBit 

88# from pygeodesy.lazily import _ALL_LAZY # from .karney 

89from pygeodesy.named import _Named, _NamedTuple 

90from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all 

91from pygeodesy.streprs import Fmt, unstr 

92from pygeodesy.units import Scalar, Scalar_ 

93from pygeodesy.utily import sincos2, sincos2d 

94 

95from math import asinh, atan, atan2, ceil, cosh, fabs, floor, sin, sqrt, tanh 

96 

97__all__ = _ALL_LAZY.elliptic 

98__version__ = '23.08.28' 

99 

100_invokation_ = 'invokation' 

101_TolRD = pow(EPS * 0.002, _0_125) # 8th root: quadquadratic?, qqrt, oqrt? 

102_TolRF = pow(EPS * 0.030, _0_125) # 4th root: biquadratic, bqrt? 

103_TolRG0 = _TolJAC * 2.7 

104_TRIPS = 21 # Max depth, 7 might be sufficient 

105 

106 

107class _CIs(object): 

108 '''(INTERAL) Hold the complete integrals. 

109 ''' 

110 def __init__(self, **kwds): 

111 self.__dict__ = kwds 

112 

113 

114class _Dsum(list): 

115 '''(INTERNAL) Deferred C{Fsum}. 

116 ''' 

117 def __call__(self, s): 

118 return Fsum(*self).fmul(s) 

119 

120 def __iadd__(self, x): 

121 list.append(self, x) 

122 return self 

123 

124 

125class Elliptic(_Named): 

126 '''Elliptic integrals and functions. 

127 

128 @see: I{Karney}'s U{Detailed Description<https://GeographicLib.SourceForge.io/ 

129 html/classGeographicLib_1_1EllipticFunction.html#details>}. 

130 ''' 

131# _alpha2 = 0 

132# _alphap2 = 0 

133# _eps = EPS 

134# _k2 = 0 

135# _kp2 = 0 

136 

137 def __init__(self, k2=0, alpha2=0, kp2=None, alphap2=None, name=NN): 

138 '''Constructor, specifying the C{modulus} and C{parameter}. 

139 

140 @kwarg name: Optional name (C{str}). 

141 

142 @see: Method L{Elliptic.reset} for further details. 

143 

144 @note: If only elliptic integrals of the first and second kinds 

145 are needed, use C{B{alpha2}=0}, the default value. In 

146 that case, we have C{Π(φ, 0, k) = F(φ, k), G(φ, 0, k) = 

147 E(φ, k)} and C{H(φ, 0, k) = F(φ, k) - D(φ, k)}. 

148 ''' 

149 self.reset(k2=k2, alpha2=alpha2, kp2=kp2, alphap2=alphap2) 

150 

151 if name: 

152 self.name = name 

153 

154 @Property_RO 

155 def alpha2(self): 

156 '''Get α^2, the square of the parameter (C{float}). 

157 ''' 

158 return self._alpha2 

159 

160 @Property_RO 

161 def alphap2(self): 

162 '''Get α'^2, the square of the complementary parameter (C{float}). 

163 ''' 

164 return self._alphap2 

165 

166 @Property_RO 

167 def cD(self): 

168 '''Get Jahnke's complete integral C{D(k)} (C{float}), 

169 U{defined<https://DLMF.NIST.gov/19.2.E6>}. 

170 ''' 

171 return self._reset_cDcEcKcKEeps.cD 

172 

173 @Property_RO 

174 def cE(self): 

175 '''Get the complete integral of the second kind C{E(k)} 

176 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E5>}. 

177 ''' 

178 return self._reset_cDcEcKcKEeps.cE 

179 

180 @Property_RO 

181 def cG(self): 

182 '''Get Legendre's complete geodesic longitude integral 

183 C{G(α^2, k)} (C{float}). 

184 ''' 

185 return self._reset_cGcHcPi.cG 

186 

187 @Property_RO 

188 def cH(self): 

189 '''Get Cayley's complete geodesic longitude difference integral 

190 C{H(α^2, k)} (C{float}). 

191 ''' 

192 return self._reset_cGcHcPi.cH 

193 

194 @Property_RO 

195 def cK(self): 

196 '''Get the complete integral of the first kind C{K(k)} 

197 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E4>}. 

198 ''' 

199 return self._reset_cDcEcKcKEeps.cK 

200 

201 @Property_RO 

202 def cKE(self): 

203 '''Get the difference between the complete integrals of the 

204 first and second kinds, C{K(k) − E(k)} (C{float}). 

205 ''' 

206 return self._reset_cDcEcKcKEeps.cKE 

207 

208 @Property_RO 

209 def cPi(self): 

210 '''Get the complete integral of the third kind C{Pi(α^2, k)} 

211 (C{float}), U{defined<https://DLMF.NIST.gov/19.2.E7>}. 

212 ''' 

213 return self._reset_cGcHcPi.cPi 

214 

215 def deltaD(self, sn, cn, dn): 

216 '''The periodic Jahnke's incomplete elliptic integral. 

217 

218 @arg sn: sin(φ). 

219 @arg cn: cos(φ). 

220 @arg dn: sqrt(1 − k2 * sin(2φ)). 

221 

222 @return: Periodic function π D(φ, k) / (2 D(k)) - φ (C{float}). 

223 

224 @raise EllipticError: Invalid invokation or no convergence. 

225 ''' 

226 return _deltaX(sn, cn, dn, self.cD, self.fD) 

227 

228 def deltaE(self, sn, cn, dn): 

229 '''The periodic incomplete integral of the second kind. 

230 

231 @arg sn: sin(φ). 

232 @arg cn: cos(φ). 

233 @arg dn: sqrt(1 − k2 * sin(2φ)). 

234 

235 @return: Periodic function π E(φ, k) / (2 E(k)) - φ (C{float}). 

236 

237 @raise EllipticError: Invalid invokation or no convergence. 

238 ''' 

239 return _deltaX(sn, cn, dn, self.cE, self.fE) 

240 

241 def deltaEinv(self, stau, ctau): 

242 '''The periodic inverse of the incomplete integral of the second kind. 

243 

244 @arg stau: sin(τ) 

245 @arg ctau: cos(τ) 

246 

247 @return: Periodic function E^−1(τ (2 E(k)/π), k) - τ (C{float}). 

248 

249 @raise EllipticError: No convergence. 

250 ''' 

251 # Function is periodic with period pi 

252 t = atan2(-stau, -ctau) if _signBit(ctau) else atan2(stau, ctau) 

253 return self.fEinv(t * self.cE / PI_2) - t 

254 

255 def deltaF(self, sn, cn, dn): 

256 '''The periodic incomplete integral of the first kind. 

257 

258 @arg sn: sin(φ). 

259 @arg cn: cos(φ). 

260 @arg dn: sqrt(1 − k2 * sin(2φ)). 

261 

262 @return: Periodic function π F(φ, k) / (2 K(k)) - φ (C{float}). 

263 

264 @raise EllipticError: Invalid invokation or no convergence. 

265 ''' 

266 return _deltaX(sn, cn, dn, self.cK, self.fF) 

267 

268 def deltaG(self, sn, cn, dn): 

269 '''Legendre's periodic geodesic longitude integral. 

270 

271 @arg sn: sin(φ). 

272 @arg cn: cos(φ). 

273 @arg dn: sqrt(1 − k2 * sin(2φ)). 

274 

275 @return: Periodic function π G(φ, k) / (2 G(k)) - φ (C{float}). 

276 

277 @raise EllipticError: Invalid invokation or no convergence. 

278 ''' 

279 return _deltaX(sn, cn, dn, self.cG, self.fG) 

280 

281 def deltaH(self, sn, cn, dn): 

282 '''Cayley's periodic geodesic longitude difference integral. 

283 

284 @arg sn: sin(φ). 

285 @arg cn: cos(φ). 

286 @arg dn: sqrt(1 − k2 * sin(2φ)). 

287 

288 @return: Periodic function π H(φ, k) / (2 H(k)) - φ (C{float}). 

289 

290 @raise EllipticError: Invalid invokation or no convergence. 

291 ''' 

292 return _deltaX(sn, cn, dn, self.cH, self.fH) 

293 

294 def deltaPi(self, sn, cn, dn): 

295 '''The periodic incomplete integral of the third kind. 

296 

297 @arg sn: sin(φ). 

298 @arg cn: cos(φ). 

299 @arg dn: sqrt(1 − k2 * sin(2φ)). 

300 

301 @return: Periodic function π Π(φ, α2, k) / (2 Π(α2, k)) - φ 

302 (C{float}). 

303 

304 @raise EllipticError: Invalid invokation or no convergence. 

305 ''' 

306 return _deltaX(sn, cn, dn, self.cPi, self.fPi) 

307 

308 @Property_RO 

309 def eps(self): 

310 '''Get epsilon (C{float}). 

311 ''' 

312 return self._reset_cDcEcKcKEeps.eps 

313 

314 def fD(self, phi_or_sn, cn=None, dn=None): 

315 '''Jahnke's incomplete elliptic integral in terms of 

316 Jacobi elliptic functions. 

317 

318 @arg phi_or_sn: φ or sin(φ). 

319 @kwarg cn: C{None} or cos(φ). 

320 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

321 

322 @return: D(φ, k) as though φ ∈ (−π, π] (C{float}), 

323 U{defined<https://DLMF.NIST.gov/19.2.E4>}. 

324 

325 @raise EllipticError: Invalid invokation or no convergence. 

326 ''' 

327 def _fD(sn, cn, dn): 

328 r = fabs(sn)**3 

329 if r: 

330 r = float(_RD(self, cn**2, dn**2, _1_0, _3_0 / r)) 

331 return r 

332 

333 return self._fXf(phi_or_sn, cn, dn, self.cD, 

334 self.deltaD, _fD) 

335 

336 def fDelta(self, sn, cn): 

337 '''The C{Delta} amplitude function. 

338 

339 @arg sn: sin(φ). 

340 @arg cn: cos(φ). 

341 

342 @return: sqrt(1 − k2 * sin(2φ)) (C{float}). 

343 ''' 

344 k2 = self.k2 

345 s = (_1_0 - sn**2 * k2) if k2 < 0 else (self.kp2 

346 + ((cn**2 * k2) if k2 > 0 else _0_0)) 

347 return sqrt(s) if s else _0_0 

348 

349 def fE(self, phi_or_sn, cn=None, dn=None): 

350 '''The incomplete integral of the second kind in terms of 

351 Jacobi elliptic functions. 

352 

353 @arg phi_or_sn: φ or sin(φ). 

354 @kwarg cn: C{None} or cos(φ). 

355 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

356 

357 @return: E(φ, k) as though φ ∈ (−π, π] (C{float}), 

358 U{defined<https://DLMF.NIST.gov/19.2.E5>}. 

359 

360 @raise EllipticError: Invalid invokation or no convergence. 

361 ''' 

362 def _fE(sn, cn, dn): 

363 '''(INTERNAL) Core of C{.fE}. 

364 ''' 

365 if sn: 

366 sn2, cn2, dn2 = sn**2, cn**2, dn**2 

367 kp2, k2 = self.kp2, self.k2 

368 if k2 <= 0: # Carlson, eq. 4.6, <https://DLMF.NIST.gov/19.25.E9> 

369 Ei = _RF3(self, cn2, dn2, _1_0) 

370 if k2: 

371 Ei -= _RD(self, cn2, dn2, _1_0, _3over(k2, sn2)) 

372 elif kp2 >= 0: # k2 > 0, <https://DLMF.NIST.gov/19.25.E10> 

373 Ei = _over(k2 * fabs(cn), dn) # float 

374 if kp2: 

375 Ei += (_RD( self, cn2, _1_0, dn2, _3over(k2, sn2)) + 

376 _RF3(self, cn2, dn2, _1_0)) * kp2 

377 else: # kp2 < 0, <https://DLMF.NIST.gov/19.25.E11> 

378 Ei = _over(dn, fabs(cn)) 

379 Ei -= _RD(self, dn2, _1_0, cn2, _3over(kp2, sn2)) 

380 Ei *= fabs(sn) 

381 ei = float(Ei) 

382 else: # PYCHOK no cover 

383 ei = _0_0 

384 return ei 

385 

386 return self._fXf(phi_or_sn, cn, dn, self.cE, 

387 self.deltaE, _fE) 

388 

389 def fEd(self, deg): 

390 '''The incomplete integral of the second kind with 

391 the argument given in degrees. 

392 

393 @arg deg: Angle (C{degrees}). 

394 

395 @return: E(π B{C{deg}}/180, k) (C{float}). 

396 

397 @raise EllipticError: No convergence. 

398 ''' 

399 if fabs(deg) < _180_0: 

400 e = _0_0 

401 else: # PYCHOK no cover 

402 e = ceil(deg / _360_0 - _0_5) 

403 deg -= e * _360_0 

404 e *= self.cE * _4_0 

405 sn, cn = sincos2d(deg) 

406 return self.fE(sn, cn, self.fDelta(sn, cn)) + e 

407 

408 def fEinv(self, x): 

409 '''The inverse of the incomplete integral of the second kind. 

410 

411 @arg x: Argument (C{float}). 

412 

413 @return: φ = 1 / E(B{C{x}}, k), such that E(φ, k) = B{C{x}} 

414 (C{float}). 

415 

416 @raise EllipticError: No convergence. 

417 ''' 

418 E2 = self.cE * _2_0 

419 n = floor(x / E2 + _0_5) 

420 r = x - E2 * n # r in [-cE, cE) 

421 # linear approximation 

422 phi = PI * r / E2 # phi in [-PI_2, PI_2) 

423 Phi = Fsum(phi) 

424 # first order correction 

425 phi = Phi.fsum_(self.eps * sin(phi * _2_0) / _N_2_0) 

426 # For kp2 close to zero use asin(r / cE) or J. P. Boyd, 

427 # Applied Math. and Computation 218, 7005-7013 (2012) 

428 # <https://DOI.org/10.1016/j.amc.2011.12.021> 

429 _Phi2, self._iteration = Phi.fsum2_, 0 # aggregate 

430 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC 

431 sn, cn, dn = self._sncndnPhi(phi) 

432 if sn: 

433 sn = self.fE(sn, cn, dn) 

434 phi, e = _Phi2((r - sn) / dn) 

435 if fabs(e) < _TolJAC: 

436 _iterations(self, i) 

437 break 

438 else: # PYCHOK no cover 

439 raise _convergenceError(e, _TolJAC, self.fEinv, x) 

440 return Phi.fsum_(n * PI) if n else phi 

441 

442 def fF(self, phi_or_sn, cn=None, dn=None): 

443 '''The incomplete integral of the first kind in terms of 

444 Jacobi elliptic functions. 

445 

446 @arg phi_or_sn: φ or sin(φ). 

447 @kwarg cn: C{None} or cos(φ). 

448 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

449 

450 @return: F(φ, k) as though φ ∈ (−π, π] (C{float}), 

451 U{defined<https://DLMF.NIST.gov/19.2.E4>}. 

452 

453 @raise EllipticError: Invalid invokation or no convergence. 

454 ''' 

455 def _fF(sn, cn, dn): 

456 r = fabs(sn) 

457 if r: 

458 r = float(_RF3(self, cn**2, dn**2, _1_0).fmul(r)) 

459 return r 

460 

461 return self._fXf(phi_or_sn, cn, dn, self.cK, 

462 self.deltaF, _fF) 

463 

464 def fG(self, phi_or_sn, cn=None, dn=None): 

465 '''Legendre's geodesic longitude integral in terms of 

466 Jacobi elliptic functions. 

467 

468 @arg phi_or_sn: φ or sin(φ). 

469 @kwarg cn: C{None} or cos(φ). 

470 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

471 

472 @return: G(φ, k) as though φ ∈ (−π, π] (C{float}). 

473 

474 @raise EllipticError: Invalid invokation or no convergence. 

475 

476 @note: Legendre expresses the longitude of a point on the 

477 geodesic in terms of this combination of elliptic 

478 integrals in U{Exercices de Calcul Intégral, Vol 1 

479 (1811), p 181<https://Books.Google.com/books?id= 

480 riIOAAAAQAAJ&pg=PA181>}. 

481 

482 @see: U{Geodesics in terms of elliptic integrals<https:// 

483 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>} 

484 for the expression for the longitude in terms of this function. 

485 ''' 

486 return self._fXa(phi_or_sn, cn, dn, self.alpha2 - self.k2, 

487 self.cG, self.deltaG) 

488 

489 def fH(self, phi_or_sn, cn=None, dn=None): 

490 '''Cayley's geodesic longitude difference integral in terms of 

491 Jacobi elliptic functions. 

492 

493 @arg phi_or_sn: φ or sin(φ). 

494 @kwarg cn: C{None} or cos(φ). 

495 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

496 

497 @return: H(φ, k) as though φ ∈ (−π, π] (C{float}). 

498 

499 @raise EllipticError: Invalid invokation or no convergence. 

500 

501 @note: Cayley expresses the longitude difference of a point 

502 on the geodesic in terms of this combination of 

503 elliptic integrals in U{Phil. Mag. B{40} (1870), p 333 

504 <https://Books.Google.com/books?id=Zk0wAAAAIAAJ&pg=PA333>}. 

505 

506 @see: U{Geodesics in terms of elliptic integrals<https:// 

507 GeographicLib.SourceForge.io/html/geodesic.html#geodellip>} 

508 for the expression for the longitude in terms of this function. 

509 ''' 

510 return self._fXa(phi_or_sn, cn, dn, -self.alphap2, 

511 self.cH, self.deltaH) 

512 

513 def fPi(self, phi_or_sn, cn=None, dn=None): 

514 '''The incomplete integral of the third kind in terms of 

515 Jacobi elliptic functions. 

516 

517 @arg phi_or_sn: φ or sin(φ). 

518 @kwarg cn: C{None} or cos(φ). 

519 @kwarg dn: C{None} or sqrt(1 − k2 * sin(2φ)). 

520 

521 @return: Π(φ, α2, k) as though φ ∈ (−π, π] (C{float}). 

522 

523 @raise EllipticError: Invalid invokation or no convergence. 

524 ''' 

525 if dn is None and cn is not None: # and isscalar(phi_or_sn) 

526 dn = self.fDelta(phi_or_sn, cn) # in .triaxial 

527 return self._fXa(phi_or_sn, cn, dn, self.alpha2, 

528 self.cPi, self.deltaPi) 

529 

530 def _fXa(self, phi_or_sn, cn, dn, aX, cX, _deltaX): 

531 '''(INTERNAL) Helper for C{.fG}, C{.fH} and C{.fPi}. 

532 ''' 

533 def _fX(sn, cn, dn): 

534 if sn: 

535 cn2, dn2 = cn**2, dn**2 

536 R = _RF3(self, cn2, dn2, _1_0) 

537 if aX: 

538 sn2 = sn**2 

539 p = sn2 * self.alphap2 + cn2 

540 R += _RJ(self, cn2, dn2, _1_0, p, _3over(aX, sn2)) 

541 r = float(R.fmul(fabs(sn))) 

542 else: # PYCHOK no cover 

543 r = _0_0 

544 return r 

545 

546 return self._fXf(phi_or_sn, cn, dn, cX, _deltaX, _fX) 

547 

548 def _fXf(self, phi_or_sn, cn, dn, cX, _deltaX, _fX): 

549 '''(INTERNAL) Helper for C{f.D}, C{.fE}, C{.fF} and C{._fXa}. 

550 ''' 

551 self._iteration = 0 # aggregate 

552 phi = sn = phi_or_sn 

553 if cn is dn is None: # fX(phi) call 

554 sn, cn, dn = self._sncndnPhi(phi) 

555 if fabs(phi) >= PI: 

556 if cX: 

557 cX *= (_deltaX(sn, cn, dn) + phi) / PI_2 

558 return cX 

559 # fall through 

560 elif cn is None or dn is None: 

561 n = NN(_f_, _deltaX.__name__[5:]) 

562 raise _invokationError(n, sn, cn, dn) 

563 

564 if _signBit(cn): # enforce usual trig-like symmetries 

565 xi = cX * _2_0 - _fX(sn, cn, dn) 

566 elif cn > 0: 

567 xi = _fX(sn, cn, dn) 

568 else: 

569 xi = cX 

570 return copysign0(xi, sn) 

571 

572 @Property_RO 

573 def k2(self): 

574 '''Get k^2, the square of the modulus (C{float}). 

575 ''' 

576 return self._k2 

577 

578 @Property_RO 

579 def kp2(self): 

580 '''Get k'^2, the square of the complementary modulus (C{float}). 

581 ''' 

582 return self._kp2 

583 

584 def reset(self, k2=0, alpha2=0, kp2=None, alphap2=None): # MCCABE 13 

585 '''Reset the modulus, parameter and the complementaries. 

586 

587 @kwarg k2: Modulus squared (C{float}, NINF <= k^2 <= 1). 

588 @kwarg alpha2: Parameter squared (C{float}, NINF <= α^2 <= 1). 

589 @kwarg kp2: Complementary modulus squared (C{float}, k'^2 >= 0). 

590 @kwarg alphap2: Complementary parameter squared (C{float}, α'^2 >= 0). 

591 

592 @raise EllipticError: Invalid B{C{k2}}, B{C{alpha2}}, B{C{kp2}} 

593 or B{C{alphap2}}. 

594 

595 @note: The arguments must satisfy C{B{k2} + B{kp2} = 1} and 

596 C{B{alpha2} + B{alphap2} = 1}. No checking is done 

597 that these conditions are met to enable accuracy to be 

598 maintained, e.g., when C{k} is very close to unity. 

599 ''' 

600 if self.__dict__: 

601 _update_all(self, _Named.iteration._uname, Base=Property_RO) 

602 

603 self._k2 = Scalar_(k2=k2, Error=EllipticError, low=None, high=_1_0) 

604 self._kp2 = Scalar_(kp2=((_1_0 - k2) if kp2 is None else kp2), Error=EllipticError) 

605 

606 self._alpha2 = Scalar_(alpha2=alpha2, Error=EllipticError, low=None, high=_1_0) 

607 self._alphap2 = Scalar_(alphap2=((_1_0 - alpha2) if alphap2 is None else alphap2), 

608 Error=EllipticError) 

609 

610 # Values of complete elliptic integrals for k = 0,1 and alpha = 0,1 

611 # K E D 

612 # k = 0: pi/2 pi/2 pi/4 

613 # k = 1: inf 1 inf 

614 # Pi G H 

615 # k = 0, alpha = 0: pi/2 pi/2 pi/4 

616 # k = 1, alpha = 0: inf 1 1 

617 # k = 0, alpha = 1: inf inf pi/2 

618 # k = 1, alpha = 1: inf inf inf 

619 # 

620 # G(0, k) = Pi(0, k) = H(1, k) = E(k) 

621 # H(0, k) = K(k) - D(k) 

622 # Pi(alpha2, 0) = G(alpha2, 0) = pi / (2 * sqrt(1 - alpha2)) 

623 # H( alpha2, 0) = pi / (2 * (sqrt(1 - alpha2) + 1)) 

624 # Pi(alpha2, 1) = inf 

625 # G( alpha2, 1) = H(alpha2, 1) = RC(1, alphap2) 

626 

627 @Property_RO 

628 def _reset_cDcEcKcKEeps(self): 

629 '''(INTERNAL) Get the complete integrals D, E, K and KE plus C{eps}. 

630 ''' 

631 k2 = self.k2 

632 if k2: 

633 kp2 = self.kp2 

634 if kp2: 

635 self._iteration = 0 

636 # D(k) = (K(k) - E(k))/k2, Carlson eq.4.3 

637 # <https://DLMF.NIST.gov/19.25.E1> 

638 D = _RD(self, _0_0, kp2, _1_0, _3_0) 

639 cD = float(D) 

640 # Complete elliptic integral E(k), Carlson eq. 4.2 

641 # <https://DLMF.NIST.gov/19.25.E1> 

642 cE = float(_RG2(self, kp2, _1_0)) 

643 # Complete elliptic integral K(k), Carlson eq. 4.1 

644 # <https://DLMF.NIST.gov/19.25.E1> 

645 cK = _rF2(self, kp2, _1_0) 

646 cKE = float(D.fmul(k2)) 

647 eps = k2 / (sqrt(kp2) + _1_0)**2 

648 else: # PYCHOK no cover 

649 cD = cK = cKE = INF 

650 cE = _1_0 

651 eps = k2 

652 else: # PYCHOK no cover 

653 cD = PI_4 

654 cE = cK = PI_2 

655 cKE = _0_0 # k2 * cD 

656 eps = EPS 

657 

658 return _CIs(cD=cD, cE=cE, cK=cK, cKE=cKE, eps=eps) 

659 

660 @Property_RO 

661 def _reset_cGcHcPi(self): 

662 '''(INTERNAL) Get the complete integrals G, H and Pi. 

663 ''' 

664 self._iteration = 0 

665 alpha2 = self.alpha2 

666 if alpha2: 

667 alphap2 = self.alphap2 

668 if alphap2: 

669 kp2 = self.kp2 

670 if kp2: # <https://DLMF.NIST.gov/19.25.E2> 

671 cK = self.cK 

672 Rj = _RJ(self, _0_0, kp2, _1_0, alphap2, _3_0) 

673 cG = float(Rj * (alpha2 - self.k2) + cK) # G(alpha2, k) 

674 cH = -float(Rj * alphap2 - cK) # H(alpha2, k) 

675 cPi = float(Rj * alpha2 + cK) # Pi(alpha2, k) 

676 else: # PYCHOK no cover 

677 cG = cH = _rC(self, _1_0, alphap2) 

678 cPi = INF # XXX or NAN? 

679 else: # PYCHOK no cover 

680 cG = cH = cPi = INF # XXX or NAN? 

681 else: 

682 cG, cPi, kp2 = self.cE, self.cK, self.kp2 

683 # H = K - D but this involves large cancellations if k2 is near 1. 

684 # So write (for alpha2 = 0) 

685 # H = int(cos(phi)**2/sqrt(1-k2*sin(phi)**2),phi,0,pi/2) 

686 # = 1/sqrt(1-k2) * int(sin(phi)**2/sqrt(1-k2/kp2*sin(phi)**2,...) 

687 # = 1/kp * D(i*k/kp) 

688 # and use D(k) = RD(0, kp2, 1) / 3 

689 # so H = 1/kp * RD(0, 1/kp2, 1) / 3 

690 # = kp2 * RD(0, 1, kp2) / 3 

691 # using <https://DLMF.NIST.gov/19.20.E18>. Equivalently 

692 # RF(x, 1) - RD(0, x, 1)/3 = x * RD(0, 1, x)/3 for x > 0 

693 # For k2 = 1 and alpha2 = 0, we have 

694 # H = int(cos(phi),...) = 1 

695 cH = float(_RD(self, _0_0, _1_0, kp2, _3_0 / kp2)) if kp2 else _1_0 

696 

697 return _CIs(cG=cG, cH=cH, cPi=cPi) 

698 

699 def sncndn(self, x): 

700 '''The Jacobi elliptic function. 

701 

702 @arg x: The argument (C{float}). 

703 

704 @return: An L{Elliptic3Tuple}C{(sn, cn, dn)} with 

705 C{*n(B{x}, k)}. 

706 

707 @raise EllipticError: No convergence. 

708 ''' 

709 self._iteration = 0 # reset 

710 # Bulirsch's sncndn routine, p 89. 

711 if self.kp2: 

712 c, d, cd, mn_ = self._sncndnBulirsch4 

713 dn = _1_0 

714 sn, cn = sincos2(x * cd) 

715 if sn: 

716 a = cn / sn 

717 c *= a 

718 for m, n in mn_: 

719 a *= c 

720 c *= dn 

721 dn = (n + a) / (m + a) 

722 a = c / m 

723 sn = copysign0(_1_0 / hypot1(c), sn) # _signBit(sn) 

724 cn = c * sn 

725 if d and _signBit(self.kp2): # PYCHOK no cover 

726 cn, dn = dn, cn 

727 sn = sn / d # /= chokes PyChecker 

728 else: 

729 sn = tanh(x) 

730 cn = dn = _1_0 / cosh(x) 

731 

732 return Elliptic3Tuple(sn, cn, dn, iteration=self._iteration) 

733 

734 @Property_RO 

735 def _sncndnBulirsch4(self): 

736 '''(INTERNAL) Get Bulirsch' 4-tuple C{(c, d, cd, mn_)}. 

737 ''' 

738 # Bulirsch's sncndn routine, p 89. 

739 d, mc = 0, self.kp2 

740 if _signBit(mc): # PYCHOK no cover 

741 d = _1_0 - mc 

742 mc = neg(mc / d) 

743 d = sqrt(d) 

744 

745 mn, a = [], _1_0 

746 for i in range(1, _TRIPS): # GEOGRAPHICLIB_PANIC 

747 mc = sqrt(mc) 

748 mn.append((a, mc)) 

749 c = (a + mc) * _0_5 

750 t = _TolJAC * a 

751 if fabs(a - mc) <= t: # 6 trips, quadratic 

752 self._iteration += i # accumulate 

753 break 

754 mc *= a 

755 a = c 

756 else: # PYCHOK no cover 

757 raise _convergenceError(a - mc, t, None, kp=self.kp, kp2=self.kp2) 

758 cd = (c * d) if d else c 

759 return c, d, cd, tuple(reversed(mn)) # mn reversed! 

760 

761 def _sncndnPhi(self, phi): 

762 '''(INTERNAL) Helper for C{.fEinv} and C{._fXf}. 

763 ''' 

764 sn, cn = sincos2(phi) 

765 return Elliptic3Tuple(sn, cn, self.fDelta(sn, cn)) 

766 

767 @staticmethod 

768 def fRC(x, y): 

769 '''Degenerate symmetric integral of the first kind C{RC(x, y)}. 

770 

771 @return: C{RC(x, y)}, equivalent to C{RF(x, y, y)}. 

772 

773 @see: U{C{RC} definition<https://DLMF.NIST.gov/19.2.E17>} and 

774 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 

775 ''' 

776 return _rC(None, x, y) 

777 

778 @staticmethod 

779 def fRD(x, y, z, *over): 

780 '''Degenerate symmetric integral of the third kind C{RD(x, y, z)}. 

781 

782 @return: C{RD(x, y, z) / over}, equivalent to C{RJ(x, y, z, z) 

783 / over} with C{over} typically 3. 

784 

785 @see: U{C{RD} definition<https://DLMF.NIST.gov/19.16.E5>} and 

786 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 

787 ''' 

788 return float(_RD(None, x, y, z, *over)) 

789 

790 @staticmethod 

791 def fRF(x, y, z=0): 

792 '''Symmetric or complete symmetric integral of the first kind 

793 C{RF(x, y, z)} respectively C{RF(x, y)}. 

794 

795 @return: C{RF(x, y, z)} or C{RF(x, y)} for missing or zero B{C{z}}. 

796 

797 @see: U{C{RF} definition<https://DLMF.NIST.gov/19.16.E1>} and 

798 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 

799 ''' 

800 return float(_RF3(None, x, y, z)) if z else _rF2(None, x, y) 

801 

802 @staticmethod 

803 def fRG(x, y, z=0): 

804 '''Symmetric or complete symmetric integral of the second kind 

805 C{RG(x, y, z)} respectively C{RG(x, y)}. 

806 

807 @return: C{RG(x, y, z)} or C{RG(x, y)} for missing or zero B{C{z}}. 

808 

809 @see: U{C{RG} definition<https://DLMF.NIST.gov/19.16.E3>} and 

810 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 

811 ''' 

812 G = _RG3(None, x, y, z) if z else (_RG2(None, x, y) * _0_5) 

813 return float(G) 

814 

815 @staticmethod 

816 def fRJ(x, y, z, p): 

817 '''Symmetric integral of the third kind C{RJ(x, y, z, p)}. 

818 

819 @return: C{RJ(x, y, z, p)}. 

820 

821 @see: U{C{RJ} definition<https://DLMF.NIST.gov/19.16.E2>} and 

822 U{Carlson<https://ArXiv.org/pdf/math/9409227.pdf>}. 

823 ''' 

824 return float(_RJ(None, x, y, z, p)) 

825 

826_allPropertiesOf_n(15, Elliptic) # # PYCHOK assert, see Elliptic.reset 

827 

828 

829class EllipticError(_ValueError): 

830 '''Elliptic integral, function, convergence or other L{Elliptic} issue. 

831 ''' 

832 pass 

833 

834 

835class Elliptic3Tuple(_NamedTuple): 

836 '''3-Tuple C{(sn, cn, dn)} all C{scalar}. 

837 ''' 

838 _Names_ = ('sn', 'cn', 'dn') 

839 _Units_ = ( Scalar, Scalar, Scalar) 

840 

841 

842class _Lxyz(list): 

843 '''(INTERNAL) Helper for C{_RD}, C{_RF3} and C{_RJ}. 

844 ''' 

845 _a = None 

846 _a0 = None 

847# _xyzp = () 

848 

849 def __init__(self, *xyzp): # x, y, z [, p] 

850 list.__init__(self, xyzp) 

851 self._xyzp = xyzp 

852 

853 def a0(self, n): 

854 '''Compute the initial C{a}. 

855 ''' 

856 t = tuple(self) 

857 m = n - len(t) 

858 if m > 0: 

859 t += t[-1:] * m 

860 try: 

861 s = Fsum(*t).fover(n) 

862 except ValueError: # Fsum(NAN) exception 

863 s = _sum(t) / n 

864 self._a0 = self._a = s 

865 return s 

866 

867 def asr3(self, a): 

868 '''Compute next C{a}, C{sqrt(xyz_)} and C{fdot(sqrt(xyz))}. 

869 ''' 

870 L = self 

871 # assert a is L._a 

872 s = map2(sqrt, L) # sqrt(x), srqt(y), sqrt(z) [, sqrt(p)] 

873 try: 

874 r = fdot(s[:3], s[1], s[2], s[0]) # sqrt(x) * sqrt(y) + ... 

875 except ValueError: # Fsum(NAN) exception 

876 r = _sum(s[i] * s[(i + 1) % 3] for i in range(3)) 

877 L[:] = [(x + r) * _0_25 for x in L] 

878 # assert L is self 

879 L._a = a = (a + r) * _0_25 

880 return a, s, r 

881 

882 def Casr4(self, inst, n, Tol, where): 

883 '''Yield Carlson 4-tuples C{(m, a, s, r)} plus sentinel, 

884 with C{s = sqrt(xyz_)} and C{r = fdot(sqrt(xyz))}. 

885 ''' 

886 L = self 

887 a = am = L.a0(n) 

888 m = 1 

889 q = max(fabs(a - x) for x in L) / Tol 

890 for i in range(_TRIPS): 

891 if fabs(am) > q: # 5-6 trips 

892 _iterations(inst, i) 

893 break 

894 a, s, r = L.asr3(a) 

895 yield m, a, s, r 

896 m *= 4 

897 am = a * m 

898 else: # PYCHOK no cover 

899 raise _convergenceError(am, q, where, *L._xyzp, thresh=True) 

900 yield m, a, (), None # sentinel: next m, same a, no s, no r 

901 

902 def rescale(self, am, *xy_): 

903 '''Rescale C{x}, C{y}, ... 

904 ''' 

905 for x in xy_: 

906 yield (self._a0 - x) / am 

907 

908 

909def _Cab3(inst, x, y, where): 

910 '''(INTERNAL) Yield C{(i, a, b)} Carlson 3-tuples. 

911 ''' 

912 a, b = sqrt(x), sqrt(y) 

913 if b > a: 

914 a, b = b, a 

915 yield a, b, 0 

916 for i in range(1, _TRIPS): 

917 t = _TolRG0 * a 

918 if fabs(a - b) <= t: # 3-4 trips 

919 _iterations(inst, i - 1) 

920 break 

921 a, b = ((a + b) * _0_5), sqrt(a * b) 

922 yield a, b, i 

923 else: # PYCHOK no cover 

924 raise _convergenceError(a - b, t, where, x, y) 

925 

926 

927def _convergenceError(d, tol, where, *args, **kwds_thresh): # PYCHOK no cover 

928 '''(INTERNAL) Return an L{EllipticError}. 

929 ''' 

930 n = Elliptic.__name__ 

931 if where: 

932 n = _DOT_(n, where.__name__) 

933 if kwds_thresh: 

934 q = _xkwds_pop(kwds_thresh, thresh=False) 

935 t = unstr(n, *args, **kwds_thresh) 

936 else: 

937 q = False 

938 t = unstr(n, *args) 

939 return EllipticError(Fmt.no_convergence(d, tol, thresh=q), txt=t) 

940 

941 

942def _deltaX(sn, cn, dn, cX, _fX): 

943 '''(INTERNAL) Helper for C{Elliptic.deltaD} thru C{.deltaPi}. 

944 ''' 

945 if cn is None or dn is None: 

946 n = NN(_delta_, _fX.__name__[1:]) 

947 raise _invokationError(n, sn, cn, dn) 

948 

949 if _signBit(cn): 

950 cn, sn = -cn, -sn 

951 return _fX(sn, cn, dn) * PI_2 / cX - atan2(sn, cn) 

952 

953 

954def _Horner(S, e1, E2, E3, E4, E5, *over): 

955 '''(INTERNAL) Horner form for C{_RD} and C{_RJ} below. 

956 ''' 

957 E22 = E2**2 

958 # Polynomial is <https://DLMF.NIST.gov/19.36.E2> 

959 # (1 - 3*E2/14 + E3/6 + 9*E2**2/88 - 3*E4/22 - 9*E2*E3/52 

960 # + 3*E5/26 - E2**3/16 + 3*E3**2/40 + 3*E2*E4/20 

961 # + 45*E2**2*E3/272 - 9*(E3*E4+E2*E5)/68) 

962 # converted to Horner-like form ... 

963 F = Fsum 

964 e = e1 * 4084080 

965 S *= e 

966 S += F(E2 * -540540, 471240).fmul(E5) 

967 S += F(E2 * 612612, E3 * -540540, -556920).fmul(E4) 

968 S += F(E2 * -706860, E22 * 675675, E3 * 306306, 680680).fmul(E3) 

969 S += F(E2 * 417690, E22 * -255255, -875160).fmul(E2) 

970 S += 4084080 

971 return S.fdiv((over[0] * e) if over else e) # Fsum 

972 

973 

974def _invokationError(name, *args): # PYCHOK no cover 

975 '''(INTERNAL) Return an L{EllipticError}. 

976 ''' 

977 n = _DOT_(Elliptic.__name__, name) 

978 n = _SPACE_(_invokation_, n) 

979 return EllipticError(NN(n, repr(args))) # unstr 

980 

981 

982def _iterations(inst, i): 

983 '''(INTERNAL) Aggregate iterations B{C{i}}. 

984 ''' 

985 if inst and i > 0: 

986 inst._iteration += i 

987 

988 

989def _3over(a, b): 

990 '''(INTERNAL) Return C{3 / (a * b)}. 

991 ''' 

992 return _over(_3_0, a * b) 

993 

994 

995def _rC(unused, x, y): 

996 '''(INTERNAL) Defined only for y != 0 and x >= 0. 

997 ''' 

998 d = x - y 

999 if d < 0: # catch _NaN 

1000 # <https://DLMF.NIST.gov/19.2.E18> 

1001 d = -d 

1002 r = atan(sqrt(d / x)) if x > 0 else PI_2 

1003 elif d == _0_0: # XXX d < EPS0? or EPS02 or _EPSmin 

1004 d, r = y, _1_0 

1005 elif y > 0: # <https://DLMF.NIST.gov/19.2.E19> 

1006 r = asinh(sqrt(d / y)) # atanh(sqrt((x - y) / x)) 

1007 elif y < 0: # <https://DLMF.NIST.gov/19.2.E20> 

1008 r = asinh(sqrt(-x / y)) # atanh(sqrt(x / (x - y))) 

1009 else: 

1010 raise _invokationError(Elliptic.fRC.__name__, x, y) 

1011 return r / sqrt(d) # float 

1012 

1013 

1014def _RD(inst, x, y, z, *over): 

1015 '''(INTERNAL) Carlson, eqs 2.28 - 2.34. 

1016 ''' 

1017 L = _Lxyz(x, y, z) 

1018 S = _Dsum() 

1019 for m, a, s, r in L.Casr4(inst, 5, _TolRF, Elliptic.fRD): 

1020 if s: 

1021 S += _over(_3_0, (z + r) * s[2] * m) 

1022 z = L[2] # s[2] = sqrt(z) 

1023 x, y = L.rescale(-a * m, x, y) 

1024 xy = x * y 

1025 z = (x + y) / _3_0 

1026 z2 = z**2 

1027 return _Horner(S(_1_0), sqrt(a) * a * m, 

1028 xy - _6_0 * z2, 

1029 (xy * _3_0 - _8_0 * z2) * z, 

1030 (xy - z2) * _3_0 * z2, 

1031 xy * z2 * z, *over) # Fsum 

1032 

1033 

1034def _rF2(inst, x, y): # 2-arg version, z=0 

1035 '''(INTERNAL) Carlson, eqs 2.36 - 2.38. 

1036 ''' 

1037 for a, b, _ in _Cab3(inst, x, y, Elliptic.fRF): # PYCHOK yield 

1038 pass 

1039 return _over(PI, a + b) # float 

1040 

1041 

1042def _RF3(inst, x, y, z): # 3-arg version 

1043 '''(INTERNAL) Carlson, eqs 2.2 - 2.7. 

1044 ''' 

1045 L = _Lxyz(x, y, z) 

1046 for m, a, _, _ in L.Casr4(inst, 3, _TolRF, Elliptic.fRF): 

1047 pass 

1048 x, y = L.rescale(a * m, x, y) 

1049 z = neg(x + y) 

1050 xy = x * y 

1051 e2 = xy - z**2 

1052 e3 = xy * z 

1053 e4 = e2**2 

1054 # Polynomial is <https://DLMF.NIST.gov/19.36.E1> 

1055 # (1 - E2/10 + E3/14 + E2**2/24 - 3*E2*E3/44 

1056 # - 5*E2**3/208 + 3*E3**2/104 + E2**2*E3/16) 

1057 # converted to Horner-like form ... 

1058 S = Fsum(e4 * 15015, e3 * 6930, e2 * -16380, 17160).fmul(e3) 

1059 S += Fsum(e4 * -5775, e2 * 10010, -24024).fmul(e2) 

1060 S += 240240 

1061 return S.fdiv(sqrt(a) * 240240) # Fsum 

1062 

1063 

1064def _RG2(inst, x, y): # 2-args and I{doubled} 

1065 '''(INTERNAL) Carlson, eqs 2.36 - 2.39. 

1066 ''' 

1067 m = -1 # neg! 

1068 S = _Dsum() 

1069 for a, b, i in _Cab3(inst, x, y, Elliptic.fRG): 

1070 if i: 

1071 S += (a - b)**2 * m 

1072 m *= 2 

1073 else: 

1074 S += (a + b)**2 * _0_5 

1075 return S(PI_2 / (a + b)) # Fsum 

1076 

1077 

1078def _RG3(inst, x, y, z): # 3-arg version 

1079 '''(INTERNAL) Never called with zero B{C{z}}, see C{.fRG}. 

1080 ''' 

1081# if not z: 

1082# y, z = z, y 

1083 R = _RF3(inst, x, y, z) 

1084 rd = (x - z) * (z - y) # - (y - z) 

1085 if rd: # Carlson, eq 1.7 

1086 R += _RD(inst, x, y, z, _3_0 * z / rd) 

1087 r = x * y 

1088 if r: 

1089 R += sqrt(r / z**3) 

1090 return R.fmul(z * _0_5) # Fsum 

1091 

1092 

1093def _RJ(inst, x, y, z, p, *over): 

1094 '''(INTERNAL) Carlson, eqs 2.17 - 2.25. 

1095 ''' 

1096 def _xyzp(x, y, z, p): 

1097 return (x + p) * (y + p) * (z + p) 

1098 

1099 L = _Lxyz(x, y, z, p) 

1100 n = neg(_xyzp(x, y, z, -p)) 

1101 S = _Dsum() 

1102 for m, a, s, _ in L.Casr4(inst, 5, _TolRD, Elliptic.fRJ): 

1103 if s: 

1104 d = _xyzp(*s) 

1105 if n: 

1106 rc = _rC(inst, _1_0, n / d**2 + _1_0) 

1107 n *= _1_64th # /= chokes PyChecker 

1108 else: 

1109 rc = _1_0 # == _rC(None, _1_0, _1_0) 

1110 S += rc / (d * m) 

1111 x, y, z = L.rescale(a * m, x, y, z) 

1112 xyz = x * y * z 

1113 p = Fsum(x, y, z).fover(_N_2_0) 

1114 p2 = p**2 

1115 p3 = p2 * p 

1116 E2 = Fsum(x * y, x * z, y * z, -p2 * _3_0) 

1117 E2p = E2 * p 

1118 return _Horner(S(_6_0), sqrt(a) * a * m, E2, 

1119 Fsum(p3 * _4_0, xyz, E2p * _2_0), 

1120 Fsum(p3 * _3_0, E2p, xyz * _2_0).fmul(p), 

1121 xyz * p2, *over) # Fsum 

1122 

1123# **) MIT License 

1124# 

1125# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

1126# 

1127# Permission is hereby granted, free of charge, to any person obtaining a 

1128# copy of this software and associated documentation files (the "Software"), 

1129# to deal in the Software without restriction, including without limitation 

1130# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

1131# and/or sell copies of the Software, and to permit persons to whom the 

1132# Software is furnished to do so, subject to the following conditions: 

1133# 

1134# The above copyright notice and this permission notice shall be included 

1135# in all copies or substantial portions of the Software. 

1136# 

1137# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

1138# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

1139# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

1140# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

1141# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

1142# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

1143# OTHER DEALINGS IN THE SOFTWARE.