Coverage for pygeodesy/fmath.py: 96%

275 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-02-23 11:26 -0500

1 

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

3 

4u'''Utilities using precision floating point summation. 

5''' 

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

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

8 

9from pygeodesy.basics import _copysign, copysign0, isint, len2 

10from pygeodesy.constants import EPS0, EPS02, EPS1, NAN, PI, PI_2, PI_4, \ 

11 _0_0, _0_125, _0_25, _0_5, _1_0, _N_1_0, \ 

12 _1_3rd, _1_5, _1_6th, _2_0, _2_3rd, _3_0, \ 

13 _isfinite, isnear0, isnear1, _over, remainder 

14from pygeodesy.errors import _IsnotError, LenError, _TypeError, _ValueError, \ 

15 _xError, _xkwds_get, _xkwds_pop2 

16from pygeodesy.fsums import _2float, _Powers, Fsum, _fsum, fsum, fsum1_, \ 

17 _pow_op_, Fmt, unstr 

18from pygeodesy.interns import MISSING, _few_, _h_, _invokation_, _negative_, \ 

19 _not_scalar_, _SPACE_, _too_ 

20from pygeodesy.lazily import _ALL_LAZY, _sys_version_info2 

21# from pygeodesy.streprs import Fmt, unstr # from .fsums 

22from pygeodesy.units import Int_, _isHeight, _isRadius, Float_ # PYCHOK for .heights 

23 

24from math import fabs, sqrt # pow 

25from operator import mul as _mul # in .triaxials 

26 

27__all__ = _ALL_LAZY.fmath 

28__version__ = '24.02.18' 

29 

30# sqrt(2) <https://WikiPedia.org/wiki/Square_root_of_2> 

31_0_4142 = 0.414213562373095 # sqrt(_2_0) - _1_0 

32 

33 

34class Fdot(Fsum): 

35 '''Precision dot product. 

36 ''' 

37 def __init__(self, a, *b, **name): 

38 '''New L{Fdot} precision dot product M{sum(a[i] * b[i] 

39 for i=0..len(a))}. 

40 

41 @arg a: Iterable, list, tuple, etc. (C{scalar}s). 

42 @arg b: Other values (C{scalar}s), all positional. 

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

44 

45 @raise OverflowError: Partial C{2sum} overflow. 

46 

47 @raise LenError: Unequal C{len(B{a})} and C{len(B{b})}. 

48 

49 @see: Function L{fdot} and method L{Fsum.fadd}. 

50 ''' 

51 Fsum.__init__(self, **name) 

52 self.fadd(_map_mul(a, b, Fdot)) 

53 

54 

55class Fhorner(Fsum): 

56 '''Precision polynomial evaluation using the Horner form. 

57 ''' 

58 def __init__(self, x, *cs, **name): 

59 '''New L{Fhorner} evaluation of the polynomial 

60 M{sum(cs[i] * x**i for i=0..len(cs))}. 

61 

62 @arg x: Polynomial argument (C{scalar}). 

63 @arg cs: Polynomial coeffients (C{scalar} or C{Fsum} 

64 instances), all positional. 

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

66 

67 @raise OverflowError: Partial C{2sum} overflow. 

68 

69 @raise TypeError: Non-scalar B{C{x}}. 

70 

71 @raise ValueError: Non-finite B{C{x}}. 

72 

73 @see: Function L{fhorner} and methods L{Fsum.fadd} and L{Fsum.fmul}. 

74 ''' 

75 Fsum.__init__(self, *cs[-1:], **name) 

76 if len(cs) > 1: 

77 x = _2float(x=x) 

78 _a = self._fadd # (other, op) 

79 _f = self._finite # (other, op) 

80 op = Fhorner.__name__ 

81 ps = self._ps 

82 for c in reversed(cs[:-1]): # multiply-accumulate 

83 ps[:] = [_f(p * x, op) for p in ps] 

84 _a(c, op) 

85 # assert self._ps is ps 

86 

87 

88class Fhypot(Fsum): 

89 '''Precision hypotenuse of summation. 

90 ''' 

91 def __init__(self, *xs, **power_name_RESIDUAL): 

92 '''New L{Fhypot} hypotenuse of (the I{power} of) several 

93 C{scalar} or C{Fsum} values. 

94 

95 @arg xs: One or more values to include (each C{scalar} 

96 or an C{Fsum} instance). 

97 @kwarg power_name_RESIDUAL: Optional exponent and root 

98 order C{B{power}=2}, C{B{name}=NN} and 

99 C{B{RESIDUAL}=None}, see L{Fsum.__init__}. 

100 ''' 

101 try: 

102 p, kwds = _xkwds_pop2(power_name_RESIDUAL, power=2) 

103 Fsum.__init__(self, **kwds) 

104 if xs: 

105 self._facc(_Powers(p, xs), up=False) # PYCHOK yield 

106 self._fset(self._fpow(_1_0 / p, _pow_op_), asis=True) 

107 except Exception as X: 

108 raise self._ErrorX(X, xs, power=p) 

109 

110 

111class Fpolynomial(Fsum): 

112 '''Precision polynomial evaluation. 

113 ''' 

114 def __init__(self, x, *cs, **name): 

115 '''New L{Fpolynomial} evaluation of the polynomial 

116 M{sum(cs[i] * x**i for i=0..len(cs))}. 

117 

118 @arg x: Polynomial argument (C{scalar}). 

119 @arg cs: Polynomial coeffients (C{scalar}s), all 

120 positional. 

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

122 

123 @raise OverflowError: Partial C{2sum} overflow. 

124 

125 @raise TypeError: Non-scalar B{C{x}}. 

126 

127 @raise ValueError: Non-finite B{C{x}}. 

128 

129 @see: Function L{fpolynomial} and method L{Fsum.fadd}. 

130 ''' 

131 Fsum.__init__(self, *cs[:1], **name) 

132 n = len(cs) - 1 

133 if n > 0: 

134 self.fadd(_map_mul(cs[1:], fpowers(x, n), Fpolynomial)) 

135 

136 

137class Fpowers(Fsum): 

138 '''Precision summation or powers, optimized for C{power=2}. 

139 ''' 

140 def __init__(self, power, *xs, **name_RESIDUAL): 

141 '''New L{Fpowers} sum of (the I{power} of) several C{scalar} 

142 or C{Fsum} values. 

143 

144 @arg power: The exponent (C{scalar} or C{Fsum}). 

145 @arg xs: One or more values to include (each C{scalar} 

146 or an C{Fsum} instance). 

147 @kwarg power_name_RESIDUAL: Optional exponent and root 

148 order C{B{power}=2}, C{B{name}=NN} and 

149 C{B{RESIDUAL}=None}, see L{Fsum.__init__}. 

150 ''' 

151 try: 

152 Fsum.__init__(self, **name_RESIDUAL) 

153 if xs: 

154 self._facc(_Powers(power, xs), up=False) # PYCHOK yield 

155 except Exception as X: 

156 raise self._ErrorX(X, xs, power=power) 

157 

158 

159class Fn_rt(Fsum): 

160 '''Precision n-th root of summation. 

161 ''' 

162 def __init__(self, root, *xs, **name_RESIDUAL): 

163 '''New L{Fn_rt} root of the precision sum of several 

164 C{scalar} or C{Fsum} values. 

165 

166 @arg root: The order (C{scalar} or C{Fsum}). 

167 @arg xs: Values to include (each C{scalar} or an 

168 C{Fsum} instance). 

169 @kwarg name_RESIDUAL: See L{Fsum.__init__}. 

170 ''' 

171 try: 

172 Fsum.__init__(self, *xs, **name_RESIDUAL) 

173 self._fset(self._fpow(_1_0 / root, _pow_op_), asis=True) 

174 except Exception as X: 

175 raise self._ErrorX(X, xs, root=root) 

176 

177 

178class Fcbrt(Fn_rt): 

179 '''Precision cubic root of summation. 

180 ''' 

181 def __init__(self, *xs, **name_RESIDUAL): 

182 '''New L{Fcbrt} cubic root of the precision sum of 

183 several C{scalar} or C{Fsum} values. 

184 

185 @arg xs: Values to include (each C{scalar} or an 

186 C{Fsum} instance). 

187 @kwarg name_RESIDUAL: See L{Fsum.__init__}. 

188 ''' 

189 Fn_rt.__init__(self, _3_0, *xs, **name_RESIDUAL) 

190 

191 

192class Fsqrt(Fn_rt): 

193 '''Precision square root of summation. 

194 ''' 

195 def __init__(self, *xs, **name_RESIDUAL): 

196 '''New L{Fsqrt} square root of the precision sum of 

197 several C{scalar} or C{Fsum} values. 

198 

199 @arg xs: Values to include (each C{scalar} or an 

200 C{Fsum} instance). 

201 @kwarg name_RESIDUAL: See L{Fsum.__init__}. 

202 ''' 

203 Fn_rt.__init__(self, _2_0, *xs, **name_RESIDUAL) 

204 

205 

206def bqrt(x): 

207 '''Return the 4-th, I{bi-quadratic} or I{quartic} root, M{x**(1 / 4)}. 

208 

209 @arg x: Value (C{scalar}). 

210 

211 @return: I{Quartic} root (C{float}). 

212 

213 @raise ValueError: Negative B{C{x}}. 

214 

215 @see: Functions L{zcrt} and L{zqrt}. 

216 ''' 

217 return _root(x, _0_25, bqrt) 

218 

219 

220try: 

221 from math import cbrt # Python 3.11+ 

222 

223 def cbrt2(x): 

224 '''Compute the cube root I{squared} M{x**(2/3)}. 

225 ''' 

226 return cbrt(x)**2 # cbrt(-0.0*2) == -0.0 

227 

228except ImportError: # Python 3.10- 

229 

230 def cbrt(x): 

231 '''Compute the cube root M{x**(1/3)}. 

232 

233 @arg x: Value (C{scalar}). 

234 

235 @return: Cubic root (C{float}). 

236 

237 @see: Functions L{cbrt2} and L{sqrt3}. 

238 ''' 

239 # <https://archive.lib.MSU.edu/crcmath/math/math/r/r021.htm> 

240 # simpler and more accurate than Ken Turkowski's CubeRoot, see 

241 # <https://People.FreeBSD.org/~lstewart/references/apple_tr_kt32_cuberoot.pdf> 

242 return _copysign(pow(fabs(x), _1_3rd), x) # cbrt(-0.0) == -0.0 

243 

244 def cbrt2(x): # PYCHOK attr 

245 '''Compute the cube root I{squared} M{x**(2/3)}. 

246 

247 @arg x: Value (C{scalar}). 

248 

249 @return: Cube root I{squared} (C{float}). 

250 

251 @see: Functions L{cbrt} and L{sqrt3}. 

252 ''' 

253 return pow(fabs(x), _2_3rd) # XXX pow(fabs(x), _1_3rd)**2 

254 

255 

256def euclid(x, y): 

257 '''I{Appoximate} the norm M{sqrt(x**2 + y**2)} by 

258 M{max(abs(x), abs(y)) + min(abs(x), abs(y)) * 0.4142...}. 

259 

260 @arg x: X component (C{scalar}). 

261 @arg y: Y component (C{scalar}). 

262 

263 @return: Appoximate norm (C{float}). 

264 

265 @see: Function L{euclid_}. 

266 ''' 

267 x, y = fabs(x), fabs(y) 

268 if x < y: 

269 x, y = y, x 

270 return x + y * _0_4142 # XXX * _0_5 before 20.10.02 

271 

272 

273def euclid_(*xs): 

274 '''I{Appoximate} the norm M{sqrt(sum(x**2 for x in xs))} 

275 by cascaded L{euclid}. 

276 

277 @arg xs: X arguments, positional (C{scalar}s). 

278 

279 @return: Appoximate norm (C{float}). 

280 

281 @see: Function L{euclid}. 

282 ''' 

283 e = _0_0 

284 for x in sorted(map(fabs, xs)): # XXX not reverse=True 

285 # e = euclid(x, e) 

286 if e < x: 

287 e, x = x, e 

288 if x: 

289 e += x * _0_4142 

290 return e 

291 

292 

293def facos1(x): 

294 '''Fast approximation of L{pygeodesy.acos1}C{(B{x})}. 

295 

296 @see: U{ShaderFastLibs.h<https://GitHub.com/michaldrobot/ 

297 ShaderFastLibs/blob/master/ShaderFastMathLib.h>}. 

298 ''' 

299 a = fabs(x) 

300 if a < EPS0: 

301 r = PI_2 

302 elif a < EPS1: 

303 H = Fhorner(-a, 1.5707288, 0.2121144, 0.0742610, 0.0187293) 

304 r = float(H * sqrt(_1_0 - a)) 

305 if x < 0: 

306 r = PI - r 

307 else: 

308 r = PI if x < 0 else _0_0 

309 return r 

310 

311 

312def fasin1(x): # PYCHOK no cover 

313 '''Fast approximation of L{pygeodesy.asin1}C{(B{x})}. 

314 

315 @see: L{facos1}. 

316 ''' 

317 return PI_2 - facos1(x) 

318 

319 

320def fatan(x): 

321 '''Fast approximation of C{atan(B{x})}. 

322 ''' 

323 a = fabs(x) 

324 if a < _1_0: 

325 r = fatan1(a) if a else _0_0 

326 elif a > _1_0: 

327 r = PI_2 - fatan1(_1_0 / a) # == fatan2(a, _1_0) 

328 else: 

329 r = PI_4 

330 if x < 0: # copysign0(r, x) 

331 r = -r 

332 return r 

333 

334 

335def fatan1(x): 

336 '''Fast approximation of C{atan(B{x})} for C{0 <= B{x} <= 1}, I{unchecked}. 

337 

338 @see: U{ShaderFastLibs.h<https://GitHub.com/michaldrobot/ShaderFastLibs/ 

339 blob/master/ShaderFastMathLib.h>} and U{Efficient approximations 

340 for the arctangent function<http://www-Labs.IRO.UMontreal.CA/ 

341 ~mignotte/IFT2425/Documents/EfficientApproximationArctgFunction.pdf>}, 

342 IEEE Signal Processing Magazine, 111, May 2006. 

343 ''' 

344 # Eq (9): PI_4 * x - x * (abs(x) - 1) * (0.2447 + 0.0663 * abs(x)), for -1 < x < 1 

345 # PI_4 * x - (x**2 - x) * (0.2447 + 0.0663 * x), for 0 < x - 1 

346 # x * (1.0300981633974482 + x * (-0.1784 - x * 0.0663)) 

347 H = Fhorner(x, _0_0, 1.0300982, -0.1784, -0.0663) 

348 return float(H) 

349 

350 

351def fatan2(y, x): 

352 '''Fast approximation of C{atan2(B{y}, B{x})}. 

353 

354 @see: U{fastApproximateAtan(x, y)<https://GitHub.com/CesiumGS/cesium/blob/ 

355 master/Source/Shaders/Builtin/Functions/fastApproximateAtan.glsl>} 

356 and L{fatan1}. 

357 ''' 

358 b, a = fabs(y), fabs(x) 

359 if a < b: 

360 r = (PI_2 - fatan1(a / b)) if a else PI_2 

361 elif b < a: 

362 r = fatan1(b / a) if b else _0_0 

363 elif a: # == b != 0 

364 r = PI_4 

365 else: # a == b == 0 

366 return _0_0 

367 if x < 0: 

368 r = PI - r 

369 if y < 0: # copysign0(r, y) 

370 r = -r 

371 return r 

372 

373 

374def favg(v1, v2, f=_0_5): 

375 '''Return the average of two values. 

376 

377 @arg v1: One value (C{scalar}). 

378 @arg v2: Other value (C{scalar}). 

379 @kwarg f: Optional fraction (C{float}). 

380 

381 @return: M{v1 + f * (v2 - v1)} (C{float}). 

382 ''' 

383# @raise ValueError: Fraction out of range. 

384# ''' 

385# if not 0 <= f <= 1: # XXX restrict fraction? 

386# raise _ValueError(fraction=f) 

387 # v1 + f * (v2 - v1) == v1 * (1 - f) + v2 * f 

388 return fsum1_(v1, -f * v1, f * v2) 

389 

390 

391def fdot(a, *b): 

392 '''Return the precision dot product M{sum(a[i] * b[i] for 

393 i=0..len(a))}. 

394 

395 @arg a: Iterable, list, tuple, etc. (C{scalar}s). 

396 @arg b: All positional arguments (C{scalar}s). 

397 

398 @return: Dot product (C{float}). 

399 

400 @raise LenError: Unequal C{len(B{a})} and C{len(B{b})}. 

401 

402 @see: Class L{Fdot} and U{Algorithm 5.10 B{DotK} 

403 <https://www.TUHH.De/ti3/paper/rump/OgRuOi05.pdf>}. 

404 ''' 

405 return fsum(_map_mul(a, b, fdot)) 

406 

407 

408def fdot3(a, b, c, start=0): 

409 '''Return the precision dot product M{start + 

410 sum(a[i] * b[i] * c[i] for i=0..len(a))}. 

411 

412 @arg a: Iterable, list, tuple, etc. (C{scalar}s). 

413 @arg b: Iterable, list, tuple, etc. (C{scalar}s). 

414 @arg c: Iterable, list, tuple, etc. (C{scalar}s). 

415 @kwarg start: Optional bias (C{scalar}). 

416 

417 @return: Dot product (C{float}). 

418 

419 @raise LenError: Unequal C{len(B{a})}, C{len(B{b})} 

420 and/or C{len(B{c})}. 

421 

422 @raise OverflowError: Partial C{2sum} overflow. 

423 ''' 

424 def _mul3(a, b, c): # map function 

425 return a * b * c 

426 

427 def _muly(a, b, c, start): 

428 yield start 

429 for abc in map(_mul3, a, b, c): 

430 yield abc 

431 

432 if not len(a) == len(b) == len(c): 

433 raise LenError(fdot3, a=len(a), b=len(b), c=len(c)) 

434 

435 return fsum(_muly(a, b, c, start) if start else map(_mul3, a, b, c)) 

436 

437 

438def fhorner(x, *cs): 

439 '''Evaluate the polynomial M{sum(cs[i] * x**i for 

440 i=0..len(cs))} using the Horner form. 

441 

442 @arg x: Polynomial argument (C{scalar}). 

443 @arg cs: Polynomial coeffients (C{scalar}s). 

444 

445 @return: Horner value (C{float}). 

446 

447 @raise OverflowError: Partial C{2sum} overflow. 

448 

449 @raise TypeError: Non-scalar B{C{x}}. 

450 

451 @raise ValueError: No B{C{cs}} coefficients or B{C{x}} is not finite. 

452 

453 @see: Function L{fpolynomial} and class L{Fhorner}. 

454 ''' 

455 H = Fhorner(x, *cs) 

456 return float(H) 

457 

458 

459def fidw(xs, ds, beta=2): 

460 '''Interpolate using U{Inverse Distance Weighting 

461 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW). 

462 

463 @arg xs: Known values (C{scalar}s). 

464 @arg ds: Non-negative distances (C{scalar}s). 

465 @kwarg beta: Inverse distance power (C{int}, 0, 1, 2, or 3). 

466 

467 @return: Interpolated value C{x} (C{float}). 

468 

469 @raise LenError: Unequal or zero C{len(B{ds})} and C{len(B{xs})}. 

470 

471 @raise ValueError: Invalid B{C{beta}}, negative B{C{ds}} value, 

472 weighted B{C{ds}} below L{EPS}. 

473 

474 @note: Using C{B{beta}=0} returns the mean of B{C{xs}}. 

475 ''' 

476 n, xs = len2(xs) 

477 d, ds = len2(ds) 

478 if n != d or n < 1: 

479 raise LenError(fidw, xs=n, ds=d) 

480 

481 d, x = min(zip(ds, xs)) 

482 if d > EPS0 and n > 1: 

483 b = -Int_(beta=beta, low=0, high=3) 

484 if b < 0: 

485 ws = tuple(float(d)**b for d in ds) 

486 t = fsum(_map_mul1(xs, ws)) # fdot(xs, *ws) 

487 x = _over(t, fsum(ws, floats=True)) 

488 else: # b == 0 

489 x = fsum(xs) / n # fmean(xs) 

490 elif d < 0: # PYCHOK no cover 

491 n = Fmt.SQUARE(ds=ds.index(d)) 

492 raise _ValueError(n, d, txt=_negative_) 

493 return x 

494 

495 

496def fmean(xs): 

497 '''Compute the accurate mean M{sum(xs[i] for 

498 i=0..len(xs)) / len(xs)}. 

499 

500 @arg xs: Values (C{scalar} or L{Fsum} instances). 

501 

502 @return: Mean value (C{float}). 

503 

504 @raise OverflowError: Partial C{2sum} overflow. 

505 

506 @raise ValueError: No B{C{xs}} values. 

507 ''' 

508 n, xs = len2(xs) 

509 if n > 0: 

510 return fsum(xs) / n # if n > 1 else _2float(index=0, xs=xs[0]) 

511 raise _ValueError(xs=xs) 

512 

513 

514def fmean_(*xs): 

515 '''Compute the accurate mean M{sum(xs[i] for 

516 i=0..len(xs)) / len(xs)}. 

517 

518 @arg xs: Values (C{scalar} or L{Fsum} instances). 

519 

520 @return: Mean value (C{float}). 

521 

522 @raise OverflowError: Partial C{2sum} overflow. 

523 

524 @raise ValueError: No B{C{xs}} values. 

525 ''' 

526 return fmean(xs) 

527 

528 

529def fpolynomial(x, *cs, **over): 

530 '''Evaluate the polynomial M{sum(cs[i] * x**i for 

531 i=0..len(cs)) [/ over]}. 

532 

533 @arg x: Polynomial argument (C{scalar}). 

534 @arg cs: Polynomial coeffients (C{scalar}s), all 

535 positional. 

536 @kwarg over: Optional, final divisor (C{scalar} 

537 

538 @return: Polynomial value (C{float}). 

539 

540 @raise OverflowError: Partial C{2sum} overflow. 

541 

542 @raise TypeError: Non-scalar B{C{x}}. 

543 

544 @raise ValueError: No B{C{cs}} coefficients or B{C{x}} is not finite. 

545 

546 @see: Function L{fhorner} and class L{Fpolynomial}. 

547 ''' 

548 P = Fpolynomial(x, *cs) 

549 d = _xkwds_get(over, over=0) if over else 0 

550 return P.fover(d) if d else float(P) 

551 

552 

553def fpowers(x, n, alts=0): 

554 '''Return a series of powers M{[x**i for i=1..n]}. 

555 

556 @arg x: Value (C{scalar}). 

557 @arg n: Highest exponent (C{int}). 

558 @kwarg alts: Only alternating powers, starting with 

559 this exponent (C{int}). 

560 

561 @return: Powers of B{C{x}} (C{float}s or C{int}s). 

562 

563 @raise TypeError: Non-scalar B{C{x}} or B{C{n}} not C{int}. 

564 

565 @raise ValueError: Non-finite B{C{x}} or non-positive B{C{n}}. 

566 ''' 

567 if not isint(n): 

568 raise _IsnotError(int.__name__, n=n) 

569 elif n < 1: 

570 raise _ValueError(n=n) 

571 

572 p = t = x if isint(x) else _2float(x=x) 

573 ps = [p] 

574 _a = ps.append 

575 for _ in range(1, n): 

576 p *= t 

577 _a(p) 

578 

579 if alts > 0: # x**2, x**4, ... 

580 # ps[alts-1::2] chokes PyChecker 

581 ps = ps[slice(alts-1, None, 2)] 

582 

583 return ps 

584 

585 

586try: 

587 from math import prod as fprod # Python 3.8 

588except ImportError: 

589 

590 def fprod(xs, start=1): 

591 '''Iterable product, like C{math.prod} or C{numpy.prod}. 

592 

593 @arg xs: Terms to be multiplied, an iterable, list, 

594 tuple, etc. (C{scalar}s). 

595 @kwarg start: Initial term, also the value returned 

596 for an empty B{C{xs}} (C{scalar}). 

597 

598 @return: The product (C{float}). 

599 

600 @see: U{NumPy.prod<https://docs.SciPy.org/doc/ 

601 numpy/reference/generated/numpy.prod.html>}. 

602 ''' 

603 return freduce(_mul, xs, start) 

604 

605 

606def frange(start, number, step=1): 

607 '''Generate a range of C{float}s. 

608 

609 @arg start: First value (C{float}). 

610 @arg number: The number of C{float}s to generate (C{int}). 

611 @kwarg step: Increment value (C{float}). 

612 

613 @return: A generator (C{float}s). 

614 

615 @see: U{NumPy.prod<https://docs.SciPy.org/doc/ 

616 numpy/reference/generated/numpy.arange.html>}. 

617 ''' 

618 if not isint(number): 

619 raise _IsnotError(int.__name__, number=number) 

620 for i in range(number): 

621 yield start + i * step 

622 

623 

624try: 

625 from functools import reduce as freduce 

626except ImportError: 

627 try: 

628 freduce = reduce # PYCHOK expected 

629 except NameError: # Python 3+ 

630 

631 def freduce(f, xs, *start): 

632 '''For missing C{functools.reduce}. 

633 ''' 

634 if start: 

635 r = v = start[0] 

636 else: 

637 r, v = 0, MISSING 

638 for v in xs: 

639 r = f(r, v) 

640 if v is MISSING: 

641 raise _TypeError(xs=(), start=MISSING) 

642 return r 

643 

644 

645def fremainder(x, y): 

646 '''Remainder in range C{[-B{y / 2}, B{y / 2}]}. 

647 

648 @arg x: Numerator (C{scalar}). 

649 @arg y: Modulus, denominator (C{scalar}). 

650 

651 @return: Remainder (C{scalar}, preserving signed 

652 0.0) or C{NAN} for any non-finite B{C{x}}. 

653 

654 @raise ValueError: Infinite or near-zero B{C{y}}. 

655 

656 @see: I{Karney}'s U{Math.remainder<https://PyPI.org/ 

657 project/geographiclib/>} and Python 3.7+ 

658 U{math.remainder<https://docs.Python.org/3/ 

659 library/math.html#math.remainder>}. 

660 ''' 

661 # with Python 2.7.16 and 3.7.3 on macOS 10.13.6 and 

662 # with Python 3.10.2 on macOS 12.2.1 M1 arm64 native 

663 # fmod( 0, 360) == 0.0 

664 # fmod( 360, 360) == 0.0 

665 # fmod(-0, 360) == 0.0 

666 # fmod(-0.0, 360) == -0.0 

667 # fmod(-360, 360) == -0.0 

668 # however, using the % operator ... 

669 # 0 % 360 == 0 

670 # 360 % 360 == 0 

671 # 360.0 % 360 == 0.0 

672 # -0 % 360 == 0 

673 # -360 % 360 == 0 == (-360) % 360 

674 # -0.0 % 360 == 0.0 == (-0.0) % 360 

675 # -360.0 % 360 == 0.0 == (-360.0) % 360 

676 

677 # On Windows 32-bit with python 2.7, math.fmod(-0.0, 360) 

678 # == +0.0. This fixes this bug. See also Math::AngNormalize 

679 # in the C++ library, Math.sincosd has a similar fix. 

680 if _isfinite(x): 

681 try: 

682 r = remainder(x, y) if x else x 

683 except Exception as e: 

684 raise _xError(e, unstr(fremainder, x, y)) 

685 else: # handle x INF and NINF as NAN 

686 r = NAN 

687 return r 

688 

689 

690if _sys_version_info2 < (3, 8): # PYCHOK no cover 

691 from math import hypot # OK in Python 3.7- 

692 

693 def hypot_(*xs): 

694 '''Compute the norm M{sqrt(sum(x**2 for x in xs))}. 

695 

696 Similar to Python 3.8+ n-dimension U{math.hypot 

697 <https://docs.Python.org/3.8/library/math.html#math.hypot>}, 

698 but exceptions, C{nan} and C{infinite} values are 

699 handled differently. 

700 

701 @arg xs: X arguments (C{scalar}s), all positional. 

702 

703 @return: Norm (C{float}). 

704 

705 @raise OverflowError: Partial C{2sum} overflow. 

706 

707 @raise ValueError: Invalid or no B{C{xs}} values. 

708 

709 @note: The Python 3.8+ Euclidian distance U{math.dist 

710 <https://docs.Python.org/3.8/library/math.html#math.dist>} 

711 between 2 I{n}-dimensional points I{p1} and I{p2} can be 

712 computed as M{hypot_(*((c1 - c2) for c1, c2 in zip(p1, p2)))}, 

713 provided I{p1} and I{p2} have the same, non-zero length I{n}. 

714 ''' 

715 h, x2 = _h_x2(xs) 

716 return (h * sqrt(x2)) if x2 else _0_0 

717 

718elif _sys_version_info2 < (3, 10): 

719 # In Python 3.8 and 3.9 C{math.hypot} is inaccurate, see 

720 # U{agdhruv<https://GitHub.com/geopy/geopy/issues/466>}, 

721 # U{cffk<https://Bugs.Python.org/issue43088>} and module 

722 # U{geomath.py<https://PyPI.org/project/geographiclib/1.52>} 

723 

724 def hypot(x, y): 

725 '''Compute the norm M{sqrt(x**2 + y**2)}. 

726 

727 @arg x: X argument (C{scalar}). 

728 @arg y: Y argument (C{scalar}). 

729 

730 @return: C{sqrt(B{x}**2 + B{y}**2)} (C{float}). 

731 ''' 

732 if x: 

733 h = sqrt(x**2 + y**2) if y else fabs(x) 

734 elif y: 

735 h = fabs(y) 

736 else: 

737 h = _0_0 

738 return h 

739 

740 from math import hypot as hypot_ # PYCHOK in Python 3.8 and 3.9 

741else: 

742 from math import hypot # PYCHOK in Python 3.10+ 

743 hypot_ = hypot 

744 

745 

746def _h_x2(xs): 

747 '''(INTERNAL) Helper for L{hypot_} and L{hypot2_}. 

748 ''' 

749 if xs: 

750 n, xs = len2(xs) 

751 if n > 0: 

752 h = float(max(map(fabs, xs))) 

753 if h < EPS0: 

754 x2 = _0_0 

755 else: # math.fsum, see C{_hypot21_} below 

756 x2 = _fsum(_x2_h2(_1_0, xs, h, _N_1_0)) 

757 return h, x2 

758 

759 raise _ValueError(xs=xs, txt=_too_(_few_)) 

760 

761 

762def hypot1(x): 

763 '''Compute the norm M{sqrt(1 + x**2)}. 

764 

765 @arg x: Argument (C{scalar}). 

766 

767 @return: Norm (C{float}). 

768 ''' 

769 return hypot(_1_0, x) if x else _1_0 

770 

771 

772def hypot2(x, y): 

773 '''Compute the I{squared} norm M{x**2 + y**2}. 

774 

775 @arg x: X argument (C{scalar}). 

776 @arg y: Y argument (C{scalar}). 

777 

778 @return: C{B{x}**2 + B{y}**2} (C{float}). 

779 ''' 

780 if x: 

781 if y: 

782 if fabs(x) < fabs(y): 

783 x, y = y, x 

784 h2 = x**2 * ((y / x)**2 + _1_0) 

785 else: 

786 h2 = x**2 

787 elif y: 

788 h2 = y**2 

789 else: 

790 h2 = _0_0 

791 return h2 

792 

793 

794def hypot2_(*xs): 

795 '''Compute the I{squared} norm C{sum(x**2 for x in B{xs})}. 

796 

797 @arg xs: X arguments (C{scalar}s), all positional. 

798 

799 @return: Squared norm (C{float}). 

800 

801 @raise OverflowError: Partial C{2sum} overflow. 

802 

803 @raise ValueError: Invalid or no B{C{xs}} value. 

804 

805 @see: Function L{hypot_}. 

806 ''' 

807 h, x2 = _h_x2(xs) 

808 return (h**2 * x2) if x2 else _0_0 

809 

810 

811def _map_mul(a, b, where): 

812 '''(INTERNAL) Yield each B{C{a * b}}. 

813 ''' 

814 n = len(b) 

815 if len(a) != n: # PYCHOK no cover 

816 raise LenError(where, a=len(a), b=n) 

817 return map(_mul, a, b) if n > 3 else _map_mul1(a, b) 

818 

819 

820def _map_mul1(a, b): 

821 '''(INTERNAL) Yield each B{C{a * b}}, 1-primed. 

822 ''' 

823 yield _1_0 

824 for ab in map(_mul, a, b): 

825 if ab: 

826 yield ab 

827 yield _N_1_0 

828 

829 

830def norm2(x, y): 

831 '''Normalize a 2-dimensional vector. 

832 

833 @arg x: X component (C{scalar}). 

834 @arg y: Y component (C{scalar}). 

835 

836 @return: 2-Tuple C{(x, y)}, normalized. 

837 

838 @raise ValueError: Invalid B{C{x}} or B{C{y}} 

839 or zero norm. 

840 ''' 

841 h = hypot(x, y) 

842 if not h: 

843 x = y = _0_0 # pass? 

844 elif not isnear1(h): 

845 try: 

846 x, y = x / h, y / h 

847 except Exception as e: 

848 raise _xError(e, x=x, y=y, h=h) 

849 return x, y 

850 

851 

852def norm_(*xs): 

853 '''Normalize all n-dimensional vector components. 

854 

855 @arg xs: Components (C{scalar}s), all positional. 

856 

857 @return: Yield each component, normalized. 

858 

859 @raise ValueError: Invalid or insufficent B{C{xs}} 

860 or zero norm. 

861 ''' 

862 h = hypot_(*xs) 

863 if h: 

864 try: 

865 for i, x in enumerate(xs): 

866 yield x / h 

867 except Exception as e: 

868 raise _xError(e, Fmt.SQUARE(xs=i), x, _h_, h) 

869 else: 

870 for _ in xs: 

871 yield _0_0 

872 

873 

874def _root(x, p, where): 

875 '''(INTERNAL) Raise C{x} to power C{0 < p < 1}. 

876 ''' 

877 if x < 0: 

878 t = _SPACE_(_invokation_, where.__name__) 

879 raise _ValueError(unstr(t, x), txt=_negative_) 

880 return pow(x, p) if x else _0_0 

881 

882 

883def sqrt0(x, Error=None): 

884 '''Return the square root iff C{B{x} >} L{EPS02}. 

885 

886 @arg x: Value (C{scalar}). 

887 @kwarg Error: Error to raise for negative B{C{x}}. 

888 

889 @return: Square root (C{float}) or C{0.0}. 

890 

891 @note: Any C{B{x} < }L{EPS02} I{including} C{B{x} < 0} 

892 returns C{0.0}. 

893 ''' 

894 if Error and x < 0: 

895 raise Error(Fmt.PAREN(sqrt=x)) 

896 return sqrt(x) if x > EPS02 else (_0_0 if x < EPS02 else EPS0) 

897 

898 

899def sqrt3(x): 

900 '''Return the square root, I{cubed} M{sqrt(x)**3} or M{sqrt(x**3)}. 

901 

902 @arg x: Value (C{scalar}). 

903 

904 @return: Square root I{cubed} (C{float}). 

905 

906 @raise ValueError: Negative B{C{x}}. 

907 

908 @see: Functions L{cbrt} and L{cbrt2}. 

909 ''' 

910 return _root(x, _1_5, sqrt3) 

911 

912 

913def sqrt_a(h, b): 

914 '''Compute C{I{a}} side of a right-angled triangle from 

915 C{sqrt(B{h}**2 - B{b}**2)}. 

916 

917 @arg h: Hypotenuse or outer annulus radius (C{scalar}). 

918 @arg b: Triangle side or inner annulus radius (C{scalar}). 

919 

920 @return: C{copysign(I{a}, B{h})} or C{unsigned 0.0} (C{float}). 

921 

922 @raise TypeError: Non-scalar B{C{h}} or B{C{b}}. 

923 

924 @raise ValueError: If C{abs(B{h}) < abs(B{b})}. 

925 

926 @see: Inner tangent chord B{I{d}} of an U{annulus 

927 <https://WikiPedia.org/wiki/Annulus_(mathematics)>} 

928 and function U{annulus_area<https://People.SC.FSU.edu/ 

929 ~jburkardt/py_src/geometry/geometry.py>}. 

930 ''' 

931 try: 

932 if not (_isHeight(h) and _isRadius(b)): 

933 raise TypeError(_not_scalar_) 

934 elif isnear0(h): # PYCHOK no cover 

935 c, b = fabs(h), fabs(b) 

936 d = c - b 

937 if d < 0: 

938 raise ValueError('abs(h) < abs(b)') 

939 a = copysign0(sqrt((c + b) * d), h) if d > 0 else _0_0 

940 else: 

941 c = float(h) 

942 s = _1_0 - (b / c)**2 

943 if s < 0: 

944 raise ValueError('abs(h) < abs(b)') 

945 a = (sqrt(s) * c) if 0 < s < 1 else (c if s else _0_0) 

946 except Exception as x: 

947 raise _xError(x, h=h, b=b) 

948 return a 

949 

950 

951def _x2_h2(s, xs, h, e): 

952 '''(INTERNAL) Yield M{(x / h)**2 for x in xs}. 

953 ''' 

954 yield s 

955 if h in (_0_0, _1_0): 

956 for x in xs: 

957 if x: 

958 yield x**2 

959 else: 

960 for x in xs: 

961 if x: 

962 yield (x / h)**2 

963 yield e 

964 

965 

966def zcrt(x): 

967 '''Return the 6-th, I{zenzi-cubic} root, M{x**(1 / 6)}. 

968 

969 @arg x: Value (C{scalar}). 

970 

971 @return: I{Zenzi-cubic} root (C{float}). 

972 

973 @see: Functions L{bqrt} and L{zqrt}. 

974 

975 @raise ValueError: Negative B{C{x}}. 

976 ''' 

977 return _root(x, _1_6th, zcrt) 

978 

979 

980def zqrt(x): 

981 '''Return the 8-th, I{zenzi-quartic} or I{squared-quartic} root, M{x**(1 / 8)}. 

982 

983 @arg x: Value (C{scalar}). 

984 

985 @return: I{Zenzi-quartic} root (C{float}). 

986 

987 @see: Functions L{bqrt} and L{zcrt}. 

988 

989 @raise ValueError: Negative B{C{x}}. 

990 ''' 

991 return _root(x, _0_125, zqrt) 

992 

993# **) MIT License 

994# 

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

996# 

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

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

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

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

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

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

1003# 

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

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

1006# 

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

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

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

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

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

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

1013# OTHER DEALINGS IN THE SOFTWARE.