Coverage for /usr/lib/python3/dist-packages/mpmath/libmp/libmpi.py: 13%

1""" 

2Computational functions for interval arithmetic. 

3 

4""" 

5 

6from .backend import xrange 

7 

8from .libmpf import ( 

9 ComplexResult, 

10 round_down, round_up, round_floor, round_ceiling, round_nearest, 

11 prec_to_dps, repr_dps, dps_to_prec, 

12 bitcount, 

13 from_float, 

14 fnan, finf, fninf, fzero, fhalf, fone, fnone, 

15 mpf_sign, mpf_lt, mpf_le, mpf_gt, mpf_ge, mpf_eq, mpf_cmp, 

16 mpf_min_max, 

17 mpf_floor, from_int, to_int, to_str, from_str, 

18 mpf_abs, mpf_neg, mpf_pos, mpf_add, mpf_sub, mpf_mul, mpf_mul_int, 

19 mpf_div, mpf_shift, mpf_pow_int, 

20 from_man_exp, MPZ_ONE) 

21 

22from .libelefun import ( 

23 mpf_log, mpf_exp, mpf_sqrt, mpf_atan, mpf_atan2, 

24 mpf_pi, mod_pi2, mpf_cos_sin 

25) 

26 

27from .gammazeta import mpf_gamma, mpf_rgamma, mpf_loggamma, mpc_loggamma 

28 

29def mpi_str(s, prec): 

30 sa, sb = s 

31 dps = prec_to_dps(prec) + 5 

32 return "[%s, %s]" % (to_str(sa, dps), to_str(sb, dps)) 

33 #dps = prec_to_dps(prec) 

34 #m = mpi_mid(s, prec) 

35 #d = mpf_shift(mpi_delta(s, 20), -1) 

36 #return "%s +/- %s" % (to_str(m, dps), to_str(d, 3)) 

37 

38mpi_zero = (fzero, fzero) 

39mpi_one = (fone, fone) 

40 

41def mpi_eq(s, t): 

42 return s == t 

43 

44def mpi_ne(s, t): 

45 return s != t 

46 

47def mpi_lt(s, t): 

48 sa, sb = s 

49 ta, tb = t 

50 if mpf_lt(sb, ta): return True 

51 if mpf_ge(sa, tb): return False 

52 return None 

53 

54def mpi_le(s, t): 

55 sa, sb = s 

56 ta, tb = t 

57 if mpf_le(sb, ta): return True 

58 if mpf_gt(sa, tb): return False 

59 return None 

60 

61def mpi_gt(s, t): return mpi_lt(t, s) 

62def mpi_ge(s, t): return mpi_le(t, s) 

63 

64def mpi_add(s, t, prec=0): 

65 sa, sb = s 

66 ta, tb = t 

67 a = mpf_add(sa, ta, prec, round_floor) 

68 b = mpf_add(sb, tb, prec, round_ceiling) 

69 if a == fnan: a = fninf 

70 if b == fnan: b = finf 

71 return a, b 

72 

73def mpi_sub(s, t, prec=0): 

74 sa, sb = s 

75 ta, tb = t 

76 a = mpf_sub(sa, tb, prec, round_floor) 

77 b = mpf_sub(sb, ta, prec, round_ceiling) 

78 if a == fnan: a = fninf 

79 if b == fnan: b = finf 

80 return a, b 

81 

82def mpi_delta(s, prec): 

83 sa, sb = s 

84 return mpf_sub(sb, sa, prec, round_up) 

85 

86def mpi_mid(s, prec): 

87 sa, sb = s 

88 return mpf_shift(mpf_add(sa, sb, prec, round_nearest), -1) 

89 

90def mpi_pos(s, prec): 

91 sa, sb = s 

92 a = mpf_pos(sa, prec, round_floor) 

93 b = mpf_pos(sb, prec, round_ceiling) 

94 return a, b 

95 

96def mpi_neg(s, prec=0): 

97 sa, sb = s 

98 a = mpf_neg(sb, prec, round_floor) 

99 b = mpf_neg(sa, prec, round_ceiling) 

100 return a, b 

101 

102def mpi_abs(s, prec=0): 

103 sa, sb = s 

104 sas = mpf_sign(sa) 

105 sbs = mpf_sign(sb) 

106 # Both points nonnegative? 

107 if sas >= 0: 

108 a = mpf_pos(sa, prec, round_floor) 

109 b = mpf_pos(sb, prec, round_ceiling) 

110 # Upper point nonnegative? 

111 elif sbs >= 0: 

112 a = fzero 

113 negsa = mpf_neg(sa) 

114 if mpf_lt(negsa, sb): 

115 b = mpf_pos(sb, prec, round_ceiling) 

116 else: 

117 b = mpf_pos(negsa, prec, round_ceiling) 

118 # Both negative? 

119 else: 

120 a = mpf_neg(sb, prec, round_floor) 

121 b = mpf_neg(sa, prec, round_ceiling) 

122 return a, b 

123 

124# TODO: optimize 

125def mpi_mul_mpf(s, t, prec): 

126 return mpi_mul(s, (t, t), prec) 

127 

128def mpi_div_mpf(s, t, prec): 

129 return mpi_div(s, (t, t), prec) 

130 

131def mpi_mul(s, t, prec=0): 

132 sa, sb = s 

133 ta, tb = t 

134 sas = mpf_sign(sa) 

135 sbs = mpf_sign(sb) 

136 tas = mpf_sign(ta) 

137 tbs = mpf_sign(tb) 

138 if sas == sbs == 0: 

139 # Should maybe be undefined 

140 if ta == fninf or tb == finf: 

141 return fninf, finf 

142 return fzero, fzero 

143 if tas == tbs == 0: 

144 # Should maybe be undefined 

145 if sa == fninf or sb == finf: 

146 return fninf, finf 

147 return fzero, fzero 

148 if sas >= 0: 

149 # positive * positive 

150 if tas >= 0: 

151 a = mpf_mul(sa, ta, prec, round_floor) 

152 b = mpf_mul(sb, tb, prec, round_ceiling) 

153 if a == fnan: a = fzero 

154 if b == fnan: b = finf 

155 # positive * negative 

156 elif tbs <= 0: 

157 a = mpf_mul(sb, ta, prec, round_floor) 

158 b = mpf_mul(sa, tb, prec, round_ceiling) 

159 if a == fnan: a = fninf 

160 if b == fnan: b = fzero 

161 # positive * both signs 

162 else: 

163 a = mpf_mul(sb, ta, prec, round_floor) 

164 b = mpf_mul(sb, tb, prec, round_ceiling) 

165 if a == fnan: a = fninf 

166 if b == fnan: b = finf 

167 elif sbs <= 0: 

168 # negative * positive 

169 if tas >= 0: 

170 a = mpf_mul(sa, tb, prec, round_floor) 

171 b = mpf_mul(sb, ta, prec, round_ceiling) 

172 if a == fnan: a = fninf 

173 if b == fnan: b = fzero 

174 # negative * negative 

175 elif tbs <= 0: 

176 a = mpf_mul(sb, tb, prec, round_floor) 

177 b = mpf_mul(sa, ta, prec, round_ceiling) 

178 if a == fnan: a = fzero 

179 if b == fnan: b = finf 

180 # negative * both signs 

181 else: 

182 a = mpf_mul(sa, tb, prec, round_floor) 

183 b = mpf_mul(sa, ta, prec, round_ceiling) 

184 if a == fnan: a = fninf 

185 if b == fnan: b = finf 

186 else: 

187 # General case: perform all cross-multiplications and compare 

188 # Since the multiplications can be done exactly, we need only 

189 # do 4 (instead of 8: two for each rounding mode) 

190 cases = [mpf_mul(sa, ta), mpf_mul(sa, tb), mpf_mul(sb, ta), mpf_mul(sb, tb)] 

191 if fnan in cases: 

192 a, b = (fninf, finf) 

193 else: 

194 a, b = mpf_min_max(cases) 

195 a = mpf_pos(a, prec, round_floor) 

196 b = mpf_pos(b, prec, round_ceiling) 

197 return a, b 

198 

199def mpi_square(s, prec=0): 

200 sa, sb = s 

201 if mpf_ge(sa, fzero): 

202 a = mpf_mul(sa, sa, prec, round_floor) 

203 b = mpf_mul(sb, sb, prec, round_ceiling) 

204 elif mpf_le(sb, fzero): 

205 a = mpf_mul(sb, sb, prec, round_floor) 

206 b = mpf_mul(sa, sa, prec, round_ceiling) 

207 else: 

208 sa = mpf_neg(sa) 

209 sa, sb = mpf_min_max([sa, sb]) 

210 a = fzero 

211 b = mpf_mul(sb, sb, prec, round_ceiling) 

212 return a, b 

213 

214def mpi_div(s, t, prec): 

215 sa, sb = s 

216 ta, tb = t 

217 sas = mpf_sign(sa) 

218 sbs = mpf_sign(sb) 

219 tas = mpf_sign(ta) 

220 tbs = mpf_sign(tb) 

221 # 0 / X 

222 if sas == sbs == 0: 

223 # 0 / <interval containing 0> 

224 if (tas < 0 and tbs > 0) or (tas == 0 or tbs == 0): 

225 return fninf, finf 

226 return fzero, fzero 

227 # Denominator contains both negative and positive numbers; 

228 # this should properly be a multi-interval, but the closest 

229 # match is the entire (extended) real line 

230 if tas < 0 and tbs > 0: 

231 return fninf, finf 

232 # Assume denominator to be nonnegative 

233 if tas < 0: 

234 return mpi_div(mpi_neg(s), mpi_neg(t), prec) 

235 # Division by zero 

236 # XXX: make sure all results make sense 

237 if tas == 0: 

238 # Numerator contains both signs? 

239 if sas < 0 and sbs > 0: 

240 return fninf, finf 

241 if tas == tbs: 

242 return fninf, finf 

243 # Numerator positive? 

244 if sas >= 0: 

245 a = mpf_div(sa, tb, prec, round_floor) 

246 b = finf 

247 if sbs <= 0: 

248 a = fninf 

249 b = mpf_div(sb, tb, prec, round_ceiling) 

250 # Division with positive denominator 

251 # We still have to handle nans resulting from inf/0 or inf/inf 

252 else: 

253 # Nonnegative numerator 

254 if sas >= 0: 

255 a = mpf_div(sa, tb, prec, round_floor) 

256 b = mpf_div(sb, ta, prec, round_ceiling) 

257 if a == fnan: a = fzero 

258 if b == fnan: b = finf 

259 # Nonpositive numerator 

260 elif sbs <= 0: 

261 a = mpf_div(sa, ta, prec, round_floor) 

262 b = mpf_div(sb, tb, prec, round_ceiling) 

263 if a == fnan: a = fninf 

264 if b == fnan: b = fzero 

265 # Numerator contains both signs? 

266 else: 

267 a = mpf_div(sa, ta, prec, round_floor) 

268 b = mpf_div(sb, ta, prec, round_ceiling) 

269 if a == fnan: a = fninf 

270 if b == fnan: b = finf 

271 return a, b 

272 

273def mpi_pi(prec): 

274 a = mpf_pi(prec, round_floor) 

275 b = mpf_pi(prec, round_ceiling) 

276 return a, b 

277 

278def mpi_exp(s, prec): 

279 sa, sb = s 

280 # exp is monotonic 

281 a = mpf_exp(sa, prec, round_floor) 

282 b = mpf_exp(sb, prec, round_ceiling) 

283 return a, b 

284 

285def mpi_log(s, prec): 

286 sa, sb = s 

287 # log is monotonic 

288 a = mpf_log(sa, prec, round_floor) 

289 b = mpf_log(sb, prec, round_ceiling) 

290 return a, b 

291 

292def mpi_sqrt(s, prec): 

293 sa, sb = s 

294 # sqrt is monotonic 

295 a = mpf_sqrt(sa, prec, round_floor) 

296 b = mpf_sqrt(sb, prec, round_ceiling) 

297 return a, b 

298 

299def mpi_atan(s, prec): 

300 sa, sb = s 

301 a = mpf_atan(sa, prec, round_floor) 

302 b = mpf_atan(sb, prec, round_ceiling) 

303 return a, b 

304 

305def mpi_pow_int(s, n, prec): 

306 sa, sb = s 

307 if n < 0: 

308 return mpi_div((fone, fone), mpi_pow_int(s, -n, prec+20), prec) 

309 if n == 0: 

310 return (fone, fone) 

311 if n == 1: 

312 return s 

313 if n == 2: 

314 return mpi_square(s, prec) 

315 # Odd -- signs are preserved 

316 if n & 1: 

317 a = mpf_pow_int(sa, n, prec, round_floor) 

318 b = mpf_pow_int(sb, n, prec, round_ceiling) 

319 # Even -- important to ensure positivity 

320 else: 

321 sas = mpf_sign(sa) 

322 sbs = mpf_sign(sb) 

323 # Nonnegative? 

324 if sas >= 0: 

325 a = mpf_pow_int(sa, n, prec, round_floor) 

326 b = mpf_pow_int(sb, n, prec, round_ceiling) 

327 # Nonpositive? 

328 elif sbs <= 0: 

329 a = mpf_pow_int(sb, n, prec, round_floor) 

330 b = mpf_pow_int(sa, n, prec, round_ceiling) 

331 # Mixed signs? 

332 else: 

333 a = fzero 

334 # max(-a,b)**n 

335 sa = mpf_neg(sa) 

336 if mpf_ge(sa, sb): 

337 b = mpf_pow_int(sa, n, prec, round_ceiling) 

338 else: 

339 b = mpf_pow_int(sb, n, prec, round_ceiling) 

340 return a, b 

341 

342def mpi_pow(s, t, prec): 

343 ta, tb = t 

344 if ta == tb and ta not in (finf, fninf): 

345 if ta == from_int(to_int(ta)): 

346 return mpi_pow_int(s, to_int(ta), prec) 

347 if ta == fhalf: 

348 return mpi_sqrt(s, prec) 

349 u = mpi_log(s, prec + 20) 

350 v = mpi_mul(u, t, prec + 20) 

351 return mpi_exp(v, prec) 

352 

353def MIN(x, y): 

354 if mpf_le(x, y): 

355 return x 

356 return y 

357 

358def MAX(x, y): 

359 if mpf_ge(x, y): 

360 return x 

361 return y 

362 

363def cos_sin_quadrant(x, wp): 

364 sign, man, exp, bc = x 

365 if x == fzero: 

366 return fone, fzero, 0 

367 # TODO: combine evaluation code to avoid duplicate modulo 

368 c, s = mpf_cos_sin(x, wp) 

369 t, n, wp_ = mod_pi2(man, exp, exp+bc, 15) 

370 if sign: 

371 n = -1-n 

372 return c, s, n 

373 

374def mpi_cos_sin(x, prec): 

375 a, b = x 

376 if a == b == fzero: 

377 return (fone, fone), (fzero, fzero) 

378 # Guaranteed to contain both -1 and 1 

379 if (finf in x) or (fninf in x): 

380 return (fnone, fone), (fnone, fone) 

381 wp = prec + 20 

382 ca, sa, na = cos_sin_quadrant(a, wp) 

383 cb, sb, nb = cos_sin_quadrant(b, wp) 

384 ca, cb = mpf_min_max([ca, cb]) 

385 sa, sb = mpf_min_max([sa, sb]) 

386 # Both functions are monotonic within one quadrant 

387 if na == nb: 

388 pass 

389 # Guaranteed to contain both -1 and 1 

390 elif nb - na >= 4: 

391 return (fnone, fone), (fnone, fone) 

392 else: 

393 # cos has maximum between a and b 

394 if na//4 != nb//4: 

395 cb = fone 

396 # cos has minimum 

397 if (na-2)//4 != (nb-2)//4: 

398 ca = fnone 

399 # sin has maximum 

400 if (na-1)//4 != (nb-1)//4: 

401 sb = fone 

402 # sin has minimum 

403 if (na-3)//4 != (nb-3)//4: 

404 sa = fnone 

405 # Perturb to force interval rounding 

406 more = from_man_exp((MPZ_ONE<<wp) + (MPZ_ONE<<10), -wp) 

407 less = from_man_exp((MPZ_ONE<<wp) - (MPZ_ONE<<10), -wp) 

408 def finalize(v, rounding): 

409 if bool(v[0]) == (rounding == round_floor): 

410 p = more 

411 else: 

412 p = less 

413 v = mpf_mul(v, p, prec, rounding) 

414 sign, man, exp, bc = v 

415 if exp+bc >= 1: 

416 if sign: 

417 return fnone 

418 return fone 

419 return v 

420 ca = finalize(ca, round_floor) 

421 cb = finalize(cb, round_ceiling) 

422 sa = finalize(sa, round_floor) 

423 sb = finalize(sb, round_ceiling) 

424 return (ca,cb), (sa,sb) 

425 

426def mpi_cos(x, prec): 

427 return mpi_cos_sin(x, prec)[0] 

428 

429def mpi_sin(x, prec): 

430 return mpi_cos_sin(x, prec)[1] 

431 

432def mpi_tan(x, prec): 

433 cos, sin = mpi_cos_sin(x, prec+20) 

434 return mpi_div(sin, cos, prec) 

435 

436def mpi_cot(x, prec): 

437 cos, sin = mpi_cos_sin(x, prec+20) 

438 return mpi_div(cos, sin, prec) 

439 

440def mpi_from_str_a_b(x, y, percent, prec): 

441 wp = prec + 20 

442 xa = from_str(x, wp, round_floor) 

443 xb = from_str(x, wp, round_ceiling) 

444 #ya = from_str(y, wp, round_floor) 

445 y = from_str(y, wp, round_ceiling) 

446 assert mpf_ge(y, fzero) 

447 if percent: 

448 y = mpf_mul(MAX(mpf_abs(xa), mpf_abs(xb)), y, wp, round_ceiling) 

449 y = mpf_div(y, from_int(100), wp, round_ceiling) 

450 a = mpf_sub(xa, y, prec, round_floor) 

451 b = mpf_add(xb, y, prec, round_ceiling) 

452 return a, b 

453 

454def mpi_from_str(s, prec): 

455 """ 

456 Parse an interval number given as a string. 

457 

458 Allowed forms are 

459 

460 "-1.23e-27" 

461 Any single decimal floating-point literal. 

462 "a +- b" or "a (b)" 

463 a is the midpoint of the interval and b is the half-width 

464 "a +- b%" or "a (b%)" 

465 a is the midpoint of the interval and the half-width 

466 is b percent of a (`a \times b / 100`). 

467 "[a, b]" 

468 The interval indicated directly. 

469 "x[y,z]e" 

470 x are shared digits, y and z are unequal digits, e is the exponent. 

471 

472 """ 

473 e = ValueError("Improperly formed interval number '%s'" % s) 

474 s = s.replace(" ", "") 

475 wp = prec + 20 

476 if "+-" in s: 

477 x, y = s.split("+-") 

478 return mpi_from_str_a_b(x, y, False, prec) 

479 # case 2 

480 elif "(" in s: 

481 # Don't confuse with a complex number (x,y) 

482 if s[0] == "(" or ")" not in s: 

483 raise e 

484 s = s.replace(")", "") 

485 percent = False 

486 if "%" in s: 

487 if s[-1] != "%": 

488 raise e 

489 percent = True 

490 s = s.replace("%", "") 

491 x, y = s.split("(") 

492 return mpi_from_str_a_b(x, y, percent, prec) 

493 elif "," in s: 

494 if ('[' not in s) or (']' not in s): 

495 raise e 

496 if s[0] == '[': 

497 # case 3 

498 s = s.replace("[", "") 

499 s = s.replace("]", "") 

500 a, b = s.split(",") 

501 a = from_str(a, prec, round_floor) 

502 b = from_str(b, prec, round_ceiling) 

503 return a, b 

504 else: 

505 # case 4 

506 x, y = s.split('[') 

507 y, z = y.split(',') 

508 if 'e' in s: 

509 z, e = z.split(']') 

510 else: 

511 z, e = z.rstrip(']'), '' 

512 a = from_str(x+y+e, prec, round_floor) 

513 b = from_str(x+z+e, prec, round_ceiling) 

514 return a, b 

515 else: 

516 a = from_str(s, prec, round_floor) 

517 b = from_str(s, prec, round_ceiling) 

518 return a, b 

519 

520def mpi_to_str(x, dps, use_spaces=True, brackets='[]', mode='brackets', error_dps=4, **kwargs): 

521 """ 

522 Convert a mpi interval to a string. 

523 

524 **Arguments** 

525 

526 *dps* 

527 decimal places to use for printing 

528 *use_spaces* 

529 use spaces for more readable output, defaults to true 

530 *brackets* 

531 pair of strings (or two-character string) giving left and right brackets 

532 *mode* 

533 mode of display: 'plusminus', 'percent', 'brackets' (default) or 'diff' 

534 *error_dps* 

535 limit the error to *error_dps* digits (mode 'plusminus and 'percent') 

536 

537 Additional keyword arguments are forwarded to the mpf-to-string conversion 

538 for the components of the output. 

539 

540 **Examples** 

541 

542 >>> from mpmath import mpi, mp 

543 >>> mp.dps = 30 

544 >>> x = mpi(1, 2)._mpi_ 

545 >>> mpi_to_str(x, 2, mode='plusminus') 

546 '1.5 +- 0.5' 

547 >>> mpi_to_str(x, 2, mode='percent') 

548 '1.5 (33.33%)' 

549 >>> mpi_to_str(x, 2, mode='brackets') 

550 '[1.0, 2.0]' 

551 >>> mpi_to_str(x, 2, mode='brackets' , brackets=('<', '>')) 

552 '<1.0, 2.0>' 

553 >>> x = mpi('5.2582327113062393041', '5.2582327113062749951')._mpi_ 

554 >>> mpi_to_str(x, 15, mode='diff') 

555 '5.2582327113062[4, 7]' 

556 >>> mpi_to_str(mpi(0)._mpi_, 2, mode='percent') 

557 '0.0 (0.0%)' 

558 

559 """ 

560 prec = dps_to_prec(dps) 

561 wp = prec + 20 

562 a, b = x 

563 mid = mpi_mid(x, prec) 

564 delta = mpi_delta(x, prec) 

565 a_str = to_str(a, dps, **kwargs) 

566 b_str = to_str(b, dps, **kwargs) 

567 mid_str = to_str(mid, dps, **kwargs) 

568 sp = "" 

569 if use_spaces: 

570 sp = " " 

571 br1, br2 = brackets 

572 if mode == 'plusminus': 

573 delta_str = to_str(mpf_shift(delta,-1), dps, **kwargs) 

574 s = mid_str + sp + "+-" + sp + delta_str 

575 elif mode == 'percent': 

576 if mid == fzero: 

577 p = fzero 

578 else: 

579 # p = 100 * delta(x) / (2*mid(x)) 

580 p = mpf_mul(delta, from_int(100)) 

581 p = mpf_div(p, mpf_mul(mid, from_int(2)), wp) 

582 s = mid_str + sp + "(" + to_str(p, error_dps) + "%)" 

583 elif mode == 'brackets': 

584 s = br1 + a_str + "," + sp + b_str + br2 

585 elif mode == 'diff': 

586 # use more digits if str(x.a) and str(x.b) are equal 

587 if a_str == b_str: 

588 a_str = to_str(a, dps+3, **kwargs) 

589 b_str = to_str(b, dps+3, **kwargs) 

590 # separate mantissa and exponent 

591 a = a_str.split('e') 

592 if len(a) == 1: 

593 a.append('') 

594 b = b_str.split('e') 

595 if len(b) == 1: 

596 b.append('') 

597 if a[1] == b[1]: 

598 if a[0] != b[0]: 

599 for i in xrange(len(a[0]) + 1): 

600 if a[0][i] != b[0][i]: 

601 break 

602 s = (a[0][:i] + br1 + a[0][i:] + ',' + sp + b[0][i:] + br2 

603 + 'e'*min(len(a[1]), 1) + a[1]) 

604 else: # no difference 

605 s = a[0] + br1 + br2 + 'e'*min(len(a[1]), 1) + a[1] 

606 else: 

607 s = br1 + 'e'.join(a) + ',' + sp + 'e'.join(b) + br2 

608 else: 

609 raise ValueError("'%s' is unknown mode for printing mpi" % mode) 

610 return s 

611 

612def mpci_add(x, y, prec): 

613 a, b = x 

614 c, d = y 

615 return mpi_add(a, c, prec), mpi_add(b, d, prec) 

616 

617def mpci_sub(x, y, prec): 

618 a, b = x 

619 c, d = y 

620 return mpi_sub(a, c, prec), mpi_sub(b, d, prec) 

621 

622def mpci_neg(x, prec=0): 

623 a, b = x 

624 return mpi_neg(a, prec), mpi_neg(b, prec) 

625 

626def mpci_pos(x, prec): 

627 a, b = x 

628 return mpi_pos(a, prec), mpi_pos(b, prec) 

629 

630def mpci_mul(x, y, prec): 

631 # TODO: optimize for real/imag cases 

632 a, b = x 

633 c, d = y 

634 r1 = mpi_mul(a,c) 

635 r2 = mpi_mul(b,d) 

636 re = mpi_sub(r1,r2,prec) 

637 i1 = mpi_mul(a,d) 

638 i2 = mpi_mul(b,c) 

639 im = mpi_add(i1,i2,prec) 

640 return re, im 

641 

642def mpci_div(x, y, prec): 

643 # TODO: optimize for real/imag cases 

644 a, b = x 

645 c, d = y 

646 wp = prec+20 

647 m1 = mpi_square(c) 

648 m2 = mpi_square(d) 

649 m = mpi_add(m1,m2,wp) 

650 re = mpi_add(mpi_mul(a,c), mpi_mul(b,d), wp) 

651 im = mpi_sub(mpi_mul(b,c), mpi_mul(a,d), wp) 

652 re = mpi_div(re, m, prec) 

653 im = mpi_div(im, m, prec) 

654 return re, im 

655 

656def mpci_exp(x, prec): 

657 a, b = x 

658 wp = prec+20 

659 r = mpi_exp(a, wp) 

660 c, s = mpi_cos_sin(b, wp) 

661 a = mpi_mul(r, c, prec) 

662 b = mpi_mul(r, s, prec) 

663 return a, b 

664 

665def mpi_shift(x, n): 

666 a, b = x 

667 return mpf_shift(a,n), mpf_shift(b,n) 

668 

669def mpi_cosh_sinh(x, prec): 

670 # TODO: accuracy for small x 

671 wp = prec+20 

672 e1 = mpi_exp(x, wp) 

673 e2 = mpi_div(mpi_one, e1, wp) 

674 c = mpi_add(e1, e2, prec) 

675 s = mpi_sub(e1, e2, prec) 

676 c = mpi_shift(c, -1) 

677 s = mpi_shift(s, -1) 

678 return c, s 

679 

680def mpci_cos(x, prec): 

681 a, b = x 

682 wp = prec+10 

683 c, s = mpi_cos_sin(a, wp) 

684 ch, sh = mpi_cosh_sinh(b, wp) 

685 re = mpi_mul(c, ch, prec) 

686 im = mpi_mul(s, sh, prec) 

687 return re, mpi_neg(im) 

688 

689def mpci_sin(x, prec): 

690 a, b = x 

691 wp = prec+10 

692 c, s = mpi_cos_sin(a, wp) 

693 ch, sh = mpi_cosh_sinh(b, wp) 

694 re = mpi_mul(s, ch, prec) 

695 im = mpi_mul(c, sh, prec) 

696 return re, im 

697 

698def mpci_abs(x, prec): 

699 a, b = x 

700 if a == mpi_zero: 

701 return mpi_abs(b) 

702 if b == mpi_zero: 

703 return mpi_abs(a) 

704 # Important: nonnegative 

705 a = mpi_square(a) 

706 b = mpi_square(b) 

707 t = mpi_add(a, b, prec+20) 

708 return mpi_sqrt(t, prec) 

709 

710def mpi_atan2(y, x, prec): 

711 ya, yb = y 

712 xa, xb = x 

713 # Constrained to the real line 

714 if ya == yb == fzero: