Coverage for pygeodesy/fmath.py: 93%

265 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-06-07 08:37 -0400

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, isscalar, len2 

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

11 _0_0, _0_5, _1_0, _N_1_0, _1_3rd, _1_5, \ 

12 _2_0, _2_3rd, _3_0, isnear0, isnear1, \ 

13 _isfinite, remainder as _remainder 

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

15 _xError, _xkwds_get, _xkwds_pop 

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

17 _pow_op_, Fmt, unstr 

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

19 _singular_, _too_ 

20from pygeodesy.lazily import _ALL_LAZY, _sys_version_info2 

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

22from pygeodesy.units import Int_, 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__ = '23.05.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_a_x_b(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 = _xkwds_pop(power_name_RESIDUAL, power=2) 

103 Fsum.__init__(self, **power_name_RESIDUAL) 

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_a_x_b(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 

206try: 

207 from math import cbrt # Python 3.11+ 

208 

209 def cbrt2(x): 

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

211 ''' 

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

213 

214except ImportError: # Python 3.10- 

215 

216 def cbrt(x): 

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

218 

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

220 

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

222 

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

224 ''' 

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

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

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

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

229 

230 def cbrt2(x): # PYCHOK attr 

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

232 

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

234 

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

236 

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

238 ''' 

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

240 

241 

242def euclid(x, y): 

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

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

245 

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

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

248 

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

250 

251 @see: Function L{euclid_}. 

252 ''' 

253 x, y = fabs(x), fabs(y) 

254 if y > x: 

255 x, y = y, x 

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

257 

258 

259def euclid_(*xs): 

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

261 by cascaded L{euclid}. 

262 

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

264 

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

266 

267 @see: Function L{euclid}. 

268 ''' 

269 e = _0_0 

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

271 # e = euclid(x, e) 

272 if x > e: 

273 e, x = x, e 

274 if x: 

275 e += x * _0_4142 

276 return e 

277 

278 

279def facos1(x): 

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

281 

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

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

284 ''' 

285 a = fabs(x) 

286 if a < EPS0: 

287 r = PI_2 

288 elif a < EPS1: 

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

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

291 if x < 0: 

292 r = PI - r 

293 else: 

294 r = PI if x < 0 else _0_0 

295 return r 

296 

297 

298def fasin1(x): # PYCHOK no cover 

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

300 

301 @see: L{facos1}. 

302 ''' 

303 return PI_2 - facos1(x) 

304 

305 

306def fatan(x): 

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

308 ''' 

309 a = fabs(x) 

310 if a < _1_0: 

311 r = fatan1(a) if a else _0_0 

312 elif a > _1_0: 

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

314 else: 

315 r = PI_4 

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

317 r = -r 

318 return r 

319 

320 

321def fatan1(x): 

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

323 

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

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

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

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

328 IEEE Signal Processing Magazine, 111, May 2006. 

329 ''' 

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

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

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

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

334 return float(H) 

335 

336 

337def fatan2(y, x): 

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

339 

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

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

342 and L{fatan1}. 

343 ''' 

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

345 if a < b: 

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

347 elif b < a: 

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

349 elif a: # == b != 0 

350 r = PI_4 

351 else: # a == b == 0 

352 return _0_0 

353 if x < 0: 

354 r = PI - r 

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

356 r = -r 

357 return r 

358 

359 

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

361 '''Return the average of two values. 

362 

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

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

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

366 

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

368 ''' 

369# @raise ValueError: Fraction out of range. 

370# ''' 

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

372# raise _ValueError(fraction=f) 

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

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

375 

376 

377def fdot(a, *b): 

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

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

380 

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

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

383 

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

385 

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

387 

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

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

390 ''' 

391 return fsum(_map_a_x_b(a, b, fdot)) 

392 

393 

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

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

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

397 

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

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

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

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

402 

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

404 

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

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

407 

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

409 ''' 

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

411 return a * b * c 

412 

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

414 yield start 

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

416 yield abc 

417 

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

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

420 

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

422 

423 

424def fhorner(x, *cs): 

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

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

427 

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

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

430 

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

432 

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

434 

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

436 

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

438 

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

440 ''' 

441 H = Fhorner(x, *cs) 

442 return float(H) 

443 

444 

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

446 '''Interpolate using U{Inverse Distance Weighting 

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

448 

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

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

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

452 

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

454 

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

456 

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

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

459 

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

461 ''' 

462 n, xs = len2(xs) 

463 d, ds = len2(ds) 

464 if n != d or n < 1: 

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

466 

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

468 if d > EPS0 and n > 1: 

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

470 if b < 0: 

471 ds = tuple(d**b for d in ds) 

472 d = fsum(ds, floats=True) 

473 if isnear0(d): # PYCHOK no cover 

474 n = Fmt.PAREN(fsum='ds') 

475 raise _ValueError(n, d, txt=_singular_) 

476 x = Fdot(xs, *ds).fover(d) 

477 else: # b == 0 

478 x = fsum(xs, floats=True) / n # fmean(xs) 

479 elif d < 0: # PYCHOK no cover 

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

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

482 return x 

483 

484 

485def fmean(xs): 

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

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

488 

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

490 

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

492 

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

494 

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

496 ''' 

497 n, xs = len2(xs) 

498 if n > 0: 

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

500 raise _ValueError(xs=xs) 

501 

502 

503def fmean_(*xs): 

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

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

506 

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

508 

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

510 

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

512 

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

514 ''' 

515 return fmean(xs) 

516 

517 

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

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

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

521 

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

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

524 positional. 

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

526 

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

528 

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

530 

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

532 

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

534 

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

536 ''' 

537 P = Fpolynomial(x, *cs) 

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

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

540 

541 

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

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

544 

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

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

547 @kwarg alts: Only alternating powers, starting with 

548 this exponent (C{int}). 

549 

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

551 

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

553 

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

555 ''' 

556 if not isint(n): 

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

558 elif n < 1: 

559 raise _ValueError(n=n) 

560 

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

562 ps = [p] 

563 a_ = ps.append 

564 for _ in range(1, n): 

565 p *= t 

566 a_(p) 

567 

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

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

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

571 

572 return ps 

573 

574 

575try: 

576 from math import prod as fprod # Python 3.8 

577except ImportError: 

578 

579 def fprod(xs, start=1): 

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

581 

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

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

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

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

586 

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

588 

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

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

591 ''' 

592 return freduce(_mul, xs, start) 

593 

594 

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

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

597 

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

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

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

601 

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

603 

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

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

606 ''' 

607 if not isint(number): 

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

609 for i in range(number): 

610 yield start + i * step 

611 

612 

613try: 

614 from functools import reduce as freduce 

615except ImportError: 

616 try: 

617 freduce = reduce # PYCHOK expected 

618 except NameError: # Python 3+ 

619 

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

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

622 ''' 

623 if start: 

624 r = v = start[0] 

625 else: 

626 r, v = 0, MISSING 

627 for v in xs: 

628 r = f(r, v) 

629 if v is MISSING: 

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

631 return r 

632 

633 

634def fremainder(x, y): 

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

636 

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

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

639 

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

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

642 

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

644 

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

646 project/geographiclib/>} and Python 3.7+ 

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

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

649 ''' 

650 # with Python 2.7.16 and 3.7.3 on macOS 10.13.6 and 

651 # with Python 3.10.2 on macOS 12.2.1 M1 arm64 native 

652 # fmod( 0, 360) == 0.0 

653 # fmod( 360, 360) == 0.0 

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

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

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

657 # however, using the % operator ... 

658 # 0 % 360 == 0 

659 # 360 % 360 == 0 

660 # 360.0 % 360 == 0.0 

661 # -0 % 360 == 0 

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

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

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

665 

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

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

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

669 if _isfinite(x): 

670 try: 

671 r = _remainder(x, y) if x else x 

672 except Exception as e: 

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

674 else: # handle x INF and NINF as NAN 

675 r = NAN 

676 return r 

677 

678 

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

680 from math import hypot # OK in Python 3.7- 

681 

682 def hypot_(*xs): 

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

684 

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

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

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

688 handled differently. 

689 

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

691 

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

693 

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

695 

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

697 

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

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

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

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

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

703 ''' 

704 h, x2 = _h_x2(xs) 

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

706 

707elif _sys_version_info2 < (3, 10): 

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

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

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

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

712 

713 def hypot(x, y): 

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

715 

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

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

718 

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

720 ''' 

721 if x: 

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

723 elif y: 

724 h = fabs(y) 

725 else: 

726 h = _0_0 

727 return h 

728 

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

730else: 

731 from math import hypot # PYCHOK in Python 3.10+ 

732 hypot_ = hypot 

733 

734 

735def _h_x2(xs): 

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

737 ''' 

738 if xs: 

739 n, xs = len2(xs) 

740 if n > 0: 

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

742 if h < EPS0: 

743 x2 = _0_0 

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

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

746 return h, x2 

747 

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

749 

750 

751def hypot1(x): 

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

753 

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

755 

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

757 ''' 

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

759 

760 

761def hypot2(x, y): 

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

763 

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

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

766 

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

768 ''' 

769 if x: 

770 if y: 

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

772 x, y = y, x 

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

774 else: 

775 h2 = x**2 

776 elif y: 

777 h2 = y**2 

778 else: 

779 h2 = _0_0 

780 return h2 

781 

782 

783def hypot2_(*xs): 

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

785 

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

787 

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

789 

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

791 

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

793 

794 @see: Function L{hypot_}. 

795 ''' 

796 h, x2 = _h_x2(xs) 

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

798 

799 

800def _map_a_x_b(a, b, where): 

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

802 ''' 

803 n = len(b) 

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

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

806 return map(_mul, a, b) if n > 3 else _map_a_x_b1(a, b) 

807 

808 

809def _map_a_x_b1(a, b): 

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

811 ''' 

812 yield _1_0 

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

814 if ab: 

815 yield ab 

816 yield _N_1_0 

817 

818 

819def norm2(x, y): 

820 '''Normalize a 2-dimensional vector. 

821 

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

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

824 

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

826 

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

828 or zero norm. 

829 ''' 

830 h = hypot(x, y) 

831 if not h: 

832 x = y = _0_0 # pass? 

833 elif not isnear1(h): 

834 try: 

835 x, y = x / h, y / h 

836 except Exception as e: 

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

838 return x, y 

839 

840 

841def norm_(*xs): 

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

843 

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

845 

846 @return: Yield each component, normalized. 

847 

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

849 or zero norm. 

850 ''' 

851 h = hypot_(*xs) 

852 if h: 

853 try: 

854 for i, x in enumerate(xs): 

855 yield x / h 

856 except Exception as e: 

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

858 else: 

859 for _ in xs: 

860 yield _0_0 

861 

862 

863def sqrt0(x): 

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

865 

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

867 

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

869 

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

871 returns C{0.0}. 

872 ''' 

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

874 

875 

876def sqrt3(x): 

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

878 

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

880 

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

882 

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

884 

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

886 ''' 

887 if x < 0: 

888 raise _ValueError(x=x, txt=_negative_) 

889 return pow(x, _1_5) if x else _0_0 

890 

891 

892def sqrt_a(h, b): 

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

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

895 

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

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

898 

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

900 

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

902 

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

904 

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

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

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

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

909 ''' 

910 try: 

911 if not (isscalar(h) and isscalar(b)): 

912 raise TypeError(_not_scalar_) 

913 elif isnear0(h): # PYCHOK no cover 

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

915 d = c - b 

916 if d < 0: 

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

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

919 else: 

920 c = float(h) 

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

922 if s < 0: 

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

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

925 except Exception as x: 

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

927 return a 

928 

929 

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

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

932 ''' 

933 yield s 

934 if h in (_0_0, _1_0): 

935 for x in xs: 

936 if x: 

937 yield x**2 

938 else: 

939 for x in xs: 

940 if x: 

941 yield (x / h)**2 

942 yield e 

943 

944# **) MIT License 

945# 

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

947# 

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

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

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

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

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

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

954# 

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

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

957# 

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

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

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

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

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

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

964# OTHER DEALINGS IN THE SOFTWARE.