# Ticket #15944: real_and_complex_numbers.rst

File real_and_complex_numbers.rst, 32.3 KB (added by , 4 years ago) |
---|

Line | |
---|---|

1 | .. escape-backslashes |

2 | |

3 | ========================================== |

4 | Tutorial: Real and complex numbers in Sage |

5 | ========================================== |

6 | |

7 | .. MODULEAUTHOR:: Thierry Monteil <sage!lma.metelu.net> |

8 | |

9 | .. contents:: |

10 | :depth: 2 |

11 | |

12 | |

13 | Even with unlimitted memory, computers can not represent all real or complex |

14 | numbers, since those form uncountable sets. Sage proposes various classes to |

15 | approach the set of real numbers. They are considered as ``Fields`` or ``Rings`` |

16 | though they do not necessarily satisfy their respective axiomatic. |

17 | |

18 | Those various classes have very different features, hence it is important to |

19 | know to which one a given number belongs to understand its behaviour. |

20 | |

21 | You should know a bit about coercion before starting this tutorial, so it is |

22 | advised to read the ``parent_element_coercion.en.rst`` worksheet first. |

23 | |

24 | |

25 | The big zoo |

26 | =========== |

27 | |

28 | Here is a rough qualitative summary of the main representations that will appear |

29 | in this tutorial: |

30 | |

31 | =========== =========================== =========== =============================== =========== =========== =========== =========== ======= ============================= |

32 | Type Name Short Complex equivalent Short Reliable Fast Expressive Exact Precision |

33 | =========== =========================== =========== =============================== =========== =========== =========== =========== ======= ============================= |

34 | Numerical ``RealDoubleField()`` ``RDF`` ``ComplexDoubleField()`` ``CDF`` ``+++`` ``+++++`` ``+`` No bounded by 53 bits |

35 | ``RealField(prec)`` ``RR`` ``ComplexField(prec)`` ``CC`` ``+++`` ``+++`` ``++`` No finite unbounded non-adaptive |

36 | ``RealIntervalField(prec)`` ``RIF`` ``ComplexIntervalField(prec)`` ``CIF`` ``+++++`` ``+++`` ``++`` No finite unbounded non-adaptive |

37 | ``RealBallField(prec)`` ``RBF`` ``ComplexBallField(prec)`` ``CBF`` ``+++++`` ``+++`` ``++`` No finite unbounded non-adaptive |

38 | |

39 | Algebraic ``AlgebraicRealField()`` ``AA`` ``AlgebraicField()`` ``QQbar`` ``+++++`` ``+`` ``++++`` Yes infinite |

40 | ``NumberField(poly)`` none itself none ``+++++`` ``++`` ``++++`` Yes infinite |

41 | ``QuadraticField(n)`` none itself none ``+++++`` ``+++`` ``+++`` Yes infinite |

42 | ``RationalField()`` ``QQ`` ``QuadraticField(-1)`` none ``+++++`` ``++++`` ``+++`` Yes infinite |

43 | |

44 | Symbolic ``SymbolicRing()`` ``SR`` itself ``SR`` ``+`` ``+`` ``+++++`` No depends on construction |

45 | =========== =========================== =========== =============================== =========== =========== =========== =========== ======= ============================= |

46 | |

47 | Numerical representations are fast but inexact, they require some special care |

48 | since the precision may become fuzzy along a computation. |

49 | |

50 | Algebraic representations are exact and reliable, but slower and able to |

51 | represent only a class of numbers. |

52 | |

53 | Symbolic representations are very versatile (they can represent numbers such as |

54 | `\sqrt{log(\pi)+sin(42)}`), but slow and sometimes not reliable. |

55 | |

56 | Let us inspect each kind of representation further. |

57 | |

58 | .. topic:: TODO |

59 | |

60 | Find a better word for "kind" that is not "type", "group", "category", |

61 | "class" (perhaps "taxon", "family", "sort" ?). |

62 | |

63 | |

64 | Numerical representations |

65 | ========================= |

66 | |

67 | Floating-point numbers |

68 | ---------------------- |

69 | |

70 | Floating-point numbers are real numbers of the form `(-1)^s 0.d_0 d_1 d_2 \dots |

71 | d_{p-1} 2^{e}` where: |

72 | |

73 | - `s \in \{0,1\}` is the *sign*, |

74 | - `d_i \in \{0,1\}` (with `d_0=1`) are the *digits*, |

75 | - the integer `d_0 d_1 d_2 \dots d_{p-1}` is the *mantissa*, |

76 | - `p` is the *precision*, |

77 | - `e` is the *exponent*. |

78 | |

79 | Floating-point numbers of given precision do not form a field in the |

80 | mathematical sense since they are not stable under the classical operations. So |

81 | when adding or multiplying two floating-point numbers, Sage has to round the |

82 | result to return a floating-point number. In particular, the results are |

83 | approximate and the error can grow along a computation. |

84 | |

85 | Some first warnings about floating-point arithmetics can be found on the guide: |

86 | http://floating-point-gui.de/ |

87 | |

88 | Error propagation is sometimes unavoidable, but it is possible to know a bound |

89 | on it by using the following secure fields: |

90 | |

91 | There exists various floating-point implementations of numbers in Sage: |

92 | |

93 | |

94 | Real Double Field |

95 | ^^^^^^^^^^^^^^^^^ |

96 | |

97 | Elements of ``RDF = RealDoubleField()`` are the floating-point numbers in double |

98 | precision from the processor (those you use when you code in C, similar to |

99 | Python `float`), see :mod:`sage.rings.real_double`. The computations using those |

100 | numbers are very fast and most optimized libraries that work with floating-point |

101 | numbers use this representation, so you will benefit from them if the number |

102 | involved in your floating-point constructions (matrices, polynomials,...) belong |

103 | to ``RDF``. The mantissa has 53 bits of precision, the exponent is encoded on 11 |

104 | bits from -1023 to 1024. During the computations, the rounding is done toward |

105 | the closest floating-point number. |

106 | |

107 | |

108 | Real Fields |

109 | ^^^^^^^^^^^ |

110 | |

111 | Elements of ``RealField(prec)`` are the floating-point numbers with ``prec`` |

112 | bits of precision, "emulated" by the MPFR library. By default, the mantissa has |

113 | `53` bits of precision (and ``RealField(prec=53)`` can be shortened as ``RR``), |

114 | the exponent is encoded on 32 bits (so that you will unlikely overflow), see |

115 | :mod:`sage.rings.real_mpfr`. The computations are slower than in ``RDF`` and |

116 | most algorithms on objects whose base ring is ``RealField(prec)`` will be pretty |

117 | naive. However, it is possible to extend the precision as well as the type of |

118 | rounding (towards `+\infty`, towards `-\infty`, towards `0`), for example |

119 | ``RealField(100, rnd='RNDZ')`` is the "field" of real floating-point numbers |

120 | with 100 bits of precision with rounding towards zero. |

121 | |

122 | |

123 | .. note:: |

124 | |

125 | While sometimes needed, increasing the precision also increases the |

126 | computation time. |

127 | |

128 | |

129 | .. warning:: |

130 | |

131 | Note that despite its generic name, ``RealField(prec)`` or ``RR`` is *not* |

132 | an abstraction of the field of real numbers, as can be seen by the following |

133 | section. |

134 | |

135 | Unreal floats |

136 | ^^^^^^^^^^^^^ |

137 | |

138 | Besides not representing all real numbers, ``RDF`` and ``RR`` contain three |

139 | additional unreal elements that allow the ``(+,-,*,/)`` operations to be always |

140 | defined (``NaN`` stands for *Not a Number*): |

141 | |

142 | :: |

143 | |

144 | sage: RDF(0) / RDF(0) |

145 | NaN |

146 | |

147 | sage: RDF(1) / RDF(0) |

148 | +infinity |

149 | |

150 | sage: RR(-1) / RR(0) |

151 | -infinity |

152 | |

153 | sage: RDF(infinity) - RDF(infinity) |

154 | NaN |

155 | |

156 | sage: NaN in RR and Infinity in RR and -Infinity in RR |

157 | True |

158 | |

159 | |

160 | Back and forth in precision |

161 | ^^^^^^^^^^^^^^^^^^^^^^^^^^^ |

162 | |

163 | Note that it is is general useless to improve the precision during a |

164 | computation. This explains why the coercion is done toward the smallest |

165 | precision: |

166 | |

167 | :: |

168 | |

169 | sage: R = RealField(200) |

170 | sage: S = RealField(100) |

171 | sage: a = R(2).sqrt() |

172 | sage: b = S(pi) |

173 | sage: a+b == R(sqrt(2)+pi) |

174 | True |

175 | sage: a+R(b) == R(sqrt(2)+pi) |

176 | False |

177 | |

178 | |

179 | .. note:: |

180 | |

181 | For each element of those two fields, there exists a |

182 | ``.sign_mantissa_exponent()`` method that returns the corresponding |

183 | triple. |

184 | |

185 | |

186 | Floating-point numbers with a control on the imprecision |

187 | -------------------------------------------------------- |

188 | |

189 | Real Interval Fields |

190 | ^^^^^^^^^^^^^^^^^^^^ |

191 | |

192 | Elements of ``RealIntervalField(prec)`` are pairs of elements of |

193 | ``RealField(prec)``, which should be interpreted as the two endpoints of an |

194 | interval containing the represented real number. ``RealIntervalField(53)`` can be |

195 | shortened as ``RIF``. |

196 | |

197 | Operations defined over the ``RealIntervalField`` implementation must offer the |

198 | *containment guaranty*: for every mathematical function `f : \mathbb{R}^d \to |

199 | \mathbb{R}` that admits an implementation ``f_rif`` over some |

200 | ``RealIntervalField(prec)``, we have the property that for every element ``x = |

201 | [(a_0, b_0), (a_1,b_1), ..., (a_(d-1),b_(d-1))]`` in |

202 | ``RealIntervalField(prec)^d``, if ``f_rif(x) = [a', b']``, then `f([a_0, b_0] |

203 | \times [a_1,b_1] \times \dots \times [a_(d-1),b_(d-1)])\subseteq [a', b']`. In |

204 | particular, the left endpoint of the interval is rounded towards `-\infty` and |

205 | the right endpoint is rounded towards `+\infty` |

206 | |

207 | |

208 | |

209 | Real Ball Fields |

210 | ^^^^^^^^^^^^^^^^ |

211 | |

212 | Elements of ``RealBallField(prec)`` are pairs ``(c,r)``, where ``c`` is a |

213 | floating-point number with ``prec`` bits of precision, and where ``r`` is a |

214 | floating-point number with 30 bits of precision. The high precision ``c`` and |

215 | the low precision number ``r`` should be understood as the center and the radius |

216 | (respectively) of a ball containing the represented real number. |

217 | ``RealBallField(53)`` can be shortened as ``RBF``. |

218 | |

219 | Operations defined over the ``RealBallField`` implementation must offer a |

220 | similar containment guaranty as for ``RealIntervalField`` (replace intervals |

221 | with balls). |

222 | |

223 | Comparison between both representations |

224 | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |

225 | |

226 | The main advantage of the ball implementation over interval arithmetics is that, |

227 | when the precision is high, the error is still encoded with a low precision |

228 | (since only the order of magnitude matters for the error bounds), so that the |

229 | loss of speed compared to non-guaranteed floating-point arithmetics with the |

230 | same precision is about `1+\varepsilon` (and not `2` as in the case of interval |

231 | arithmetics, where two precise numbers have to be computed). |

232 | |

233 | Also, it is closer to the mathematical `\varepsilon-\delta` definitions and |

234 | distance-based error estimations. |

235 | |

236 | However, for large intervals and for evaluating monotonic functions, |

237 | interval arithmetics is usually more precise and easier to implement (just |

238 | round the images of the endpoints towards the right direction). It is also |

239 | well-adapted for implementing dichotomy-based precision refinements. |

240 | |

241 | |

242 | Which numerical representation to choose |

243 | ---------------------------------------- |

244 | |

245 | Sage interfaces various libraries written in C, C++, Fortran, etc. which |

246 | use the processor double precision floats. Hence, the use of ``RDF`` has |

247 | the advantage of being faster and allowing the use of those libraries. |

248 | |

249 | If you use other floating-point representations, Sage will mostly use |

250 | generic algorithms, which are usually slower and numerically less stable. |

251 | However, those representations allow you to express real numbers with a |

252 | higher precision. |

253 | |

254 | Regarding speed, ``RealBallField(prec)`` and ``RealField(prec)`` are |

255 | similar, while ``RealIntervalField(prec)`` is twice slower. However, the |

256 | ``RealBallField(prec)`` is a new feature of Sage, hence some methods might |

257 | raise a ``NotImplementedError``. |

258 | |

259 | .. topic:: Summary |

260 | |

261 | So, if you need speed and the power of Sage libraries: use ``RDF``. |

262 | |

263 | If you need more precision, start with ``RealBallField(prec)``, and |

264 | try ``RealIntervalField(prec)`` in case the imprecision becomes too |

265 | large (alternatively: increase the precision or try to think about the |

266 | order and the expanding nature of your operations). Use |

267 | ``RealField(prec)`` if the methods you wanted to use with |

268 | ``RealBallField(prec)`` are not implemented. |

269 | |

270 | |

271 | Exercises |

272 | --------- |

273 | |

274 | Binomial |

275 | ^^^^^^^^ |

276 | |

277 | From the three following values: |

278 | |

279 | :: |

280 | |

281 | sage: a = 1.0 |

282 | sage: b = 10.0^4 |

283 | sage: c = 1.0 |

284 | |

285 | solve the `ax^2+bx+c = 0` equation by the method you learned in highschool |

286 | involving discriminant, and compare the sum and product of the solution |

287 | you found to their theoretical values `-b/a` and `c/a`. |

288 | |

289 | :: |

290 | |

291 | sage: # edit here |

292 | |

293 | Same question by letting Sage solve this equation for you. |

294 | |

295 | :: |

296 | |

297 | sage: # edit here |

298 | |

299 | Even on such a simple computation, you have been victim of a phenomenon of |

300 | *error amplification*. Another striking example is the following: |

301 | |

302 | :: |

303 | |

304 | sage: a = 10000.0 |

305 | sage: b = -9999.5 |

306 | sage: c = 0.1 |

307 | sage: a+c+b |

308 | 0.600000000000364 |

309 | sage: a+b+c |

310 | 0.600000000000000 |

311 | |

312 | Note that, by construction and following the IEEE 754 specification, the |

313 | ``+`` (and ``*``) law is commutative in ``RDF`` and ``RR``:: |

314 | |

315 | sage: all([a+b==b+a, a+c==c+a, b+c==c+b, (a+b)+c==c+(a+b), (a+c)+b==b+(a+c), (b+c)+a==a+(b+c)]) |

316 | True |

317 | |

318 | Actually the associativity is not satisfied here: |

319 | |

320 | :: |

321 | |

322 | sage: (a+b)+c == a+(b+c) |

323 | False |

324 | |

325 | Try to find, by random sampling, 3 elements `a,b,c` in ``RDF`` such that |

326 | ``(a+b)+c != a+(b+c)``: |

327 | |

328 | :: |

329 | |

330 | sage: # edit here |

331 | |

332 | Why is it so hard ? What makes the previous example work ? |

333 | |

334 | |

335 | Interval arithmetics |

336 | ^^^^^^^^^^^^^^^^^^^^ |

337 | |

338 | "Prove" using ``RealIntervalField`` that `cos(2\pi)=1` with a high precision |

339 | (say `10^{-100}`). |

340 | |

341 | :: |

342 | |

343 | sage: # edit here |

344 | |

345 | Explain the following behaviour: |

346 | |

347 | :: |

348 | |

349 | sage: sqrt(0.1)^2 == 0.1 |

350 | True |

351 | sage: sqrt(RIF(0.1))^2 == RIF(0.1) |

352 | False |

353 | |

354 | Provide a relation between ``sqrt(RIF(0.1))^2`` and ``RIF(0.1)`` |

355 | |

356 | :: |

357 | |

358 | sage: # edit here |

359 | |

360 | |

361 | Note that the ordering of the computations has an influence on the result, |

362 | but also on its precision. With the previous example, we have: |

363 | |

364 | :: |

365 | |

366 | sage: a = RIF(10000.0) |

367 | sage: b = RIF(9999.5) |

368 | sage: c = RIF(0.1) |

369 | sage: a+c-b |

370 | 0.60000000000? |

371 | sage: a-b+c |

372 | 0.6000000000000000? |

373 | |

374 | sage: (a+c-b).endpoints() |

375 | (0.599999999998544, 0.600000000000364) |

376 | sage: (a+c-b).absolute_diameter() |

377 | 1.81898940354586e-12 |

378 | |

379 | sage: (a-b+c).endpoints() |

380 | (0.599999999999999, 0.600000000000001) |

381 | sage: (a-b+c).absolute_diameter() |

382 | 1.11022302462516e-16 |

383 | |

384 | |

385 | Before evaluating them, could you predict and explain the behaviour of the |

386 | following commands ? |

387 | |

388 | :: |

389 | |

390 | sage: # 1/2 in RIF |

391 | |

392 | sage: # 1/3 in RIF |

393 | |

394 | *Hint*: you can look at the documentation of the ``.__contains__()`` method of |

395 | ``RIF``. |

396 | |

397 | |

398 | A Sage developer wants to write the real function `f` that maps `0.01` to |

399 | `1` and that is the identity elsewhere, and suggests the following |

400 | implementation: |

401 | |

402 | :: |

403 | |

404 | sage: def f_unstable(x): |

405 | ....: if x == 0.01: |

406 | ....: return 1 |

407 | ....: else: |

408 | ....: return x |

409 | |

410 | Prove that this implementation does not offers the containment guaranty |

411 | when used on the ``RealIntervalField``. |

412 | |

413 | Could you help this lost developer to find a correct implementation so |

414 | that it could be included in Sage source code ? |

415 | |

416 | :: |

417 | |

418 | sage: def f_correct(x): |

419 | sage: # edit here |

420 | |

421 | |

422 | Divergent series ? |

423 | ^^^^^^^^^^^^^^^^^^ |

424 | |

425 | The series `u_n := 1 + 1/2 + 1/3 + 1/4 + \dots + 1/n = \sum_{k=1}^{n} 1/k` is |

426 | known to diverge on `\mathbb{R}`. |

427 | |

428 | Explore its convergence on some `Real Fields` (start with small precision). |

429 | Explore the difference with the (mathematically equal) series `v_n := 1/n + |

430 | 1/(n-1) + ... + 1/2 + 1`. |

431 | |

432 | Explain, formulate conjectures and try to prove them. |

433 | |

434 | |

435 | A picture |

436 | ^^^^^^^^^ |

437 | |

438 | Draw (by hand) a picture of the elements of ``RealField(3)``. |

439 | |

440 | |

441 | Convergence of a Markov Chain |

442 | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |

443 | |

444 | Let $M$ be the rational bistochastic matrix defined by:: |

445 | |

446 | sage: M = matrix(QQ, [[1/6, 1/2, 1/3], [1/2, 1/4, 1/4], [1/3, 1/4, 5/12]]) |

447 | sage: M |

448 | [ 1/6 1/2 1/3] |

449 | [ 1/2 1/4 1/4] |

450 | [ 1/3 1/4 5/12] |

451 | |

452 | Could you guess the limit of $M^n$ when $n$ goes to $+\infty$ ? |

453 | |

454 | :: |

455 | |

456 | sage: # edit here |

457 | |

458 | Could you prove that this limit exists and give its exact value ? |

459 | |

460 | :: |

461 | |

462 | sage: # edit here |

463 | |

464 | |

465 | Repulsive and attractive fixed points |

466 | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |

467 | |

468 | Define the maps `f:x\to 4x-1` and `g:x\to x/3+1`. |

469 | |

470 | :: |

471 | |

472 | sage: # edit here |

473 | |

474 | Let `u_n` and `v_n` be the sequences defined recursively by |

475 | `u_0 = RDF(1/3)`, `u_{n+1} = f(u_n)` and |

476 | `v_0 = RDF(3/2)`, `v_{n+1} = g(v_n)`. |

477 | |

478 | Compute the first 100 values of the sequences `u_n` and `v_n`. |

479 | |

480 | :: |

481 | |

482 | sage: # edit here |

483 | |

484 | For which reason you should not worry to find the fixed point of the function `h:x\to sin(x + 1/2)` ? Compute the value of this fixed point with 50 bits of precision. |

485 | |

486 | :: |

487 | |

488 | sage: # edit here |

489 | |

490 | |

491 | Decimals of pi |

492 | ^^^^^^^^^^^^^^ |

493 | |

494 | What is the 150th decimal of `\pi` ? |

495 | |

496 | :: |

497 | |

498 | sage: # edit here |

499 | |

500 | |

501 | Money, money, money |

502 | ^^^^^^^^^^^^^^^^^^^ |

503 | |

504 | :: |

505 | |

506 | sage: one_cent = 0.01 |

507 | sage: my_money = 0 |

508 | sage: for i in range(100): |

509 | sage: my_money += one_cent |

510 | sage: my_money > 1 |

511 | |

512 | What is the origin of this behaviour ? |

513 | |

514 | How to fix the issue in Sage (respecively un pure Python, compare both) ? |

515 | |

516 | :: |

517 | |

518 | sage: # edit here |

519 | |

520 | Which risk does computer science take if finance and high-frequency |

521 | trading stays on power for a while ? |

522 | |

523 | Provide the first 30 bits of the binary expansion of `0.01`. |

524 | |

525 | :: |

526 | |

527 | sage: # edit here |

528 | |

529 | |

530 | Huge matrix computation |

531 | ^^^^^^^^^^^^^^^^^^^^^^^ |

532 | |

533 | Construct a dense matrix of size ``500*500`` with real coefficients |

534 | randomly distributed in `[0,1]`. |

535 | |

536 | :: |

537 | |

538 | sage: # edit here |

539 | |

540 | Take the 100th power of this matrix in less than one second. |

541 | |

542 | :: |

543 | |

544 | sage: # edit here |

545 | |

546 | Comment. |

547 | |

548 | |

549 | Exponents |

550 | ^^^^^^^^^ |

551 | |

552 | :: |

553 | |

554 | sage: RDF(1/2^10000) == 0 |

555 | True |

556 | sage: RR(1/2^10000) == 0 |

557 | False |

558 | |

559 | Comment. |

560 | |

561 | A strange impression |

562 | ^^^^^^^^^^^^^^^^^^^^ |

563 | :: |

564 | |

565 | sage: RR(1/pi) |

566 | 0.318309886183791 |

567 | sage: RDF(1/pi) |

568 | 0.3183098861837907 |

569 | |

570 | |

571 | Before evaluating the next command, could you tell what will be the answer of |

572 | the following command ? |

573 | |

574 | :: |

575 | |

576 | sage: # RR(1/pi) == RDF(1/pi) |

577 | |

578 | What is the cause of the behaviour of the first two commands ? |

579 | |

580 | Rounding |

581 | ^^^^^^^^ |

582 | |

583 | .. https://ask.sagemath.org/question/32169/why-5000000000000000-is-not-equal-to-5/ |

584 | |

585 | Let us define:: |

586 | |

587 | sage: n = 5.0000000000000001 |

588 | |

589 | Explain each of the following behaviours:: |

590 | |

591 | sage: n |

592 | 5.000000000000000 |

593 | |

594 | :: |

595 | |

596 | sage: n.str(truncate=False) |

597 | '5.000000000000000111' |

598 | |

599 | :: |

600 | |

601 | sage: 5.000000000000000111 == 5.0000000000000001 |

602 | |

603 | |

604 | |

605 | Algebraic (exact) representations |

606 | ================================= |

607 | |

608 | Some subfields of the real numbers can be represented in an exact way in |

609 | Sage. They are pretty safe (no rounding problem, well defined |

610 | semantics) through sometimes slow. |

611 | |

612 | As a general rule, the smallest the field is (for the inclusion), the |

613 | fastest it is, which means that, in principle, working within rational |

614 | field is faster than within quadratic fields, which is faster than within |

615 | cyclotomic fields, which is faster than within number fields, which is |

616 | faster that within the universal cyclotomic field, that is faster than |

617 | within the whole algebraic field. |

618 | |

619 | Here is a short overview: |

620 | |

621 | Rational numbers |

622 | ---------------- |

623 | |

624 | ``QQ = RationalField()`` stands for the field of rational numbers. |

625 | Computations are exact and only limited by the available memory. You |

626 | should note that the numerators/denominators could become very large along |

627 | a computation, which might slow down this latter. |

628 | |

629 | :: |

630 | |

631 | sage: (123/321)^10000 |

632 | |

633 | It is possible to "approximate" a floating-point real number by rational |

634 | numbers with the methods ``nearby_rational()``, ``simplest_rational()``, |

635 | ``exact_rational()``. |

636 | |

637 | :: |

638 | |

639 | sage: # edit here |

640 | |

641 | If ``x`` belongs to ``RR``, what is the result of ``QQ(x)`` ? Is it a |

642 | coercion ? |

643 | |

644 | See the page |

645 | http://doc.sagemath.org/html/en/reference/rings_standard/sage/rings/rational_field.html |

646 | |

647 | |

648 | Algebraic number field |

649 | ---------------------- |

650 | |

651 | ``AA = AlgebraicRealField()`` stands for te field of real algebraic numbers. |

652 | |

653 | :: |

654 | |

655 | sage: a = AA(2).sqrt() + AA(7).sqrt() |

656 | sage: a |

657 | 4.059964873437685? |

658 | |

659 | Here, the number ``a`` is *printed* as its approximated numerical value, |

660 | but the question mark at the end indicates that it is not internally a |

661 | floating-point number. |

662 | |

663 | Basically, an algebraic number is represented as the root of an integer |

664 | polynomial, and a small interval surrounding the appropriate root: |

665 | |

666 | :: |

667 | |

668 | sage: a.minpoly() |

669 | x^4 - 18*x^2 + 25 |

670 | |

671 | sage: a.n(digits=50) |

672 | 4.0599648734376856393033044778489585042799310584594 |

673 | sage: a.interval_diameter(1e-40) |

674 | 4.0599648734376856393033044778489585042799310584593982535450141971918013016924? |

675 | |

676 | Computing such a representation is costly, hence Sage is lazy and computes |

677 | it only when needed: |

678 | |

679 | :: |

680 | |

681 | sage: b = QQbar(2).sqrt() |

682 | sage: c = a^2 |

683 | sage: c |

684 | 2.000000000000000? |

685 | |

686 | Here, Sage did not compute the minimal polynomial nor an enclosing |

687 | interval of ``c``, it just knows ``c`` as "the square of ``b``", in |

688 | particular, it does not know yet that ``c`` is the integer ``2``. However, |

689 | such an "exactification" is computed when required (and kept internally): |

690 | |

691 | :: |

692 | |

693 | sage: c.is_integer() |

694 | True |

695 | sage: b |

696 | 2 |

697 | |

698 | You can force the exactification with the ``exactify`` method. |

699 | |

700 | For more details, see the page |

701 | http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/qqbar.html |

702 | |

703 | |

704 | Number fields |

705 | ------------- |

706 | |

707 | See the pages: |

708 | |

709 | - http://doc.sagemath.org/html/en/reference/number_fields/index.html |

710 | - http://doc.sagemath.org/html/en/constructions/number_fields.html |

711 | - http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/number_field/number_field.html |

712 | |

713 | .. warning:: |

714 | |

715 | You should be aware that number field are purely algebraic objects that are |

716 | *not* embedded into the complex plane by default. Do to this, you have to |

717 | define the embedding you want to work with (if you need one). |

718 | |

719 | Note that the quadratic number fields have a dedicated faster |

720 | implementation. |

721 | |

722 | What is the difference between those two objects ? |

723 | |

724 | :: |

725 | |

726 | sage: q = QuadraticField(-1,'I') |

727 | sage: qq = QQ[I] |

728 | |

729 | |

730 | Cyclotomic fields (and the universal cyclotomic field) |

731 | ------------------------------------------------------ |

732 | |

733 | See the pages: |

734 | |

735 | - http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/number_field/number_field.html#sage.rings.number_field.number_field.CyclotomicFieldFactory |

736 | - http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/universal_cyclotomic_field.html |

737 | |

738 | |

739 | Symbolic representations |

740 | ======================== |

741 | |

742 | ``SR = sage.symbolic.ring.SymbolicRing()`` is the set of symbolic expressions. |

743 | |

744 | ``SR`` is pretty fuzzy but quite handy (here expressivity is prefered over |

745 | standardization and the number of real numbers representable within ``SR`` |

746 | is constantly increasing). This "ring" contains most mathematical |

747 | constants such as `e`, `\pi`, `cos(12)`, `\sqrt{2}`, `\gamma` (Euler |

748 | constant), ... and allows some computations such as: |

749 | |

750 | :: |

751 | |

752 | sage: e^(I*pi) |

753 | -1 |

754 | |

755 | .. warning:: |

756 | |

757 | ``SR`` also contains symbolic expressions that are not real or complex numbers: |

758 | |

759 | :: |

760 | |

761 | sage: e = cos(x)^2 |

762 | sage: e.derivative() |

763 | -2*cos(x)*sin(x) |

764 | sage: _ in SR |

765 | True |

766 | |

767 | ``SR`` represents numbers as symbolic expressions, whose length (as a |

768 | formula) can grow along a computation. This phenomenon can slow down the |

769 | computations. If our aim is to end with a numerical value, you should use |

770 | a floating-point representation, as they keep a constant length. |

771 | |

772 | :: |

773 | |

774 | sage: a = exp(2) |

775 | sage: f = lambda x : sqrt(x+1)^x |

776 | sage: for i in range(10): |

777 | ....: a = f(a) |

778 | sage: a |

779 | |

780 | Of course, using floating-point arithmetics will not necessarily solve all |

781 | your problems since you will have to ensure the numerical stability of |

782 | such iterated computation (see the sections above). |

783 | |

784 | |

785 | ``SR`` always answer, even when it does not know the answer, for example: |

786 | |

787 | :: |

788 | |

789 | sage: # missing example |

790 | |

791 | .. sage: alpha = 2*pi/7 |

792 | .. sage: cos(alpha)^(1/3)-(-cos(2*alpha))^(1/3)-(-cos(4*alpha))^(1/3) + ((-5+3*7^(1/3))/2)^(1/3) |

793 | .. sage: cos(alpha)^(1/3)-(-cos(2*alpha))^(1/3)-(-cos(4*alpha))^(1/3) == - ((-5+3*7^(1/3))/2)^(1/3) |

794 | .. sage: %time bool(_) |

795 | .. CPU times: user 12min 49s, sys: 33.6 s, total: 13min 22s |

796 | .. Wall time: 14min 17s |

797 | .. True |

798 | .. https://cs.uwaterloo.ca/journals/JIS/VOL12/Witula/witula17.pdf |

799 | |

800 | |

801 | Hence, computations involving the symbolic ring should not be considered |

802 | as reliable. In particular, ``SR`` is a dead-end for the coercion among |

803 | numbers (and can add pears and apples formally): |

804 | |

805 | :: |

806 | |

807 | sage: 0.1 + pi + QQbar(2).sqrt() |

808 | pi + 1.51421356237310 |

809 | sage: _.parent() |

810 | Symbolic Ring |

811 | |

812 | Here are some issues that can appear with the symbolic ring, when you can |

813 | represent your numbers in a different way, we strongly encourage you to do |

814 | so: |

815 | |

816 | - https://ask.sagemath.org/questions/tags:symbolic_issue/ |

817 | - https://trac.sagemath.org/wiki/symbolics |

818 | - https://ask.sagemath.org/question/2469/quadratic-number-fields |

819 | - https://ask.sagemath.org/question/2561/is-it-a-bug |

820 | - https://ask.sagemath.org/question/2518/get-variants-of-complex-cube-root |

821 | - https://ask.sagemath.org/question/2734/computing-the-digist-of-pi |

822 | |

823 | |

824 | |

825 | Conversions, bridges, exactification |

826 | ==================================== |

827 | |

828 | .. topic:: TODO |

829 | |

830 | Should this be a separate section, or should it appear in the |

831 | rational/algebraic sections ? |

832 | |

833 | Along a computation, coercion goes towards the less precise numbers: |

834 | |

835 | :: |

836 | |

837 | sage: a = 1/10 ; a |

838 | 1/10 |

839 | sage: a.parent() |

840 | Rational Field |

841 | sage: b = QQbar(2).sqrt() ; b |

842 | 1.414213562373095? |

843 | sage: b.parent() |

844 | Algebraic Field |

845 | sage: c = a+b ; c |

846 | 1.514213562373096? |

847 | sage: c.parent() |

848 | Algebraic Field |

849 | sage: d = RDF(0.1) ; d |

850 | 0.1 |

851 | sage: d.parent() |

852 | Real Double Field |

853 | sage: e = c+d ; e |

854 | 1.614213562373095 |

855 | sage: e.parent() |

856 | Complex Double Field |

857 | sage: f = pi ; f |

858 | pi |

859 | sage: f.parent() |

860 | Symbolic Ring |

861 | sage: g = e+f ; g |

862 | pi + 1.614213562373095 |

863 | sage: g.parent() |

864 | Symbolic Ring |

865 | |

866 | |

867 | But you might want to swim upstream. Given a floating point approximation of a |

868 | real number, there are various ways to recover an exact representation it *may* |

869 | come from. |

870 | |

871 | |

872 | From numerical to rational |

873 | -------------------------- |

874 | |

875 | An element ``x`` of ``RealField(prec)`` can be transformed into a rational in |

876 | different ways: |

877 | |

878 | |

879 | - ``x.exact_rational()`` returns the exact rational representation of the |

880 | floating-point number ``x``. In particular its denominator is a power of |

881 | 2. It contains the same information as ``.sign_mantissa_exponent``: |

882 | |

883 | :: |

884 | |

885 | sage: a = RR(pi) ; a |

886 | 3.14159265358979 |

887 | sage: a.sign_mantissa_exponent() |

888 | (1, 7074237752028440, -51) |

889 | sage: s,m,e = _ |

890 | sage: s*m*2^(e) |

891 | 884279719003555/281474976710656 |

892 | sage: b = a.exact_rational() |

893 | 884279719003555/281474976710656 |

894 | |

895 | |

896 | - ``x.simplest_rational()`` returns the simplest rational (i.e. with |

897 | smallest denominator) that will become equal to ``x`` when converted |

898 | into the parent of ``x``: |

899 | |

900 | :: |

901 | |

902 | sage: a = RR(pi) ; a |

903 | 3.14159265358979 |

904 | sage: b = a.simplest_rational() ; b |

905 | 245850922/78256779 |

906 | sage: a == b |

907 | True |

908 | sage: RR(b) |

909 | 3.14159265358979 |

910 | |

911 | |

912 | |

913 | - ``nearby_rational(max_error=e)`` finds the simplest rational that is at |

914 | distance at most ``e`` of ``x``: |

915 | |

916 | :: |

917 | |

918 | sage: RR(pi).nearby_rational(max_error=0.01) |

919 | 22/7 |

920 | |

921 | - ``nearby_rational(max_denominator=d)`` finds the closest rational to |

922 | ``x`` with denominator at most ``d``: |

923 | |

924 | :: |

925 | |

926 | sage: RR(pi).nearby_rational(max_denominator=120) |

927 | 355/113 |

928 | |

929 | |

930 | From numerical to algebraic |

931 | --------------------------- |

932 | |

933 | If ``x`` is an element of ``RealField(prec)``, |

934 | ``x.algebraic_dependency(d)`` returns an integer polynomial ``P`` of |

935 | degree at most ``d`` such that ``P(x)`` is almost zero: |

936 | |

937 | :: |

938 | |

939 | sage: a = 3.0.nth_root(3) |

940 | sage: a += 5*a.ulp() |

941 | sage: a.algebraic_dependency(3) |

942 | x^3 - 3 |

943 | |

944 | |

945 | From algebraic field to (smaller, faster) embedded number fields |

946 | ---------------------------------------------------------------- |

947 | |

948 | Algebraic numbers have an ``as_number_field_element`` method: |

949 | |

950 | :: |

951 | |

952 | sage: a = QQbar(2).sqrt() + QQbar(3).nth_root(3) |

953 | sage: a.as_number_field_element() |

954 | |

955 | (Number Field in a with defining polynomial y^6 - 6*y^4 + 6*y^3 + 12*y^2 + 36*y + 1, |

956 | 96/755*a^5 - 54/755*a^4 - 128/151*a^3 + 936/755*a^2 + 1003/755*a + 2184/755, |

957 | Ring morphism: |

958 | From: Number Field in a with defining polynomial y^6 - 6*y^4 + 6*y^3 + 12*y^2 + 36*y + 1 |

959 | To: Algebraic Real Field |

960 | Defn: a |--> -0.02803600793431334?) |

961 | |

962 | sage: field, number, morphism = _ |

963 | sage: morphism(number) == a |

964 | True |

965 | |

966 | For more than one algebraic number, the previous method generalizes with |

967 | the ``number_field_elements_from_algebraics`` function: |

968 | |

969 | :: |

970 | |

971 | sage: from sage.rings.qqbar import number_field_elements_from_algebraics |

972 | sage: sqrt2 = AA(2).sqrt() |

973 | sage: sqrt3 = AA(3).sqrt() |

974 | sage: number_field_elements_from_algebraics([sqrt2, sqrt3]) |

975 | (Number Field in a with defining polynomial y^4 - 4*y^2 + 1, |

976 | [-a^3 + 3*a, -a^2 + 2], |

977 | Ring morphism: |

978 | From: Number Field in a with defining polynomial y^4 - 4*y^2 + 1 |

979 | To: Algebraic Real Field |

980 | Defn: a |--> 0.5176380902050415?) |

981 | sage: field, numbers, morphism = _ |

982 | sage: morphism(numbers[0]) == sqrt2 |

983 | True |

984 | |

985 | |

986 | From algebraic to symbolic |

987 | -------------------------- |

988 | |

989 | Algebraic numbers have a ``radical_expression`` method: |

990 | |

991 | :: |

992 | |

993 | sage: x = polygen(QQ,'x') |

994 | sage: P = x^3-x^2-x-1 |

995 | sage: P |

996 | x^3 - x^2 - x - 1 |

997 | sage: P.roots(QQbar)[0][0] |

998 | 1.839286755214161? |

999 | sage: _.radical_expression() |

1000 | (1/9*sqrt(11)*sqrt(3) + 19/27)^(1/3) + 4/9/(1/9*sqrt(11)*sqrt(3) + |

1001 | 19/27)^(1/3) + 1/3 |

1002 | |

1003 | |

1004 | .. note:: |

1005 | |

1006 | Note that the ``radical_expression`` method does not implement the |

1007 | whole Galois theory yet, meaning that some algebraic numbers that can |

1008 | be represented by radicals, might not be modified that way: |

1009 | |

1010 | :: |

1011 | |

1012 | sage: a = QQbar(2).sqrt() + QQbar(3).nth_root(3) |

1013 | sage: a.radical_expression() |

1014 | 2.856463132680504? |

1015 | |

1016 | |

1017 | Exercise: reverse engineering |

1018 | ----------------------------- |

1019 | |

1020 | Suppose that, with a numerical experiment, you have approached an |

1021 | irrational algebraic number by: |

1022 | |

1023 | :: |

1024 | |

1025 | sage: a = 3.14626436994197234 |

1026 | |

1027 | What is this algebraic number likely to be ? |

1028 | |

1029 | :: |

1030 | |

1031 | sage: # edit here |

1032 | |

1033 | |

1034 | |

1035 | Complex numbers |

1036 | =============== |

1037 | |

1038 | Most real number representation have their complex equivalent. |

1039 | |

1040 | - ``RR`` becomes ``CC=ComplexField()`` |

1041 | - ``RDF`` becomes ``CDF=ComplexDoubleField()`` |

1042 | - ``RIF`` becomes ``CIF=ComplexIntervalField()`` |

1043 | - ``RBF`` becomes ``CBF=ComplexBallField()`` |

1044 | - ``QQ`` becomes ``QQ[I]`` (Gauss integers) |

1045 | - ``AA`` becomes ``QQbar`` (the algebraic closure of ``QQ``) |

1046 | - ``SR`` remains ``SR`` (which contains much more than just complex numbers) |

1047 | |

1048 | Note that by default, ``I`` and ``i`` belong to the ``Symbolic Ring``, so you |

1049 | will have to convert it, e.g. ``QQbar(I)``, ``CDF(I)``, ... or obtain it as |

1050 | ``QQbar.gen()``, ``QQbar.gen()``, ... . If you overwrite it, you can recover it |

1051 | with: |

1052 | |

1053 | :: |

1054 | |

1055 | sage: from sage.symbolic.all import i as I |

1056 | |

1057 | |

1058 | Diamonds are not round |

1059 | ---------------------- |

1060 | |

1061 | We saw previously that contracting makes computations safe. Consider the |

1062 | following example: |

1063 | |

1064 | :: |

1065 | |

1066 | sage: a = CBF(1/10+I*4/5) |

1067 | sage: b = CBF((1+I)*2/3) |

1068 | sage: a.diameter() |

1069 | 4.4408921e-16 |

1070 | sage: b.diameter() |

1071 | 4.4408921e-16 |

1072 | |

1073 | :: |

1074 | |

1075 | sage: b.abs() < 1 |

1076 | True |

1077 | |

1078 | :: |

1079 | |

1080 | sage: (a*b).diameter() |

1081 | 1.1768364e-15 |

1082 | |

1083 | Explain this loss of precision. |

1084 | |

1085 | *Hint*: draw a picture ! |

1086 | |

1087 | Complex ball elements are pairs of real ball elements, see |

1088 | http://arblib.org/acb.html |

1089 | |

1090 | |

1091 | |

1092 | Transitional representations (under the hood) |

1093 | ============================================= |

1094 | |

1095 | Some "transitional representations" were not included in the big zoo for |

1096 | simplicity, and because they exist mostly as wrappers for internal stuff. |

1097 | |

1098 | Real literals: between your fingers and Sage |

1099 | --------------------------------------------- |

1100 | |

1101 | When the user writes: |

1102 | |

1103 | :: |

1104 | |

1105 | sage: a = 0.1 |

1106 | |

1107 | this defines a floating-point number though the precision is not |

1108 | explicitely set. |

1109 | |

1110 | Sage defines some default precision (depending on the |

1111 | number of digits that are provided, 53 in the current case): |

1112 | |

1113 | :: |

1114 | |

1115 | sage: type(a) |

1116 | <type 'sage.rings.real_mpfr.RealLiteral'> |

1117 | sage: a.parent() |

1118 | Real Field with 53 bits of precision |

1119 | |

1120 | |

1121 | but it allows casting into higher precision rings, without losing the |

1122 | precision from the default choice: |

1123 | |

1124 | :: |

1125 | |

1126 | sage: RealField(200)(0.1) |

1127 | 0.10000000000000000000000000000000000000000000000000000000000 |

1128 | sage: RealField(200)(RDF(0.1)) |

1129 | 0.10000000000000000555111512312578270211815834045410156250000 |

1130 | |

1131 | |

1132 | |

1133 | Real Lazy Field: between exact and numerical representations |

1134 | ------------------------------------------------------------ |

1135 | |

1136 | Elements of ``RLF = RealLazyField()`` behaves like floating-point numbers, |

1137 | This represen |

1138 | |

1139 | while |

1140 | |

1141 | ``RLF = RealLazyField()`` is lazy, in the sense that it doesn't really do |

1142 | anything but simply sits between exact rings of characteristic 0 and the |

1143 | real numbers. The values are actually computed when they are cast into a |

1144 | field of fixed precision. |

1145 | |

1146 | |

1147 | ``RLF = RealLazyField()`` like ``RealField()`` but can wait that the user fixes the |

1148 | precision to get one. |

1149 | http://www.sagemath.org/doc/reference/rings_numerical/sage/rings/real_lazy.html |

1150 | |

1151 | :: |

1152 | |

1153 | sage: a = RLF(pi) + 2 |

1154 | sage: a |

1155 | 5.141592653589794? |

1156 | sage: RealField(100)(a) |

1157 | 5.1415926535897932384626433833 |

1158 | sage: RealField(150)(a) |

1159 | 5.1415926535897932384626433832795028841971694 |

1160 | |

1161 | The precision of ``a`` can be increased on demand since, internally, ``a`` |

1162 | keeps the exact representation it comes from: |

1163 | |

1164 | :: |

1165 | |

1166 | sage: a._value |

1167 | pi + 2 |

1168 | |

1169 | |

1170 | Its use is sometimes tricky, and shares some weirdnesses of the symbolic |

1171 | ring since it also warps its elements: |

1172 | |

1173 | :: |

1174 | |

1175 | sage: RLF(pi) + RLF(QQbar(2).sqrt()) |

1176 | 4.555806215962889? |

1177 | sage: _._value |

1178 | pi + 1.414213562373095? |

1179 | |

1180 |