1

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

3

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

5

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

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

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

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

10

11 Python 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

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

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

25 file C{Header}.

26

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

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

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

30

31 B{Elliptic integrals and functions.}

32

33 This provides the elliptic functions and integrals needed for

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

35 categories 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

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

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

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

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

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

49

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

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

52

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

54

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

56 period I{C{pi}}.

57

58 The computation of the elliptic integrals uses the algorithms given

59 in 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,

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

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

64

65 The computation of the Jacobi elliptic functions uses the algorithm

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

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

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

69

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

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

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

73'''

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

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

76

77 from pygeodesy.basics import copysign0, map2, neg, neg_

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

79 _EPStol as _TolJAC, _0_0, _1_64th, \

80 _0_25, _0_5, _1_0, _2_0, _N_2_0, \

81 _3_0, _4_0, _6_0, _8_0, _180_0, \

82 _360_0, _over

83# from pygeodesy.errors import _ValueError # from .fmath

84 from pygeodesy.fmath import fdot, hypot1, zqrt, _ValueError

85 from pygeodesy.fsums import Fsum, _sum

86# from pygeodesy.internals import _dunder_nameof # from .lazily

87 from pygeodesy.interns import NN, _delta_, _DOT_, _f_, _invalid_, \

88 _invokation_, _negative_, _SPACE_

89 from pygeodesy.karney import _K_2_0, _norm180, _signBit, _sincos2

90 from pygeodesy.lazily import _ALL_LAZY, _dunder_nameof

91 from pygeodesy.named import _Named, _NamedTuple, Fmt, unstr

92 from pygeodesy.props import _allPropertiesOf_n, Property_RO, _update_all

93# from pygeodesy.streprs import Fmt, unstr # from .named

94 from pygeodesy.units import Scalar, Scalar_

95# from pygeodesy.utily import sincos2 as _sincos2 # from .karney

96

97 from math import asinh, atan, atan2, ceil, cosh, fabs, floor, \

98 radians, sin, sqrt, tanh

99

100__all__ = _ALL_LAZY.elliptic

101__version__ = '24.05.13'

102

103_TolRD = zqrt(EPS * 0.002)

104_TolRF = zqrt(EPS * 0.030)

105_TolRG0 = _TolJAC * 2.7

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

107

108

109 class _Cs(object):

110 '''(INTERAL) Complete integrals cache.

111 '''

112 def __init__(self, **kwds):

113 self.__dict__ = kwds

114

115

116 class _Dsum(list):

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

118 '''

119 def __call__(self, s):

120 try: # Fsum *= s

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

122 except ValueError: # Fsum(NAN) exception

123 return _sum(self) * s

124

125 def __iadd__(self, x):

126 list.append(self, x)

127 return self

128

129

130 class Elliptic(_Named):

131 '''Elliptic integrals and functions.

132

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

134 C++/doc/classGeographicLib_1_1EllipticFunction.html#details>}.

135 '''

136# _alpha2 = 0

137# _alphap2 = 0

138# _eps = EPS

139# _k2 = 0

140# _kp2 = 0

141

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

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

144

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

146

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

148

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

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

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

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

153 '''

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

155

156 if name:

157 self.name = name

158

159 @Property_RO

160 def alpha2(self):

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

162 '''

163 return self._alpha2

164

165 @Property_RO

166 def alphap2(self):

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

168 '''

169 return self._alphap2

170

171 @Property_RO

172 def cD(self):

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

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

175 '''

176 return self._cDEKEeps.cD

177

178 @Property_RO

179 def _cDEKEeps(self):

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

181 '''

182 k2, kp2 = self.k2, self.kp2

183 if k2:

184 if kp2:

185 try:

186 self._iteration = 0

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

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

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

190 cD = float(D)

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

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

193 cE = _rG2(self, kp2, _1_0, PI_=PI_2)

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

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

196 cK = _rF2(self, kp2, _1_0)

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

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

199

200 except Exception as e:

201 raise _ellipticError(self.reset, k2=k2, kp2=kp2, cause=e)

202 else:

203 cD = cK = cKE = INF

204 cE = _1_0

205 eps = k2

206 else:

207 cD = PI_4

208 cE = cK = PI_2

209 cKE = _0_0 # k2 * cD

210 eps = EPS

211

212 return _Cs(cD=cD, cE=cE, cK=cK, cKE=cKE, eps=eps)

213

214 @Property_RO

215 def cE(self):

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

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

218 '''

219 return self._cDEKEeps.cE

220

221 @Property_RO

222 def cG(self):

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

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

225 '''

226 return self._cGHPi.cG

227

228 @Property_RO

229 def _cGHPi(self):

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

231 '''

232 alpha2, alphap2, kp2 = self.alpha2, self.alphap2, self.kp2

233 try:

234 self._iteration = 0

235 if alpha2:

236 if alphap2:

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

238 cK = self.cK

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

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

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

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

243 else:

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

245 cPi = INF # XXX or NAN?

246 else:

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

248 else:

249 cG, cPi = self.cE, self.cK

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

251 # So write (for alpha2 = 0)

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

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

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

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

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

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

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

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

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

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

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

263

264 except Exception as e:

265 raise _ellipticError(self.reset, kp2=kp2, alpha2 =alpha2,

266 alphap2=alphap2, cause=e)

267 return _Cs(cG=cG, cH=cH, cPi=cPi)

268

269 @Property_RO

270 def cH(self):

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

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

273 '''

274 return self._cGHPi.cH

275

276 @Property_RO

277 def cK(self):

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

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

280 '''

281 return self._cDEKEeps.cK

282

283 @Property_RO

284 def cKE(self):

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

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

287 '''

288 return self._cDEKEeps.cKE

289

290 @Property_RO

291 def cPi(self):

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

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

294 '''

295 return self._cGHPi.cPi

296

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

298 '''Jahnke's periodic incomplete elliptic integral.

299

300 @arg sn: sin(φ).

301 @arg cn: cos(φ).

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

303

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

305

306 @raise EllipticError: Invalid invokation or no convergence.

307 '''

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

309

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

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

312

313 @arg sn: sin(φ).

314 @arg cn: cos(φ).

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

316

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

318

319 @raise EllipticError: Invalid invokation or no convergence.

320 '''

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

322

323 def deltaEinv(self, stau, ctau):

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

325

326 @arg stau: sin(τ)

327 @arg ctau: cos(τ)

328

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

330

331 @raise EllipticError: No convergence.

332 '''

333 try:

334 if _signBit(ctau): # pi periodic

335 stau, ctau = neg_(stau, ctau)

336 t = atan2(stau, ctau)

337 return self._Einv(t * self.cE / PI_2) - t

338

339 except Exception as e:

340 raise _ellipticError(self.deltaEinv, stau, ctau, cause=e)

341

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

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

344

345 @arg sn: sin(φ).

346 @arg cn: cos(φ).

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

348

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

350

351 @raise EllipticError: Invalid invokation or no convergence.

352 '''

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

354

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

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

357

358 @arg sn: sin(φ).

359 @arg cn: cos(φ).

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

361

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

363

364 @raise EllipticError: Invalid invokation or no convergence.

365 '''

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

367

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

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

370

371 @arg sn: sin(φ).

372 @arg cn: cos(φ).

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

374

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

376

377 @raise EllipticError: Invalid invokation or no convergence.

378 '''

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

380

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

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

383

384 @arg sn: sin(φ).

385 @arg cn: cos(φ).

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

387

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

389 (C{float}).

390

391 @raise EllipticError: Invalid invokation or no convergence.

392 '''

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

394

395 def _Einv(self, x):

396 '''(INTERNAL) Helper for C{.deltaEinv} and C{.fEinv}.

397 '''

398 E2 = self.cE * _2_0

399 n = floor(x / E2 + _0_5)

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

401 # linear approximation

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

403 Phi = Fsum(phi)

404 # first order correction

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

406 self._iteration = 0

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

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

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

410 _Phi2 = Phi.fsum2f_ # aggregate

411 _abs = fabs

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

413 sn, cn, dn = self._sncndn3(phi)

414 if dn:

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

416 phi, d = _Phi2((r - sn) / dn)

417 else: # PYCHOK no cover

418 d = _0_0 # XXX continue?

419 if _abs(d) < _TolJAC: # 3-4 trips

420 _iterations(self, i)

421 break

422 else: # PYCHOK no cover

423 raise _convergenceError(d, _TolJAC)

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

425

426 @Property_RO

427 def eps(self):

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

429 '''

430 return self._cDEKEeps.eps

431

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

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

434 Jacobi elliptic functions.

435

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

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

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

439

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

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

442

443 @raise EllipticError: Invalid invokation or no convergence.

444 '''

445 def _fD(sn, cn, dn):

446 r = fabs(sn)**3

447 if r:

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

449 return r

450

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

452 self.deltaD, _fD)

453

454 def fDelta(self, sn, cn):

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

456

457 @arg sn: sin(φ).

458 @arg cn: cos(φ).

459

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

461 '''

462 try:

463 k2 = self.k2

464 s = (self.kp2 + cn**2 * k2) if k2 > 0 else (

465 (_1_0 - sn**2 * k2) if k2 < 0 else self.kp2)

466 return sqrt(s) if s else _0_0

467

468 except Exception as e:

469 raise _ellipticError(self.fDelta, sn, cn, k2=k2, cause=e)

470

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

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

473 Jacobi elliptic functions.

474

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

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

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

478

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

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

481

482 @raise EllipticError: Invalid invokation or no convergence.

483 '''

484 def _fE(sn, cn, dn):

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

486 '''

487 if sn:

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

489 kp2, k2 = self.kp2, self.k2

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

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

492 if k2:

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

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

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

496 if kp2:

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

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

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

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

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

502 Ei *= fabs(sn)

503 ei = float(Ei)

504 else: # PYCHOK no cover

505 ei = _0_0

506 return ei

507

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

509 self.deltaE, _fE)

510

511 def fEd(self, deg):

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

513 the argument given in C{degrees}.

514

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

516

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

518

519 @raise EllipticError: No convergence.

520 '''

521 if _K_2_0:

522 e = round((deg - _norm180(deg)) / _360_0)

523 elif fabs(deg) < _180_0:

524 e = _0_0

525 else:

526 e = ceil(deg / _360_0 - _0_5)

527 deg -= e * _360_0

528 return self.fE(radians(deg)) + e * self.cE * _4_0

529

530 def fEinv(self, x):

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

532

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

534

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

536 (C{float}).

537

538 @raise EllipticError: No convergence.

539 '''

540 try:

541 return self._Einv(x)

542 except Exception as e:

543 raise _ellipticError(self.fEinv, x, cause=e)

544

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

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

547 Jacobi elliptic functions.

548

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

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

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

552

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

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

555

556 @raise EllipticError: Invalid invokation or no convergence.

557 '''

558 def _fF(sn, cn, dn):

559 r = fabs(sn)

560 if r:

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

562 return r

563

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

565 self.deltaF, _fF)

566

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

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

569 Jacobi elliptic functions.

570

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

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

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

574

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

576

577 @raise EllipticError: Invalid invokation or no convergence.

578

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

580 geodesic in terms of this combination of elliptic

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

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

583 riIOAAAAQAAJ&pg=PA181>}.

584

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

586 GeographicLib.SourceForge.io/C++/doc/geodesic.html#geodellip>}

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

588 '''

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

590 self.cG, self.deltaG)

591

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

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

594 Jacobi elliptic functions.

595

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

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

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

599

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

601

602 @raise EllipticError: Invalid invokation or no convergence.

603

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

605 on the geodesic in terms of this combination of

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

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

608

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

610 GeographicLib.SourceForge.io/C++/doc/geodesic.html#geodellip>}

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

612 '''

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

614 self.cH, self.deltaH)

615

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

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

618 Jacobi elliptic functions.

619

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

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

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

623

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

625

626 @raise EllipticError: Invalid invokation or no convergence.

627 '''

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

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

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

631 self.cPi, self.deltaPi)

632

633 def _fXa(self, phi_or_sn, cn, dn, aX, cX, deltaX):

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

635 '''

636 def _fX(sn, cn, dn):

637 if sn:

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

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

640 if aX:

641 sn2 = sn**2

642 p = sn2 * self.alphap2 + cn2

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

644 R *= fabs(sn)

645 r = float(R)

646 else: # PYCHOK no cover

647 r = _0_0

648 return r

649

650 return self._fXf(phi_or_sn, cn, dn, cX, deltaX, _fX)

651

652 def _fXf(self, phi_or_sn, cn, dn, cX, deltaX, fX):

653 '''(INTERNAL) Helper for C{.fD}, C{.fE}, C{.fF} and C{._fXa}.

654 '''

655 self._iteration = 0 # aggregate

656 phi = sn = phi_or_sn

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

658 sn, cn, dn = self._sncndn3(phi)

659 if fabs(phi) >= PI:

660 return (deltaX(sn, cn, dn) + phi) * cX / PI_2

661 # fall through

662 elif cn is None or dn is None:

663 n = NN(_f_, deltaX.__name__[5:])

664 raise _ellipticError(n, sn, cn, dn)

665

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

667 xi = cX * _2_0 - fX(sn, cn, dn)

668 else:

669 xi = fX(sn, cn, dn) if cn > 0 else cX

670 return copysign0(xi, sn)

671

672 @Property_RO

673 def k2(self):

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

675 '''

676 return self._k2

677

678 @Property_RO

679 def kp2(self):

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

681 '''

682 return self._kp2

683

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

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

686

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

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

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

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

691

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

693 or B{C{alphap2}}.

694

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

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

697 that these conditions are met to enable accuracy to be

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

699 '''

700 if self.__dict__:

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

702

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

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

705

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

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

708 Error=EllipticError)

709

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

711 # K E D

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

713 # k = 1: inf 1 inf

714 # Pi G H

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

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

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

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

719 #

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

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

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

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

724 # Pi(alpha2, 1) = inf

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

726

727 def sncndn(self, x):

728 '''The Jacobi elliptic function.

729

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

731

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

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

734

735 @raise EllipticError: No convergence.

736 '''

737 self._iteration = 0 # reset

738 try: # Bulirsch's sncndn routine, p 89.

739 if self.kp2:

740 c, d, cd, mn = self._sncndn4

741 dn = _1_0

742 sn, cn = _sincos2(x * cd)

743 if sn:

744 a = cn / sn

745 c *= a

746 for m, n in reversed(mn):

747 a *= c

748 c *= dn

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

750 a = c / m

751 a = _1_0 / hypot1(c)

752 sn = neg(a) if _signBit(sn) else a

753 cn = c * sn

754 if d and _signBit(self.kp2):

755 cn, dn = dn, cn

756 sn = sn / d # /= chokes PyChecker

757 else:

758 sn = tanh(x)

759 cn = dn = _1_0 / cosh(x)

760

761 except Exception as e:

762 raise _ellipticError(self.sncndn, x, kp2=self.kp2, cause=e)

763

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

765

766 def _sncndn3(self, phi):

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

768 '''

769 sn, cn = _sincos2(phi)

770 return sn, cn, self.fDelta(sn, cn)

771

772 @Property_RO

773 def _sncndn4(self):

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

775 '''

776 # Bulirsch's sncndn routine, p 89.

777 d, mc = 0, self.kp2

778 if _signBit(mc):

779 d = _1_0 - mc

780 mc = neg(mc / d)

781 d = sqrt(d)

782

783 mn, a = [], _1_0

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

785 mc = sqrt(mc)

786 mn.append((a, mc))

787 c = (a + mc) * _0_5

788 r = fabs(mc - a)

789 t = _TolJAC * a

790 if r <= t: # 6 trips, quadratic

791 _iterations(self, i)

792 break

793 mc *= a

794 a = c

795 else: # PYCHOK no cover

796 raise _convergenceError(r, t)

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

798 return c, d, cd, mn

799

800 @staticmethod

801 def fRC(x, y):

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

803

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

805

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

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

808 '''

809 return _rC(None, x, y)

810

811 @staticmethod

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

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

814

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

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

817

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

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

820 '''

821 try:

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

823 except Exception as e:

824 raise _ellipticError(Elliptic.fRD, x, y, z, *over, cause=e)

825

826 @staticmethod

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

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

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

830

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

832

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

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

835 '''

836 try:

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

838 except Exception as e:

839 raise _ellipticError(Elliptic.fRF, x, y, z, cause=e)

840

841 @staticmethod

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

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

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

845

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

847

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

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

850 U{RG<https://GeographicLib.SourceForge.io/C++/doc/

851 EllipticFunction_8cpp_source.html#l00096>} version 2.3.

852 '''

853 try:

854 return _rG2(None, x, y) if z == 0 else (

855 _rG2(None, z, x) if y == 0 else (

856 _rG2(None, y, z) if x == 0 else _rG3(None, x, y, z)))

857 except Exception as e:

858 t = _negative_ if min(x, y, z) < 0 else NN

859 raise _ellipticError(Elliptic.fRG, x, y, z, cause=e, txt=t)

860

861 @staticmethod

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

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

864

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

866

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

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

869 '''

870 try:

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

872 except Exception as e:

873 raise _ellipticError(Elliptic.fRJ, x, y, z, p, cause=e)

874

875 @staticmethod

876 def _RFRD(x, y, z, m):

877 # in .auxilats.AuxDLat.DE, .auxilats.AuxLat.Rectifying

878 try: # float(RF(x, y, z) - RD(x, y, z, 3 / m))

879 R = _RF3(None, x, y, z)

880 if m:

881 R -= _RD(None, x, y, z, _3_0 / m)

882 except Exception as e:

883 raise _ellipticError(Elliptic._RFRD, x, y, z, m, cause=e)

884 return float(R)

885

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

887

888

889 class EllipticError(_ValueError):

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

891 '''

892 pass

893

894

895 class Elliptic3Tuple(_NamedTuple):

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

897 '''

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

899 _Units_ = ( Scalar, Scalar, Scalar)

900

901

902 class _List(list):

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

904 '''

905 _a0 = None

906# _xyzp = ()

907

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

909 list.__init__(self, xyzp)

910 self._xyzp = xyzp

911

912 def a0(self, n):

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

914 '''

915 t = tuple(self)

916 m = n - len(t)

917 if m > 0:

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

919 try:

920 a = Fsum(*t).fover(n)

921 except ValueError: # Fsum(NAN) exception

922 a = _sum(t) / n

923 self._a0 = a

924 return a

925

926 def amrs4(self, inst, y, Tol):

927 '''Yield Carlson 4-tuples C{(An, mul, lam, s)} plus sentinel, with

928 C{lam = fdot(sqrt(x), ... (z))} and C{s = (sqrt(x), ... (p))}.

929 '''

930 L = self

931 a = L.a0(5 if y else 3)

932 m = 1

933 t = max(fabs(a - _) for _ in L) / Tol

934 for i in range(_TRIPS):

935 d = fabs(a * m)

936 if d > t: # 3-6 trips

937 _iterations(inst, i)

938 break

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

940 try:

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

942 except ValueError: # Fsum(NAN) exception

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

944 L[:] = ((r + _) * _0_25 for _ in L)

945 a = (r + a) * _0_25

946 if y: # yield only if used

947 yield a, m, r, s # L[2] is next z

948 m *= 4

949 else: # PYCHOK no cover

950 raise _convergenceError(d, t, thresh=True)

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

952

953 def rescale(self, am, *xs):

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

955 '''

956 # assert am

957 a0 = self._a0

958 _am = _1_0 / am

959 for x in xs:

960 yield (a0 - x) * _am

961

962

963 def _ab2(inst, x, y):

964 '''(INTERNAL) Yield Carlson 2-tuples C{(xn, yn)}.

965 '''

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

967 if b > a:

968 b, a = a, b

969 for i in range(_TRIPS):

970 yield a, b # xi, yi

971 d = fabs(a - b)

972 t = _TolRG0 * a

973 if d <= t: # 3-4 trips

974 _iterations(inst, i)

975 break

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

977 else: # PYCHOK no cover

978 raise _convergenceError(d, t)

979

980

981 def _convergenceError(d, tol, **thresh):

982 '''(INTERNAL) Format a no-convergence Error.

983 '''

984 t = Fmt.no_convergence(d, tol, **thresh)

985 return ValueError(t) # txt only

986

987

988 def _deltaX(sn, cn, dn, cX, fX):

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

990 '''

991 try:

992 if cn is None or dn is None:

993 raise ValueError(_invalid_)

994

995 if _signBit(cn):

996 sn, cn = neg_(sn, cn)

997 r = fX(sn, cn, dn) * PI_2 / cX

998 return r - atan2(sn, cn)

999

1000 except Exception as e:

1001 n = NN(_delta_, fX.__name__[1:])

1002 raise _ellipticError(n, sn, cn, dn, cause=e)

1003

1004

1005 def _ellipticError(where, *args, **kwds_cause_txt):

1006 '''(INTERNAL) Format an L{EllipticError}.

1007 '''

1008 def _x_t_kwds(cause=None, txt=NN, **kwds):

1009 return cause, txt, kwds

1010

1011 x, t, kwds = _x_t_kwds(**kwds_cause_txt)

1012

1013 n = _dunder_nameof(where, where)

1014 n = _DOT_(Elliptic.__name__, n)

1015 n = _SPACE_(_invokation_, n)

1016 u = unstr(n, *args, **kwds)

1017 return EllipticError(u, cause=x, txt=t)

1018

1019

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

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

1022 '''

1023 E22 = E2**2

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

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

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

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

1028 # converted to Horner-like form ...

1029 e = e1 * 4084080

1030 S *= e

1031 S += Fsum(E2 * -540540, 471240).fmul(E5)

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

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

1034 S += Fsum(E2 * 417690, E22 * -255255, -875160).fmul(E2)

1035 S += 4084080

1036 if over:

1037 e *= over[0]

1038 return S.fdiv(e) # Fsum

1039

1040

1041 def _iterations(inst, i):

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

1043 '''

1044 if inst and i > 0:

1045 inst._iteration += i

1046

1047

1048 def _3over(a, b):

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

1050 '''

1051 return _over(_3_0, a * b)

1052

1053

1054 def _rC(unused, x, y):

1055 '''(INTERNAL) Defined only for C{y != 0} and C{x >= 0}.

1056 '''

1057 d = x - y

1058 if d < 0: # catch NaN

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

1060 d = -d

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

1062 elif d == 0: # XXX d < EPS0? or EPS02 or _EPSmin

1063 d, r = y, _1_0

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

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

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

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

1068 else: # PYCHOK no cover

1069 raise _ellipticError(Elliptic.fRC, x, y)

1070 return r / sqrt(d) # float

1071

1072

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

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

1075 '''

1076 L = _List(x, y, z)

1077 S = _Dsum()

1078 for a, m, r, s in L.amrs4(inst, True, _TolRF):

1079 if s:

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

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

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

1083 xy = x * y

1084 z = (x + y) / _3_0

1085 z2 = z**2

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

1087 xy - _6_0 * z2,

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

1089 (xy - z2) * _3_0 * z2,

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

1091

1092

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

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

1095 '''

1096 for a, b in _ab2(inst, x, y): # PYCHOK yield

1097 pass

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

1099

1100

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

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

1103 '''

1104 L = _List(x, y, z)

1105 for a, m, _, _ in L.amrs4(inst, False, _TolRF):

1106 pass

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

1108 z = neg(x + y)

1109 xy = x * y

1110 e2 = xy - z**2

1111 e3 = xy * z

1112 e4 = e2**2

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

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

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

1116 # converted to Horner-like form ...

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

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

1119 S += 240240

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

1121

1122

1123 def _rG2(inst, x, y, PI_=PI_4): # 2-args

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

1125 '''

1126 m = -1 # neg!

1127 S = None

1128 for a, b in _ab2(inst, x, y): # PYCHOK yield

1129 if S is None: # initial

1130 S = _Dsum()

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

1132 else:

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

1134 m *= 2

1135 return S(PI_).fover(a + b)

1136

1137

1138 def _rG3(inst, x, y, z): # 3-arg version

1139 '''(INTERNAL) C{x}, C{y} and C{z} all non-zero, see C{.fRG}.

1140 '''

1141 R = _RF3(inst, x, y, z) * z

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

1143 if rd: # Carlson, eq 1.7

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

1145 R += sqrt(x * y / z)

1146 return R.fover(_2_0)

1147

1148

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

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

1151 '''

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

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

1154

1155 L = _List(x, y, z, p)

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

1157 S = _Dsum()

1158 for a, m, _, s in L.amrs4(inst, True, _TolRD):

1159 if s:

1160 d = _xyzp(*s)

1161 if d:

1162 if n:

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

1164 n *= _1_64th # /= chokes PyChecker

1165 else:

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

1167 S += rc / (d * m)

1168 else: # PYCHOK no cover

1169 return NAN

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

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

1172 p2 = p**2

1173 p3 = p2 * p

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

1175 E2p = E2 * p

1176 xyz = x * y * z

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

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

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

1180 xyz * p2, *over) # Fsum

1181

1182# **) MIT License

1183#

1184# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1185#

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

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

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

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

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

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

1192#

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

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

1195#

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

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

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

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

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

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

1202# OTHER DEALINGS IN THE SOFTWARE.