Coverage for /usr/lib/python3/dist-packages/mpmath/ctx_iv.py: 49%

1import operator 

2 

3from . import libmp 

4 

5from .libmp.backend import basestring 

6 

7from .libmp import ( 

8 int_types, MPZ_ONE, 

9 prec_to_dps, dps_to_prec, repr_dps, 

10 round_floor, round_ceiling, 

11 fzero, finf, fninf, fnan, 

12 mpf_le, mpf_neg, 

13 from_int, from_float, from_str, from_rational, 

14 mpi_mid, mpi_delta, mpi_str, 

15 mpi_abs, mpi_pos, mpi_neg, mpi_add, mpi_sub, 

16 mpi_mul, mpi_div, mpi_pow_int, mpi_pow, 

17 mpi_from_str, 

18 mpci_pos, mpci_neg, mpci_add, mpci_sub, mpci_mul, mpci_div, mpci_pow, 

19 mpci_abs, mpci_pow, mpci_exp, mpci_log, 

20 ComplexResult, 

21 mpf_hash, mpc_hash) 

22from .matrices.matrices import _matrix 

23 

24mpi_zero = (fzero, fzero) 

25 

26from .ctx_base import StandardBaseContext 

27 

28new = object.__new__ 

29 

30def convert_mpf_(x, prec, rounding): 

31 if hasattr(x, "_mpf_"): return x._mpf_ 

32 if isinstance(x, int_types): return from_int(x, prec, rounding) 

33 if isinstance(x, float): return from_float(x, prec, rounding) 

34 if isinstance(x, basestring): return from_str(x, prec, rounding) 

35 raise NotImplementedError 

36 

37 

38class ivmpf(object): 

39 """ 

40 Interval arithmetic class. Precision is controlled by iv.prec. 

41 """ 

42 

43 def __new__(cls, x=0): 

44 return cls.ctx.convert(x) 

45 

46 def cast(self, cls, f_convert): 

47 a, b = self._mpi_ 

48 if a == b: 

49 return cls(f_convert(a)) 

50 raise ValueError 

51 

52 def __int__(self): 

53 return self.cast(int, libmp.to_int) 

54 

55 def __float__(self): 

56 return self.cast(float, libmp.to_float) 

57 

58 def __complex__(self): 

59 return self.cast(complex, libmp.to_float) 

60 

61 def __hash__(self): 

62 a, b = self._mpi_ 

63 if a == b: 

64 return mpf_hash(a) 

65 else: 

66 return hash(self._mpi_) 

67 

68 @property 

69 def real(self): return self 

70 

71 @property 

72 def imag(self): return self.ctx.zero 

73 

74 def conjugate(self): return self 

75 

76 @property 

77 def a(self): 

78 a, b = self._mpi_ 

79 return self.ctx.make_mpf((a, a)) 

80 

81 @property 

82 def b(self): 

83 a, b = self._mpi_ 

84 return self.ctx.make_mpf((b, b)) 

85 

86 @property 

87 def mid(self): 

88 ctx = self.ctx 

89 v = mpi_mid(self._mpi_, ctx.prec) 

90 return ctx.make_mpf((v, v)) 

91 

92 @property 

93 def delta(self): 

94 ctx = self.ctx 

95 v = mpi_delta(self._mpi_, ctx.prec) 

96 return ctx.make_mpf((v,v)) 

97 

98 @property 

99 def _mpci_(self): 

100 return self._mpi_, mpi_zero 

101 

102 def _compare(*args): 

103 raise TypeError("no ordering relation is defined for intervals") 

104 

105 __gt__ = _compare 

106 __le__ = _compare 

107 __gt__ = _compare 

108 __ge__ = _compare 

109 

110 def __contains__(self, t): 

111 t = self.ctx.mpf(t) 

112 return (self.a <= t.a) and (t.b <= self.b) 

113 

114 def __str__(self): 

115 return mpi_str(self._mpi_, self.ctx.prec) 

116 

117 def __repr__(self): 

118 if self.ctx.pretty: 

119 return str(self) 

120 a, b = self._mpi_ 

121 n = repr_dps(self.ctx.prec) 

122 a = libmp.to_str(a, n) 

123 b = libmp.to_str(b, n) 

124 return "mpi(%r, %r)" % (a, b) 

125 

126 def _compare(s, t, cmpfun): 

127 if not hasattr(t, "_mpi_"): 

128 try: 

129 t = s.ctx.convert(t) 

130 except: 

131 return NotImplemented 

132 return cmpfun(s._mpi_, t._mpi_) 

133 

134 def __eq__(s, t): return s._compare(t, libmp.mpi_eq) 

135 def __ne__(s, t): return s._compare(t, libmp.mpi_ne) 

136 def __lt__(s, t): return s._compare(t, libmp.mpi_lt) 

137 def __le__(s, t): return s._compare(t, libmp.mpi_le) 

138 def __gt__(s, t): return s._compare(t, libmp.mpi_gt) 

139 def __ge__(s, t): return s._compare(t, libmp.mpi_ge) 

140 

141 def __abs__(self): 

142 return self.ctx.make_mpf(mpi_abs(self._mpi_, self.ctx.prec)) 

143 def __pos__(self): 

144 return self.ctx.make_mpf(mpi_pos(self._mpi_, self.ctx.prec)) 

145 def __neg__(self): 

146 return self.ctx.make_mpf(mpi_neg(self._mpi_, self.ctx.prec)) 

147 

148 def ae(s, t, rel_eps=None, abs_eps=None): 

149 return s.ctx.almosteq(s, t, rel_eps, abs_eps) 

150 

151class ivmpc(object): 

152 

153 def __new__(cls, re=0, im=0): 

154 re = cls.ctx.convert(re) 

155 im = cls.ctx.convert(im) 

156 y = new(cls) 

157 y._mpci_ = re._mpi_, im._mpi_ 

158 return y 

159 

160 def __hash__(self): 

161 (a, b), (c,d) = self._mpci_ 

162 if a == b and c == d: 

163 return mpc_hash((a, c)) 

164 else: 

165 return hash(self._mpci_) 

166 

167 def __repr__(s): 

168 if s.ctx.pretty: 

169 return str(s) 

170 return "iv.mpc(%s, %s)" % (repr(s.real), repr(s.imag)) 

171 

172 def __str__(s): 

173 return "(%s + %s*j)" % (str(s.real), str(s.imag)) 

174 

175 @property 

176 def a(self): 

177 (a, b), (c,d) = self._mpci_ 

178 return self.ctx.make_mpf((a, a)) 

179 

180 @property 

181 def b(self): 

182 (a, b), (c,d) = self._mpci_ 

183 return self.ctx.make_mpf((b, b)) 

184 

185 @property 

186 def c(self): 

187 (a, b), (c,d) = self._mpci_ 

188 return self.ctx.make_mpf((c, c)) 

189 

190 @property 

191 def d(self): 

192 (a, b), (c,d) = self._mpci_ 

193 return self.ctx.make_mpf((d, d)) 

194 

195 @property 

196 def real(s): 

197 return s.ctx.make_mpf(s._mpci_[0]) 

198 

199 @property 

200 def imag(s): 

201 return s.ctx.make_mpf(s._mpci_[1]) 

202 

203 def conjugate(s): 

204 a, b = s._mpci_ 

205 return s.ctx.make_mpc((a, mpf_neg(b))) 

206 

207 def overlap(s, t): 

208 t = s.ctx.convert(t) 

209 real_overlap = (s.a <= t.a <= s.b) or (s.a <= t.b <= s.b) or (t.a <= s.a <= t.b) or (t.a <= s.b <= t.b) 

210 imag_overlap = (s.c <= t.c <= s.d) or (s.c <= t.d <= s.d) or (t.c <= s.c <= t.d) or (t.c <= s.d <= t.d) 

211 return real_overlap and imag_overlap 

212 

213 def __contains__(s, t): 

214 t = s.ctx.convert(t) 

215 return t.real in s.real and t.imag in s.imag 

216 

217 def _compare(s, t, ne=False): 

218 if not isinstance(t, s.ctx._types): 

219 try: 

220 t = s.ctx.convert(t) 

221 except: 

222 return NotImplemented 

223 if hasattr(t, '_mpi_'): 

224 tval = t._mpi_, mpi_zero 

225 elif hasattr(t, '_mpci_'): 

226 tval = t._mpci_ 

227 if ne: 

228 return s._mpci_ != tval 

229 return s._mpci_ == tval 

230 

231 def __eq__(s, t): return s._compare(t) 

232 def __ne__(s, t): return s._compare(t, True) 

233 

234 def __lt__(s, t): raise TypeError("complex intervals cannot be ordered") 

235 __le__ = __gt__ = __ge__ = __lt__ 

236 

237 def __neg__(s): return s.ctx.make_mpc(mpci_neg(s._mpci_, s.ctx.prec)) 

238 def __pos__(s): return s.ctx.make_mpc(mpci_pos(s._mpci_, s.ctx.prec)) 

239 def __abs__(s): return s.ctx.make_mpf(mpci_abs(s._mpci_, s.ctx.prec)) 

240 

241 def ae(s, t, rel_eps=None, abs_eps=None): 

242 return s.ctx.almosteq(s, t, rel_eps, abs_eps) 

243 

244def _binary_op(f_real, f_complex): 

245 def g_complex(ctx, sval, tval): 

246 return ctx.make_mpc(f_complex(sval, tval, ctx.prec)) 

247 def g_real(ctx, sval, tval): 

248 try: 

249 return ctx.make_mpf(f_real(sval, tval, ctx.prec)) 

250 except ComplexResult: 

251 sval = (sval, mpi_zero) 

252 tval = (tval, mpi_zero) 

253 return g_complex(ctx, sval, tval) 

254 def lop_real(s, t): 

255 if isinstance(t, _matrix): return NotImplemented 

256 ctx = s.ctx 

257 if not isinstance(t, ctx._types): t = ctx.convert(t) 

258 if hasattr(t, "_mpi_"): return g_real(ctx, s._mpi_, t._mpi_) 

259 if hasattr(t, "_mpci_"): return g_complex(ctx, (s._mpi_, mpi_zero), t._mpci_) 

260 return NotImplemented 

261 def rop_real(s, t): 

262 ctx = s.ctx 

263 if not isinstance(t, ctx._types): t = ctx.convert(t) 

264 if hasattr(t, "_mpi_"): return g_real(ctx, t._mpi_, s._mpi_) 

265 if hasattr(t, "_mpci_"): return g_complex(ctx, t._mpci_, (s._mpi_, mpi_zero)) 

266 return NotImplemented 

267 def lop_complex(s, t): 

268 if isinstance(t, _matrix): return NotImplemented 

269 ctx = s.ctx 

270 if not isinstance(t, s.ctx._types): 

271 try: 

272 t = s.ctx.convert(t) 

273 except (ValueError, TypeError): 

274 return NotImplemented 

275 return g_complex(ctx, s._mpci_, t._mpci_) 

276 def rop_complex(s, t): 

277 ctx = s.ctx 

278 if not isinstance(t, s.ctx._types): 

279 t = s.ctx.convert(t) 

280 return g_complex(ctx, t._mpci_, s._mpci_) 

281 return lop_real, rop_real, lop_complex, rop_complex 

282 

283ivmpf.__add__, ivmpf.__radd__, ivmpc.__add__, ivmpc.__radd__ = _binary_op(mpi_add, mpci_add) 

284ivmpf.__sub__, ivmpf.__rsub__, ivmpc.__sub__, ivmpc.__rsub__ = _binary_op(mpi_sub, mpci_sub) 

285ivmpf.__mul__, ivmpf.__rmul__, ivmpc.__mul__, ivmpc.__rmul__ = _binary_op(mpi_mul, mpci_mul) 

286ivmpf.__div__, ivmpf.__rdiv__, ivmpc.__div__, ivmpc.__rdiv__ = _binary_op(mpi_div, mpci_div) 

287ivmpf.__pow__, ivmpf.__rpow__, ivmpc.__pow__, ivmpc.__rpow__ = _binary_op(mpi_pow, mpci_pow) 

288 

289ivmpf.__truediv__ = ivmpf.__div__; ivmpf.__rtruediv__ = ivmpf.__rdiv__ 

290ivmpc.__truediv__ = ivmpc.__div__; ivmpc.__rtruediv__ = ivmpc.__rdiv__ 

291 

292class ivmpf_constant(ivmpf): 

293 def __new__(cls, f): 

294 self = new(cls) 

295 self._f = f 

296 return self 

297 def _get_mpi_(self): 

298 prec = self.ctx._prec[0] 

299 a = self._f(prec, round_floor) 

300 b = self._f(prec, round_ceiling) 

301 return a, b 

302 _mpi_ = property(_get_mpi_) 

303 

304class MPIntervalContext(StandardBaseContext): 

305 

306 def __init__(ctx): 

307 ctx.mpf = type('ivmpf', (ivmpf,), {}) 

308 ctx.mpc = type('ivmpc', (ivmpc,), {}) 

309 ctx._types = (ctx.mpf, ctx.mpc) 

310 ctx._constant = type('ivmpf_constant', (ivmpf_constant,), {}) 

311 ctx._prec = [53] 

312 ctx._set_prec(53) 

313 ctx._constant._ctxdata = ctx.mpf._ctxdata = ctx.mpc._ctxdata = [ctx.mpf, new, ctx._prec] 

314 ctx._constant.ctx = ctx.mpf.ctx = ctx.mpc.ctx = ctx 

315 ctx.pretty = False 

316 StandardBaseContext.__init__(ctx) 

317 ctx._init_builtins() 

318 

319 def _mpi(ctx, a, b=None): 

320 if b is None: 

321 return ctx.mpf(a) 

322 return ctx.mpf((a,b)) 

323 

324 def _init_builtins(ctx): 

325 ctx.one = ctx.mpf(1) 

326 ctx.zero = ctx.mpf(0) 

327 ctx.inf = ctx.mpf('inf') 

328 ctx.ninf = -ctx.inf 

329 ctx.nan = ctx.mpf('nan') 

330 ctx.j = ctx.mpc(0,1) 

331 ctx.exp = ctx._wrap_mpi_function(libmp.mpi_exp, libmp.mpci_exp) 

332 ctx.sqrt = ctx._wrap_mpi_function(libmp.mpi_sqrt) 

333 ctx.ln = ctx._wrap_mpi_function(libmp.mpi_log, libmp.mpci_log) 

334 ctx.cos = ctx._wrap_mpi_function(libmp.mpi_cos, libmp.mpci_cos) 

335 ctx.sin = ctx._wrap_mpi_function(libmp.mpi_sin, libmp.mpci_sin) 

336 ctx.tan = ctx._wrap_mpi_function(libmp.mpi_tan) 

337 ctx.gamma = ctx._wrap_mpi_function(libmp.mpi_gamma, libmp.mpci_gamma) 

338 ctx.loggamma = ctx._wrap_mpi_function(libmp.mpi_loggamma, libmp.mpci_loggamma) 

339 ctx.rgamma = ctx._wrap_mpi_function(libmp.mpi_rgamma, libmp.mpci_rgamma) 

340 ctx.factorial = ctx._wrap_mpi_function(libmp.mpi_factorial, libmp.mpci_factorial) 

341 ctx.fac = ctx.factorial 

342 

343 ctx.eps = ctx._constant(lambda prec, rnd: (0, MPZ_ONE, 1-prec, 1)) 

344 ctx.pi = ctx._constant(libmp.mpf_pi) 

345 ctx.e = ctx._constant(libmp.mpf_e) 

346 ctx.ln2 = ctx._constant(libmp.mpf_ln2) 

347 ctx.ln10 = ctx._constant(libmp.mpf_ln10) 

348 ctx.phi = ctx._constant(libmp.mpf_phi) 

349 ctx.euler = ctx._constant(libmp.mpf_euler) 

350 ctx.catalan = ctx._constant(libmp.mpf_catalan) 

351 ctx.glaisher = ctx._constant(libmp.mpf_glaisher) 

352 ctx.khinchin = ctx._constant(libmp.mpf_khinchin) 

353 ctx.twinprime = ctx._constant(libmp.mpf_twinprime) 

354 

355 def _wrap_mpi_function(ctx, f_real, f_complex=None): 

356 def g(x, **kwargs): 

357 if kwargs: 

358 prec = kwargs.get('prec', ctx._prec[0]) 

359 else: 

360 prec = ctx._prec[0] 

361 x = ctx.convert(x) 

362 if hasattr(x, "_mpi_"): 

363 return ctx.make_mpf(f_real(x._mpi_, prec)) 

364 if hasattr(x, "_mpci_"): 

365 return ctx.make_mpc(f_complex(x._mpci_, prec)) 

366 raise ValueError 

367 return g 

368 

369 @classmethod 

370 def _wrap_specfun(cls, name, f, wrap): 

371 if wrap: 

372 def f_wrapped(ctx, *args, **kwargs): 

373 convert = ctx.convert 

374 args = [convert(a) for a in args] 

375 prec = ctx.prec 

376 try: 

377 ctx.prec += 10 

378 retval = f(ctx, *args, **kwargs) 

379 finally: 

380 ctx.prec = prec 

381 return +retval 

382 else: 

383 f_wrapped = f 

384 setattr(cls, name, f_wrapped) 

385 

386 def _set_prec(ctx, n): 

387 ctx._prec[0] = max(1, int(n)) 

388 ctx._dps = prec_to_dps(n) 

389 

390 def _set_dps(ctx, n): 

391 ctx._prec[0] = dps_to_prec(n) 

392 ctx._dps = max(1, int(n)) 

393 

394 prec = property(lambda ctx: ctx._prec[0], _set_prec) 

395 dps = property(lambda ctx: ctx._dps, _set_dps) 

396 

397 def make_mpf(ctx, v): 

398 a = new(ctx.mpf) 

399 a._mpi_ = v 

400 return a 

401 

402 def make_mpc(ctx, v): 

403 a = new(ctx.mpc) 

404 a._mpci_ = v 

405 return a 

406 

407 def _mpq(ctx, pq): 

408 p, q = pq 

409 a = libmp.from_rational(p, q, ctx.prec, round_floor) 

410 b = libmp.from_rational(p, q, ctx.prec, round_ceiling) 

411 return ctx.make_mpf((a, b)) 

412 

413 def convert(ctx, x): 

414 if isinstance(x, (ctx.mpf, ctx.mpc)): 

415 return x 

416 if isinstance(x, ctx._constant): 

417 return +x 

418 if isinstance(x, complex) or hasattr(x, "_mpc_"): 

419 re = ctx.convert(x.real) 

420 im = ctx.convert(x.imag) 

421 return ctx.mpc(re,im) 

422 if isinstance(x, basestring): 

423 v = mpi_from_str(x, ctx.prec) 

424 return ctx.make_mpf(v) 

425 if hasattr(x, "_mpi_"): 

426 a, b = x._mpi_ 

427 else: 

428 try: 

429 a, b = x 

430 except (TypeError, ValueError): 

431 a = b = x 

432 if hasattr(a, "_mpi_"): 

433 a = a._mpi_[0] 

434 else: 

435 a = convert_mpf_(a, ctx.prec, round_floor) 

436 if hasattr(b, "_mpi_"): 

437 b = b._mpi_[1] 

438 else: 

439 b = convert_mpf_(b, ctx.prec, round_ceiling) 

440 if a == fnan or b == fnan: 

441 a = fninf 

442 b = finf 

443 assert mpf_le(a, b), "endpoints must be properly ordered" 

444 return ctx.make_mpf((a, b)) 

445 

446 def nstr(ctx, x, n=5, **kwargs): 

447 x = ctx.convert(x) 

448 if hasattr(x, "_mpi_"): 

449 return libmp.mpi_to_str(x._mpi_, n, **kwargs) 

450 if hasattr(x, "_mpci_"): 

451 re = libmp.mpi_to_str(x._mpci_[0], n, **kwargs) 

452 im = libmp.mpi_to_str(x._mpci_[1], n, **kwargs) 

453 return "(%s + %s*j)" % (re, im) 

454 

455 def mag(ctx, x): 

456 x = ctx.convert(x) 

457 if isinstance(x, ctx.mpc): 

458 return max(ctx.mag(x.real), ctx.mag(x.imag)) + 1 

459 a, b = libmp.mpi_abs(x._mpi_) 

460 sign, man, exp, bc = b 

461 if man: 

462 return exp+bc 

463 if b == fzero: 

464 return ctx.ninf 

465 if b == fnan: 

466 return ctx.nan 

467 return ctx.inf 

468 

469 def isnan(ctx, x): 

470 return False 

471 

472 def isinf(ctx, x): 

473 return x == ctx.inf 

474 

475 def isint(ctx, x): 

476 x = ctx.convert(x) 

477 a, b = x._mpi_ 

478 if a == b: 

479 sign, man, exp, bc = a 

480 if man: 

481 return exp >= 0 

482 return a == fzero 

483 return None 

484 

485 def ldexp(ctx, x, n): 

486 a, b = ctx.convert(x)._mpi_ 

487 a = libmp.mpf_shift(a, n) 

488 b = libmp.mpf_shift(b, n) 

489 return ctx.make_mpf((a,b)) 

490 

491 def absmin(ctx, x): 

492 return abs(ctx.convert(x)).a 

493 

494 def absmax(ctx, x): 

495 return abs(ctx.convert(x)).b 

496 

497 def atan2(ctx, y, x): 

498 y = ctx.convert(y)._mpi_ 

499 x = ctx.convert(x)._mpi_ 

500 return ctx.make_mpf(libmp.mpi_atan2(y,x,ctx.prec)) 

501 

502 def _convert_param(ctx, x): 

503 if isinstance(x, libmp.int_types): 

504 return x, 'Z' 

505 if isinstance(x, tuple): 

506 p, q = x 

507 return (ctx.mpf(p) / ctx.mpf(q), 'R') 

508 x = ctx.convert(x) 

509 if isinstance(x, ctx.mpf): 

510 return x, 'R' 

511 if isinstance(x, ctx.mpc): 

512 return x, 'C' 

513 raise ValueError 

514 

515 def _is_real_type(ctx, z): 

516 return isinstance(z, ctx.mpf) or isinstance(z, int_types) 

517 

518 def _is_complex_type(ctx, z): 

519 return isinstance(z, ctx.mpc) 

520 

521 def hypsum(ctx, p, q, types, coeffs, z, maxterms=6000, **kwargs): 

522 coeffs = list(coeffs) 

523 num = range(p) 

524 den = range(p,p+q) 

525 #tol = ctx.eps 

526 s = t = ctx.one 

527 k = 0 

528 while 1: 

529 for i in num: t *= (coeffs[i]+k) 

530 for i in den: t /= (coeffs[i]+k) 

531 k += 1; t /= k; t *= z; s += t 

532 if t == 0: 

533 return s 

534 #if abs(t) < tol: 

535 # return s 

536 if k > maxterms: 

537 raise ctx.NoConvergence 

538 

539 

540# Register with "numbers" ABC 

541# We do not subclass, hence we do not use the @abstractmethod checks. While 

542# this is less invasive it may turn out that we do not actually support 

543# parts of the expected interfaces. See 

544# http://docs.python.org/2/library/numbers.html for list of abstract 

545# methods. 

546try: 

547 import numbers 

548 numbers.Complex.register(ivmpc) 

549 numbers.Real.register(ivmpf) 

550except ImportError: 

551 pass