# Ticket #15944: real_and_complex_numbers.ipynb

File real_and_complex_numbers.ipynb, 65.3 KB (added by tmonteil, 4 years ago)
Line
1{
2 "nbformat_minor": 2,
3 "nbformat": 4,
4 "cells": [
5  {
6   "source": [
7    "$$\n", 8 "\\def\\CC{\\bf C}\n", 9 "\\def\\QQ{\\bf Q}\n", 10 "\\def\\RR{\\bf R}\n", 11 "\\def\\ZZ{\\bf Z}\n", 12 "\\def\\NN{\\bf N}\n", 13 "$$\n",
14    "# Tutorial: Real and complex numbers in Sage\n",
15    "\n",
16    "Even with unlimitted memory, computers can not represent all real or complex numbers, since those form uncountable sets. Sage proposes various classes to approach the set of real numbers. They are considered as Fields or Rings though they do not necessarily satisfy their respective axiomatic.\n",
17    "\n",
18    "Those various classes have very different features, hence it is important to know to which one a given number belongs to understand its behaviour.\n",
19    "\n",
20    "You should know a bit about coercion before starting this tutorial, so it is advised to read the parent_element_coercion.en.rst worksheet first.\n",
21    "\n",
22    "## The big zoo\n",
23    "\n",
24    "Here is a rough qualitative summary of the main representations that will appear in this tutorial:\n",
25    "\n",
26    "| Type      | Name                      | Short | Complex equivalent           | Short   | Reliable | Fast    | Expressive | Exact | Precision                     |\n",
27    "|-----------|---------------------------|-------|------------------------------|---------|----------|---------|------------|-------|-------------------------------|\n",
28    "| Numerical | RealDoubleField()       | RDF | ComplexDoubleField()       | CDF   | +++    | +++++ | +        | No    | bounded by 53 bits            |\n",
29    "|           | RealField(prec)         | RR  | ComplexField(prec)         | CC    | +++    | +++   | ++       | No    | finite unbounded non-adaptive |\n",
30    "|           | RealIntervalField(prec) | RIF | ComplexIntervalField(prec) | CIF   | +++++  | +++   | ++       | No    | finite unbounded non-adaptive |\n",
31    "|           | RealBallField(prec)     | RBF | ComplexBallField(prec)     | CBF   | +++++  | +++   | ++       | No    | finite unbounded non-adaptive |\n",
32    "| Algebraic | AlgebraicRealField()    | AA  | AlgebraicField()           | QQbar | +++++  | +     | ++++     | Yes   | infinite                      |\n",
33    "|           | NumberField(poly)       | none  | itself                       | none    | +++++  | ++    | ++++     | Yes   | infinite                      |\n",
34    "|           | QuadraticField(n)       | none  | itself                       | none    | +++++  | +++   | +++      | Yes   | infinite                      |\n",
35    "|           | RationalField()         | QQ  | QuadraticField(-1)         | none    | +++++  | ++++  | +++      | Yes   | infinite                      |\n",
36    "| Symbolic  | SymbolicRing()          | SR  | itself                       | SR    | +      | +     | +++++    | No    | depends on construction       |\n",
37    "\n",
38    "Numerical representations are fast but inexact, they require some special care since the precision may become fuzzy along a computation.\n",
39    "\n",
40    "Algebraic representations are exact and reliable, but slower and able to represent only a class of numbers.\n",
41    "\n",
42    "Symbolic representations are very versatile (they can represent numbers such as $\\sqrt{log(\\pi)+sin(42)}$), but slow and sometimes not reliable.\n",
43    "\n",
44    "Let us inspect each kind of representation further.\n",
45    "\n",
46    "**TODO**\n",
47    "\n",
48    "Find a better word for \"kind\" that is not \"type\", \"group\", \"category\", \"class\" (perhaps \"taxon\", \"family\", \"sort\" ?).\n",
49    "\n",
50    "## Numerical representations\n",
51    "\n",
52    "### Floating-point numbers\n",
53    "\n",
54    "Floating-point numbers are real numbers of the form $(-1)^s 0.d_0 d_1 d_2 \\dots\n", 55 "d_{p-1} 2^{e}$ where:\n",
56    "\n",
57    "-   $s \\in \\{0,1\\}$ is the *sign*,\n",
58    "-   $d_i \\in \\{0,1\\}$ (with $d_0=1$) are the *digits*,\n",
59    "-   the integer $d_0 d_1 d_2 \\dots d_{p-1}$ is the *mantissa*,\n",
60    "-   $p$ is the *precision*,\n",
61    "-   $e$ is the *exponent*.\n",
62    "\n",
63    "Floating-point numbers of given precision do not form a field in the mathematical sense since they are not stable under the classical operations. So when adding or multiplying two floating-point numbers, Sage has to round the result to return a floating-point number. In particular, the results are approximate and the error can grow along a computation.\n",
64    "\n",
65    "Some first warnings about floating-point arithmetics can be found on the guide: <http://floating-point-gui.de/>\n",
66    "\n",
67    "Error propagation is sometimes unavoidable, but it is possible to know a bound on it by using the following secure fields:\n",
68    "\n",
69    "There exists various floating-point implementations of numbers in Sage:\n",
70    "\n",
71    "#### Real Double Field\n",
72    "\n",
73    "Elements of RDF = RealDoubleField() are the floating-point numbers in double precision from the processor (those you use when you code in C, similar to Python $float$), see sage.rings.real\\_double. The computations using those numbers are very fast and most optimized libraries that work with floating-point numbers use this representation, so you will benefit from them if the number involved in your floating-point constructions (matrices, polynomials,...) belong to RDF. The mantissa has 53 bits of precision, the exponent is encoded on 11 bits from -1023 to 1024. During the computations, the rounding is done toward the closest floating-point number.\n",
74    "\n",
75    "#### Real Fields\n",
76    "\n",
77    "Elements of RealField(prec) are the floating-point numbers with prec bits of precision, \"emulated\" by the MPFR library. By default, the mantissa has $53$ bits of precision (and RealField(prec=53) can be shortened as RR), the exponent is encoded on 32 bits (so that you will unlikely overflow), see sage.rings.real\\_mpfr. The computations are slower than in RDF and most algorithms on objects whose base ring is RealField(prec) will be pretty naive. However, it is possible to extend the precision as well as the type of rounding (towards $+\\infty$, towards $-\\infty$, towards $0$), for example RealField(100, rnd='RNDZ') is the \"field\" of real floating-point numbers with 100 bits of precision with rounding towards zero.\n",
78    "\n",
79    "> **note**\n",
80    ">\n",
81    "> While sometimes needed, increasing the precision also increases the computation time.\n",
82    "\n",
83    "> **warning**\n",
84    ">\n",
85    "> Note that despite its generic name, RealField(prec) or RR is *not* an abstraction of the field of real numbers, as can be seen by the following section.\n",
86    "\n",
87    "#### Unreal floats\n",
88    "\n",
89    "Besides not representing all real numbers, RDF and RR contain three additional unreal elements that allow the (+,-,*,/) operations to be always defined (NaN stands for *Not a Number*):"
90   ],
91   "cell_type": "markdown",
93  },
94  {
95   "execution_count": null,
96   "cell_type": "code",
97   "source": [
98    "RDF(0) / RDF(0)"
99   ],
100   "outputs": [
101    {
102     "execution_count": 1,
103     "output_type": "execute_result",
104     "data": {
105      "text/plain": [
106       "NaN\n"
107      ]
108     },
110    }
111   ],
113  },
114  {
115   "execution_count": null,
116   "cell_type": "code",
117   "source": [
118    "RDF(1) / RDF(0)"
119   ],
120   "outputs": [
121    {
122     "execution_count": 1,
123     "output_type": "execute_result",
124     "data": {
125      "text/plain": [
126       "+infinity\n"
127      ]
128     },
130    }
131   ],
133  },
134  {
135   "execution_count": null,
136   "cell_type": "code",
137   "source": [
138    "RR(-1) / RR(0)"
139   ],
140   "outputs": [
141    {
142     "execution_count": 1,
143     "output_type": "execute_result",
144     "data": {
145      "text/plain": [
146       "-infinity\n"
147      ]
148     },
150    }
151   ],
153  },
154  {
155   "execution_count": null,
156   "cell_type": "code",
157   "source": [
158    "RDF(infinity) - RDF(infinity)"
159   ],
160   "outputs": [
161    {
162     "execution_count": 1,
163     "output_type": "execute_result",
164     "data": {
165      "text/plain": [
166       "NaN\n"
167      ]
168     },
170    }
171   ],
173  },
174  {
175   "execution_count": null,
176   "cell_type": "code",
177   "source": [
178    "NaN in RR and Infinity in RR and -Infinity in RR"
179   ],
180   "outputs": [
181    {
182     "execution_count": 1,
183     "output_type": "execute_result",
184     "data": {
185      "text/plain": [
186       "True"
187      ]
188     },
190    }
191   ],
193  },
194  {
195   "source": [
196    "#### Back and forth in precision\n",
197    "\n",
198    "Note that it is is general useless to improve the precision during a computation. This explains why the coercion is done toward the smallest precision:"
199   ],
200   "cell_type": "markdown",
202  },
203  {
204   "execution_count": null,
205   "cell_type": "code",
206   "source": [
207    "R = RealField(200)\n",
208    "S = RealField(100)\n",
209    "a = R(2).sqrt()\n",
210    "b = S(pi)\n",
211    "a+b == R(sqrt(2)+pi)"
212   ],
213   "outputs": [
214    {
215     "execution_count": 1,
216     "output_type": "execute_result",
217     "data": {
218      "text/plain": [
219       "True"
220      ]
221     },
223    }
224   ],
226  },
227  {
228   "execution_count": null,
229   "cell_type": "code",
230   "source": [
231    "a+R(b) == R(sqrt(2)+pi)"
232   ],
233   "outputs": [
234    {
235     "execution_count": 1,
236     "output_type": "execute_result",
237     "data": {
238      "text/plain": [
239       "False"
240      ]
241     },
243    }
244   ],
246  },
247  {
248   "source": [
249    "> **note**\n",
250    ">\n",
251    "> For each element of those two fields, there exists a .sign_mantissa_exponent() method that returns the corresponding triple.\n",
252    "\n",
253    "### Floating-point numbers with a control on the imprecision\n",
254    "\n",
255    "#### Real Interval Fields\n",
256    "\n",
257    "Elements of RealIntervalField(prec) are pairs of elements of RealField(prec), which should be interpreted as the two endpoints of an interval containing the represented real number. RealIntervalField(53) can be shortened as RIF.\n",
258    "\n",
259    "Operations defined over the RealIntervalField implementation must offer the *containment guaranty*: for every mathematical function $f : \\mathbb{R}^d \\to\n", 260 "\\mathbb{R}$ that admits an implementation f_rif over some RealIntervalField(prec), we have the property that for every element x = [(a_0, b_0), (a_1,b_1), ..., (a_(d-1),b_(d-1))] in RealIntervalField(prec)^d, if f_rif(x) = [a', b'], then $f([a_0, b_0]\n", 261 "\\times [a_1,b_1] \\times \\dots \\times [a_(d-1),b_(d-1)])\\subseteq [a', b']$. In particular, the left endpoint of the interval is rounded towards $-\\infty$ and the right endpoint is rounded towards $+\\infty$\n",
262    "\n",
263    "#### Real Ball Fields\n",
264    "\n",
265    "Elements of RealBallField(prec) are pairs (c,r), where c is a floating-point number with prec bits of precision, and where r is a floating-point number with 30 bits of precision. The high precision c and the low precision number r should be understood as the center and the radius (respectively) of a ball containing the represented real number. RealBallField(53) can be shortened as RBF.\n",
266    "\n",
267    "Operations defined over the RealBallField implementation must offer a similar containment guaranty as for RealIntervalField (replace intervals with balls).\n",
268    "\n",
269    "#### Comparison between both representations\n",
270    "\n",
271    "The main advantage of the ball implementation over interval arithmetics is that, when the precision is high, the error is still encoded with a low precision (since only the order of magnitude matters for the error bounds), so that the loss of speed compared to non-guaranteed floating-point arithmetics with the same precision is about $1+\\varepsilon$ (and not $2$ as in the case of interval arithmetics, where two precise numbers have to be computed).\n",
272    "\n",
273    "Also, it is closer to the mathematical $\\varepsilon-\\delta$ definitions and distance-based error estimations.\n",
274    "\n",
275    "However, for large intervals and for evaluating monotonic functions, interval arithmetics is usually more precise and easier to implement (just round the images of the endpoints towards the right direction). It is also well-adapted for implementing dichotomy-based precision refinements.\n",
276    "\n",
277    "### Which numerical representation to choose\n",
278    "\n",
279    "Sage interfaces various libraries written in C, C++, Fortran, etc. which use the processor double precision floats. Hence, the use of RDF has the advantage of being faster and allowing the use of those libraries.\n",
280    "\n",
281    "If you use other floating-point representations, Sage will mostly use generic algorithms, which are usually slower and numerically less stable. However, those representations allow you to express real numbers with a higher precision.\n",
282    "\n",
283    "Regarding speed, RealBallField(prec) and RealField(prec) are similar, while RealIntervalField(prec) is twice slower. However, the RealBallField(prec) is a new feature of Sage, hence some methods might raise a NotImplementedError.\n",
284    "\n",
285    "**Summary**\n",
286    "\n",
287    "So, if you need speed and the power of Sage libraries: use RDF.\n",
288    "\n",
289    "If you need more precision, start with RealBallField(prec), and try RealIntervalField(prec) in case the imprecision becomes too large (alternatively: increase the precision or try to think about the order and the expanding nature of your operations). Use RealField(prec) if the methods you wanted to use with RealBallField(prec) are not implemented.\n",
290    "\n",
291    "### Exercises\n",
292    "\n",
293    "#### Binomial\n",
294    "\n",
295    "From the three following values:"
296   ],
297   "cell_type": "markdown",
299  },
300  {
301   "execution_count": null,
302   "cell_type": "code",
303   "source": [
304    "a = 1.0\n",
305    "b = 10.0^4\n",
306    "c = 1.0"
307   ],
308   "outputs": [],
310  },
311  {
312   "source": [
313    "solve the $ax^2+bx+c = 0$ equation by the method you learned in highschool involving discriminant, and compare the sum and product of the solution you found to their theoretical values $-b/a$ and $c/a$."
314   ],
315   "cell_type": "markdown",
317  },
318  {
319   "execution_count": null,
320   "cell_type": "code",
321   "source": [
322    "# edit here"
323   ],
324   "outputs": [],
326  },
327  {
328   "source": [
329    "Same question by letting Sage solve this equation for you."
330   ],
331   "cell_type": "markdown",
333  },
334  {
335   "execution_count": null,
336   "cell_type": "code",
337   "source": [
338    "# edit here"
339   ],
340   "outputs": [],
342  },
343  {
344   "source": [
345    "Even on such a simple computation, you have been victim of a phenomenon of *error amplification*. Another striking example is the following:"
346   ],
347   "cell_type": "markdown",
349  },
350  {
351   "execution_count": null,
352   "cell_type": "code",
353   "source": [
354    "a = 10000.0\n",
355    "b = -9999.5\n",
356    "c = 0.1\n",
357    "a+c+b"
358   ],
359   "outputs": [
360    {
361     "execution_count": 1,
362     "output_type": "execute_result",
363     "data": {
364      "text/plain": [
365       "0.600000000000364"
366      ]
367     },
369    }
370   ],
372  },
373  {
374   "execution_count": null,
375   "cell_type": "code",
376   "source": [
377    "a+b+c"
378   ],
379   "outputs": [
380    {
381     "execution_count": 1,
382     "output_type": "execute_result",
383     "data": {
384      "text/plain": [
385       "0.600000000000000"
386      ]
387     },
389    }
390   ],
392  },
393  {
394   "source": [
395    "Note that, by construction and following the IEEE 754 specification, the + (and *) law is commutative in RDF and RR :"
396   ],
397   "cell_type": "markdown",
399  },
400  {
401   "execution_count": null,
402   "cell_type": "code",
403   "source": [
404    "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)])"
405   ],
406   "outputs": [
407    {
408     "execution_count": 1,
409     "output_type": "execute_result",
410     "data": {
411      "text/plain": [
412       "True"
413      ]
414     },
416    }
417   ],
419  },
420  {
421   "source": [
422    "Actually the associativity is not satisfied here:"
423   ],
424   "cell_type": "markdown",
426  },
427  {
428   "execution_count": null,
429   "cell_type": "code",
430   "source": [
431    "(a+b)+c == a+(b+c)"
432   ],
433   "outputs": [
434    {
435     "execution_count": 1,
436     "output_type": "execute_result",
437     "data": {
438      "text/plain": [
439       "False"
440      ]
441     },
443    }
444   ],
446  },
447  {
448   "source": [
449    "Try to find, by random sampling, 3 elements $a,b,c$ in RDF such that (a+b)+c != a+(b+c) :"
450   ],
451   "cell_type": "markdown",
453  },
454  {
455   "execution_count": null,
456   "cell_type": "code",
457   "source": [
458    "# edit here"
459   ],
460   "outputs": [],
462  },
463  {
464   "source": [
465    "Why is it so hard ? What makes the previous example work ?\n",
466    "\n",
467    "#### Interval arithmetics\n",
468    "\n",
469    "\"Prove\" using RealIntervalField that $cos(2\\pi)=1$ with a high precision (say $10^{-100}$)."
470   ],
471   "cell_type": "markdown",
473  },
474  {
475   "execution_count": null,
476   "cell_type": "code",
477   "source": [
478    "# edit here "
479   ],
480   "outputs": [],
482  },
483  {
484   "source": [
485    "Explain the following behaviour:"
486   ],
487   "cell_type": "markdown",
489  },
490  {
491   "execution_count": null,
492   "cell_type": "code",
493   "source": [
494    "sqrt(0.1)^2 == 0.1"
495   ],
496   "outputs": [
497    {
498     "execution_count": 1,
499     "output_type": "execute_result",
500     "data": {
501      "text/plain": [
502       "True"
503      ]
504     },
506    }
507   ],
509  },
510  {
511   "execution_count": null,
512   "cell_type": "code",
513   "source": [
514    "sqrt(RIF(0.1))^2 == RIF(0.1)"
515   ],
516   "outputs": [
517    {
518     "execution_count": 1,
519     "output_type": "execute_result",
520     "data": {
521      "text/plain": [
522       "False"
523      ]
524     },
526    }
527   ],
529  },
530  {
531   "source": [
532    "Provide a relation between sqrt(RIF(0.1))^2 and RIF(0.1)"
533   ],
534   "cell_type": "markdown",
536  },
537  {
538   "execution_count": null,
539   "cell_type": "code",
540   "source": [
541    "# edit here"
542   ],
543   "outputs": [],
545  },
546  {
547   "source": [
548    "Note that the ordering of the computations has an influence on the result, but also on its precision. With the previous example, we have:"
549   ],
550   "cell_type": "markdown",
552  },
553  {
554   "execution_count": null,
555   "cell_type": "code",
556   "source": [
557    "a = RIF(10000.0)\n",
558    "b = RIF(9999.5)\n",
559    "c = RIF(0.1)\n",
560    "a+c-b"
561   ],
562   "outputs": [
563    {
564     "execution_count": 1,
565     "output_type": "execute_result",
566     "data": {
567      "text/plain": [
568       "0.60000000000?"
569      ]
570     },
572    }
573   ],
575  },
576  {
577   "execution_count": null,
578   "cell_type": "code",
579   "source": [
580    "a-b+c"
581   ],
582   "outputs": [
583    {
584     "execution_count": 1,
585     "output_type": "execute_result",
586     "data": {
587      "text/plain": [
588       "0.6000000000000000?\n"
589      ]
590     },
592    }
593   ],
595  },
596  {
597   "execution_count": null,
598   "cell_type": "code",
599   "source": [
600    "(a+c-b).endpoints()"
601   ],
602   "outputs": [
603    {
604     "execution_count": 1,
605     "output_type": "execute_result",
606     "data": {
607      "text/plain": [
608       "(0.599999999998544, 0.600000000000364)"
609      ]
610     },
612    }
613   ],
615  },
616  {
617   "execution_count": null,
618   "cell_type": "code",
619   "source": [
620    "(a+c-b).absolute_diameter()"
621   ],
622   "outputs": [
623    {
624     "execution_count": 1,
625     "output_type": "execute_result",
626     "data": {
627      "text/plain": [
628       "1.81898940354586e-12\n"
629      ]
630     },
632    }
633   ],
635  },
636  {
637   "execution_count": null,
638   "cell_type": "code",
639   "source": [
640    "(a-b+c).endpoints()"
641   ],
642   "outputs": [
643    {
644     "execution_count": 1,
645     "output_type": "execute_result",
646     "data": {
647      "text/plain": [
648       "(0.599999999999999, 0.600000000000001)"
649      ]
650     },
652    }
653   ],
655  },
656  {
657   "execution_count": null,
658   "cell_type": "code",
659   "source": [
660    "(a-b+c).absolute_diameter()"
661   ],
662   "outputs": [
663    {
664     "execution_count": 1,
665     "output_type": "execute_result",
666     "data": {
667      "text/plain": [
668       "1.11022302462516e-16"
669      ]
670     },
672    }
673   ],
675  },
676  {
677   "source": [
678    "Before evaluating them, could you predict and explain the behaviour of the following commands ?"
679   ],
680   "cell_type": "markdown",
682  },
683  {
684   "execution_count": null,
685   "cell_type": "code",
686   "source": [
687    "# 1/2 in RIF"
688   ],
689   "outputs": [],
691  },
692  {
693   "execution_count": null,
694   "cell_type": "code",
695   "source": [
696    "# 1/3 in RIF"
697   ],
698   "outputs": [],
700  },
701  {
702   "source": [
703    "*Hint*: you can look at the documentation of the .__contains__() method of RIF.\n",
704    "\n",
705    "A Sage developer wants to write the real function $f$ that maps $0.01$ to $1$ and that is the identity elsewhere, and suggests the following implementation:"
706   ],
707   "cell_type": "markdown",
709  },
710  {
711   "execution_count": null,
712   "cell_type": "code",
713   "source": [
714    "def f_unstable(x):\n",
715    "    if x == 0.01:\n",
716    "        return 1\n",
717    "    else: \n",
718    "        return x"
719   ],
720   "outputs": [],
722  },
723  {
724   "source": [
725    "Prove that this implementation does not offers the containment guaranty when used on the RealIntervalField.\n",
726    "\n",
727    "Could you help this lost developer to find a correct implementation so that it could be included in Sage source code ?"
728   ],
729   "cell_type": "markdown",
731  },
732  {
733   "execution_count": null,
734   "cell_type": "code",
735   "source": [
736    "def f_correct(x):\n",
737    "    # edit here"
738   ],
739   "outputs": [],
741  },
742  {
743   "source": [
744    "#### Divergent series ?\n",
745    "\n",
746    "The series $u_n := 1 + 1/2 + 1/3 + 1/4 + \\dots + 1/n = \\sum_{k=1}^{n} 1/k$ is known to diverge on $\\mathbb{R}$.\n",
747    "\n",
748    "Explore its convergence on some $Real Fields$ (start with small precision). Explore the difference with the (mathematically equal) series $v_n := 1/n +\n", 749 "1/(n-1) + ... + 1/2 + 1$.\n",
750    "\n",
751    "Explain, formulate conjectures and try to prove them.\n",
752    "\n",
753    "#### A picture\n",
754    "\n",
755    "Draw (by hand) a picture of the elements of RealField(3).\n",
756    "\n",
757    "#### Convergence of a Markov Chain\n",
758    "\n",
759    "Let \\$M\\$ be the rational bistochastic matrix defined by:"
760   ],
761   "cell_type": "markdown",
763  },
764  {
765   "execution_count": null,
766   "cell_type": "code",
767   "source": [
768    "M = matrix(QQ, [[1/6, 1/2, 1/3], [1/2, 1/4, 1/4], [1/3, 1/4, 5/12]])\n",
769    "M"
770   ],
771   "outputs": [
772    {
773     "execution_count": 1,
774     "output_type": "execute_result",
775     "data": {
776      "text/plain": [
777       "[ 1/6  1/2  1/3]\n",
778       "[ 1/2  1/4  1/4]\n",
779       "[ 1/3  1/4 5/12]"
780      ]
781     },
783    }
784   ],
786  },
787  {
788   "source": [
789    "Could you guess the limit of \\$M^n\\$ when \\$n\\$ goes to \\$+infty\\$ ?"
790   ],
791   "cell_type": "markdown",
793  },
794  {
795   "execution_count": null,
796   "cell_type": "code",
797   "source": [
798    "# edit here"
799   ],
800   "outputs": [],
802  },
803  {
804   "source": [
805    "Could you prove that this limit exists and give its exact value ?"
806   ],
807   "cell_type": "markdown",
809  },
810  {
811   "execution_count": null,
812   "cell_type": "code",
813   "source": [
814    "# edit here"
815   ],
816   "outputs": [],
818  },
819  {
820   "source": [
821    "#### Repulsive and attractive fixed points\n",
822    "\n",
823    "Define the maps $f:x\\to 4x-1$ and $g:x\\to x/3+1$."
824   ],
825   "cell_type": "markdown",
827  },
828  {
829   "execution_count": null,
830   "cell_type": "code",
831   "source": [
832    "# edit here"
833   ],
834   "outputs": [],
836  },
837  {
838   "source": [
839    "Let $u_n$ and $v_n$ be the sequences defined recursively by $u_0 = RDF(1/3)$, $u_{n+1} = f(u_n)$ and $v_0 = RDF(3/2)$, $v_{n+1} = g(v_n)$.\n",
840    "\n",
841    "Compute the first 100 values of the sequences $u_n$ and $v_n$."
842   ],
843   "cell_type": "markdown",
845  },
846  {
847   "execution_count": null,
848   "cell_type": "code",
849   "source": [
850    "# edit here"
851   ],
852   "outputs": [],
854  },
855  {
856   "source": [
857    "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."
858   ],
859   "cell_type": "markdown",
861  },
862  {
863   "execution_count": null,
864   "cell_type": "code",
865   "source": [
866    "# edit here"
867   ],
868   "outputs": [],
870  },
871  {
872   "source": [
873    "#### Decimals of pi\n",
874    "\n",
875    "What is the 150th decimal of $\\pi$ ?"
876   ],
877   "cell_type": "markdown",
879  },
880  {
881   "execution_count": null,
882   "cell_type": "code",
883   "source": [
884    "# edit here"
885   ],
886   "outputs": [],
888  },
889  {
890   "source": [
891    "#### Money, money, money"
892   ],
893   "cell_type": "markdown",
895  },
896  {
897   "execution_count": null,
898   "cell_type": "code",
899   "source": [
900    "one_cent = 0.01\n",
901    "my_money = 0\n",
902    "for i in range(100):\n",
903    "    my_money += one_cent\n",
904    "my_money > 1"
905   ],
906   "outputs": [],
908  },
909  {
910   "source": [
911    "What is the origin of this behaviour ?\n",
912    "\n",
913    "How to fix the issue in Sage (respecively un pure Python, compare both) ?"
914   ],
915   "cell_type": "markdown",
917  },
918  {
919   "execution_count": null,
920   "cell_type": "code",
921   "source": [
922    "# edit here"
923   ],
924   "outputs": [],
926  },
927  {
928   "source": [
929    "Which risk does computer science take if finance and high-frequency trading stays on power for a while ?\n",
930    "\n",
931    "Provide the first 30 bits of the binary expansion of $0.01$."
932   ],
933   "cell_type": "markdown",
935  },
936  {
937   "execution_count": null,
938   "cell_type": "code",
939   "source": [
940    "# edit here"
941   ],
942   "outputs": [],
944  },
945  {
946   "source": [
947    "#### Huge matrix computation\n",
948    "\n",
949    "Construct a dense matrix of size 500*500 with real coefficients randomly distributed in $[0,1]$."
950   ],
951   "cell_type": "markdown",
953  },
954  {
955   "execution_count": null,
956   "cell_type": "code",
957   "source": [
958    "# edit here"
959   ],
960   "outputs": [],
962  },
963  {
964   "source": [
965    "Take the 100th power of this matrix in less than one second."
966   ],
967   "cell_type": "markdown",
969  },
970  {
971   "execution_count": null,
972   "cell_type": "code",
973   "source": [
974    "# edit here"
975   ],
976   "outputs": [],
978  },
979  {
980   "source": [
981    "Comment.\n",
982    "\n",
983    "#### Exponents"
984   ],
985   "cell_type": "markdown",
987  },
988  {
989   "execution_count": null,
990   "cell_type": "code",
991   "source": [
992    "RDF(1/2^10000) == 0"
993   ],
994   "outputs": [
995    {
996     "execution_count": 1,
997     "output_type": "execute_result",
998     "data": {
999      "text/plain": [
1000       "True"
1001      ]
1002     },
1004    }
1005   ],
1007  },
1008  {
1009   "execution_count": null,
1010   "cell_type": "code",
1011   "source": [
1012    "RR(1/2^10000) == 0"
1013   ],
1014   "outputs": [
1015    {
1016     "execution_count": 1,
1017     "output_type": "execute_result",
1018     "data": {
1019      "text/plain": [
1020       "False"
1021      ]
1022     },
1024    }
1025   ],
1027  },
1028  {
1029   "source": [
1030    "Comment.\n",
1031    "\n",
1032    "#### A strange impression"
1033   ],
1034   "cell_type": "markdown",
1036  },
1037  {
1038   "execution_count": null,
1039   "cell_type": "code",
1040   "source": [
1041    "RR(1/pi)"
1042   ],
1043   "outputs": [
1044    {
1045     "execution_count": 1,
1046     "output_type": "execute_result",
1047     "data": {
1048      "text/plain": [
1049       "0.318309886183791"
1050      ]
1051     },
1053    }
1054   ],
1056  },
1057  {
1058   "execution_count": null,
1059   "cell_type": "code",
1060   "source": [
1061    "RDF(1/pi)"
1062   ],
1063   "outputs": [
1064    {
1065     "execution_count": 1,
1066     "output_type": "execute_result",
1067     "data": {
1068      "text/plain": [
1069       "0.3183098861837907"
1070      ]
1071     },
1073    }
1074   ],
1076  },
1077  {
1078   "source": [
1079    "Before evaluating the next command, could you tell what will be the answer of the following command ?"
1080   ],
1081   "cell_type": "markdown",
1083  },
1084  {
1085   "execution_count": null,
1086   "cell_type": "code",
1087   "source": [
1088    "# RR(1/pi) == RDF(1/pi)"
1089   ],
1090   "outputs": [],
1092  },
1093  {
1094   "source": [
1095    "What is the cause of the behaviour of the first two commands ?\n",
1096    "\n",
1097    "#### Rounding\n",
1098    "\n",
1099    "Let us define:"
1100   ],
1101   "cell_type": "markdown",
1103  },
1104  {
1105   "execution_count": null,
1106   "cell_type": "code",
1107   "source": [
1108    "n = 5.0000000000000001"
1109   ],
1110   "outputs": [],
1112  },
1113  {
1114   "source": [
1115    "Explain each of the following behaviours:"
1116   ],
1117   "cell_type": "markdown",
1119  },
1120  {
1121   "execution_count": null,
1122   "cell_type": "code",
1123   "source": [
1124    "n"
1125   ],
1126   "outputs": [
1127    {
1128     "execution_count": 1,
1129     "output_type": "execute_result",
1130     "data": {
1131      "text/plain": [
1132       "5.000000000000000"
1133      ]
1134     },
1136    }
1137   ],
1139  },
1140  {
1141   "execution_count": null,
1142   "cell_type": "code",
1143   "source": [
1144    "n.str(truncate=False)"
1145   ],
1146   "outputs": [
1147    {
1148     "execution_count": 1,
1149     "output_type": "execute_result",
1150     "data": {
1151      "text/plain": [
1152       "'5.000000000000000111'"
1153      ]
1154     },
1156    }
1157   ],
1159  },
1160  {
1161   "execution_count": null,
1162   "cell_type": "code",
1163   "source": [
1164    "5.000000000000000111 == 5.0000000000000001"
1165   ],
1166   "outputs": [],
1168  },
1169  {
1170   "source": [
1171    "## Algebraic (exact) representations\n",
1172    "\n",
1173    "Some subfields of the real numbers can be represented in an exact way in Sage. They are pretty safe (no rounding problem, well defined semantics) through sometimes slow.\n",
1174    "\n",
1175    "As a general rule, the smallest the field is (for the inclusion), the fastest it is, which means that, in principle, working within rational field is faster than within quadratic fields, which is faster than within cyclotomic fields, which is faster than within number fields, which is faster that within the universal cyclotomic field, that is faster than within the whole algebraic field.\n",
1176    "\n",
1177    "Here is a short overview:\n",
1178    "\n",
1179    "### Rational numbers\n",
1180    "\n",
1181    "QQ = RationalField() stands for the field of rational numbers. Computations are exact and only limited by the available memory. You should note that the numerators/denominators could become very large along a computation, which might slow down this latter."
1182   ],
1183   "cell_type": "markdown",
1185  },
1186  {
1187   "execution_count": null,
1188   "cell_type": "code",
1189   "source": [
1190    "(123/321)^10000"
1191   ],
1192   "outputs": [],
1194  },
1195  {
1196   "source": [
1197    "It is possible to \"approximate\" a floating-point real number by rational numbers with the methods nearby_rational(), simplest_rational(), exact_rational()."
1198   ],
1199   "cell_type": "markdown",
1201  },
1202  {
1203   "execution_count": null,
1204   "cell_type": "code",
1205   "source": [
1206    "# edit here"
1207   ],
1208   "outputs": [],
1210  },
1211  {
1212   "source": [
1213    "If x belongs to RR, what is the result of QQ(x) ? Is it a coercion ?\n",
1214    "\n",
1215    "See the page <http://doc.sagemath.org/html/en/reference/rings_standard/sage/rings/rational_field.html>\n",
1216    "\n",
1217    "### Algebraic number field\n",
1218    "\n",
1219    "AA = AlgebraicRealField() stands for te field of real algebraic numbers."
1220   ],
1221   "cell_type": "markdown",
1223  },
1224  {
1225   "execution_count": null,
1226   "cell_type": "code",
1227   "source": [
1228    "a = AA(2).sqrt() + AA(7).sqrt()\n",
1229    "a"
1230   ],
1231   "outputs": [
1232    {
1233     "execution_count": 1,
1234     "output_type": "execute_result",
1235     "data": {
1236      "text/plain": [
1237       "4.059964873437685?"
1238      ]
1239     },
1241    }
1242   ],
1244  },
1245  {
1246   "source": [
1247    "Here, the number a is *printed* as its approximated numerical value, but the question mark at the end indicates that it is not internally a floating-point number.\n",
1248    "\n",
1249    "Basically, an algebraic number is represented as the root of an integer polynomial, and a small interval surrounding the appropriate root:"
1250   ],
1251   "cell_type": "markdown",
1253  },
1254  {
1255   "execution_count": null,
1256   "cell_type": "code",
1257   "source": [
1258    "a.minpoly()"
1259   ],
1260   "outputs": [
1261    {
1262     "execution_count": 1,
1263     "output_type": "execute_result",
1264     "data": {
1265      "text/plain": [
1266       "x^4 - 18*x^2 + 25\n"
1267      ]
1268     },
1270    }
1271   ],
1273  },
1274  {
1275   "execution_count": null,
1276   "cell_type": "code",
1277   "source": [
1278    "a.n(digits=50)"
1279   ],
1280   "outputs": [
1281    {
1282     "execution_count": 1,
1283     "output_type": "execute_result",
1284     "data": {
1285      "text/plain": [
1286       "4.0599648734376856393033044778489585042799310584594"
1287      ]
1288     },
1290    }
1291   ],
1293  },
1294  {
1295   "execution_count": null,
1296   "cell_type": "code",
1297   "source": [
1298    "a.interval_diameter(1e-40)"
1299   ],
1300   "outputs": [
1301    {
1302     "execution_count": 1,
1303     "output_type": "execute_result",
1304     "data": {
1305      "text/plain": [
1306       "4.0599648734376856393033044778489585042799310584593982535450141971918013016924?"
1307      ]
1308     },
1310    }
1311   ],
1313  },
1314  {
1315   "source": [
1316    "Computing such a representation is costly, hence Sage is lazy and computes it only when needed:"
1317   ],
1318   "cell_type": "markdown",
1320  },
1321  {
1322   "execution_count": null,
1323   "cell_type": "code",
1324   "source": [
1325    "b = QQbar(2).sqrt()\n",
1326    "c = a^2\n",
1327    "c"
1328   ],
1329   "outputs": [
1330    {
1331     "execution_count": 1,
1332     "output_type": "execute_result",
1333     "data": {
1334      "text/plain": [
1335       "2.000000000000000?"
1336      ]
1337     },
1339    }
1340   ],
1342  },
1343  {
1344   "source": [
1345    "Here, Sage did not compute the minimal polynomial nor an enclosing interval of c, it just knows c as \"the square of b\", in particular, it does not know yet that c is the integer 2. However, such an \"exactification\" is computed when required (and kept internally):"
1346   ],
1347   "cell_type": "markdown",
1349  },
1350  {
1351   "execution_count": null,
1352   "cell_type": "code",
1353   "source": [
1354    "c.is_integer()"
1355   ],
1356   "outputs": [
1357    {
1358     "execution_count": 1,
1359     "output_type": "execute_result",
1360     "data": {
1361      "text/plain": [
1362       "True"
1363      ]
1364     },
1366    }
1367   ],
1369  },
1370  {
1371   "execution_count": null,
1372   "cell_type": "code",
1373   "source": [
1374    "b"
1375   ],
1376   "outputs": [
1377    {
1378     "execution_count": 1,
1379     "output_type": "execute_result",
1380     "data": {
1381      "text/plain": [
1382       "2"
1383      ]
1384     },
1386    }
1387   ],
1389  },
1390  {
1391   "source": [
1392    "You can force the exactification with the exactify method.\n",
1393    "\n",
1394    "For more details, see the page <http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/qqbar.html>\n",
1395    "\n",
1396    "### Number fields\n",
1397    "\n",
1398    "See the pages:\n",
1399    "\n",
1400    "-   <http://doc.sagemath.org/html/en/reference/number_fields/index.html>\n",
1401    "-   <http://doc.sagemath.org/html/en/constructions/number_fields.html>\n",
1402    "-   <http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/number_field/number_field.html>\n",
1403    "\n",
1404    "> **warning**\n",
1405    ">\n",
1406    "> You should be aware that number field are purely algebraic objects that are *not* embedded into the complex plane by default. Do to this, you have to define the embedding you want to work with (if you need one).\n",
1407    "\n",
1408    "Note that the quadratic number fields have a dedicated faster implementation.\n",
1409    "\n",
1410    "What is the difference between those two objects ?"
1411   ],
1412   "cell_type": "markdown",
1414  },
1415  {
1416   "execution_count": null,
1417   "cell_type": "code",
1418   "source": [
1420    "qq = QQ[I]"
1421   ],
1422   "outputs": [],
1424  },
1425  {
1426   "source": [
1427    "### Cyclotomic fields (and the universal cyclotomic field)\n",
1428    "\n",
1429    "See the pages:\n",
1430    "\n",
1431    "-   <http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/number_field/number_field.html#sage.rings.number_field.number_field.CyclotomicFieldFactory>\n",
1432    "-   <http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/universal_cyclotomic_field.html>\n",
1433    "\n",
1434    "## Symbolic representations\n",
1435    "\n",
1436    "SR = sage.symbolic.ring.SymbolicRing() is the set of symbolic expressions.\n",
1437    "\n",
1438    "SR is pretty fuzzy but quite handy (here expressivity is prefered over standardization and the number of real numbers representable within SR is constantly increasing). This \"ring\" contains most mathematical constants such as $e$, $\\pi$, $cos(12)$, $\\sqrt{2}$, $\\gamma$ (Euler constant), ... and allows some computations such as:"
1439   ],
1440   "cell_type": "markdown",
1442  },
1443  {
1444   "execution_count": null,
1445   "cell_type": "code",
1446   "source": [
1447    "e^(I*pi)"
1448   ],
1449   "outputs": [
1450    {
1451     "execution_count": 1,
1452     "output_type": "execute_result",
1453     "data": {
1454      "text/plain": [
1455       "-1"
1456      ]
1457     },
1459    }
1460   ],
1462  },
1463  {
1464   "source": [
1465    "> **warning**\n",
1466    ">\n",
1467    "> SR also contains symbolic expressions that are not real or complex numbers:\n",
1468    ">\n",
1469    "> {.python .input}\n",
1470    "> e = cos(x)^2\n",
1471    "> e.derivative()\n",
1472    "> \n",
1473    ">\n",
1474    "> {.json .output}\n",
1475    "> [{\"execution_count\": 1, \"output_type\": \"execute_result\", \"data\": {\"text/plain\": \"-2*cos(x)*sin(x)\"}, \"metadata\": {}}]\n",
1476    "> \n",
1477    ">\n",
1478    "> {.python .input}\n",
1479    "> _ in SR\n",
1480    "> \n",
1481    ">\n",
1482    "> {.json .output}\n",
1483    "> [{\"execution_count\": 1, \"output_type\": \"execute_result\", \"data\": {\"text/plain\": \"True\"}, \"metadata\": {}}]\n",
1484    "> \n",
1485    "\n",
1486    "SR represents numbers as symbolic expressions, whose length (as a formula) can grow along a computation. This phenomenon can slow down the computations. If our aim is to end with a numerical value, you should use a floating-point representation, as they keep a constant length."
1487   ],
1488   "cell_type": "markdown",
1490  },
1491  {
1492   "execution_count": null,
1493   "cell_type": "code",
1494   "source": [
1495    "a = exp(2)\n",
1496    "f = lambda x : sqrt(x+1)^x\n",
1497    "for i in range(10):   \n",
1498    "    a = f(a)\n",
1499    "a"
1500   ],
1501   "outputs": [],
1503  },
1504  {
1505   "source": [
1506    "Of course, using floating-point arithmetics will not necessarily solve all your problems since you will have to ensure the numerical stability of such iterated computation (see the sections above).\n",
1507    "\n",
1508    "SR always answer, even when it does not know the answer, for example:"
1509   ],
1510   "cell_type": "markdown",
1512  },
1513  {
1514   "execution_count": null,
1515   "cell_type": "code",
1516   "source": [
1517    "# missing example"
1518   ],
1519   "outputs": [],
1521  },
1522  {
1523   "source": [
1524    "Hence, computations involving the symbolic ring should not be considered as reliable. In particular, SR is a dead-end for the coercion among numbers (and can add pears and apples formally):"
1525   ],
1526   "cell_type": "markdown",
1528  },
1529  {
1530   "execution_count": null,
1531   "cell_type": "code",
1532   "source": [
1533    "0.1 + pi + QQbar(2).sqrt()"
1534   ],
1535   "outputs": [
1536    {
1537     "execution_count": 1,
1538     "output_type": "execute_result",
1539     "data": {
1540      "text/plain": [
1541       "pi + 1.51421356237310"
1542      ]
1543     },
1545    }
1546   ],
1548  },
1549  {
1550   "execution_count": null,
1551   "cell_type": "code",
1552   "source": [
1553    "_.parent()"
1554   ],
1555   "outputs": [
1556    {
1557     "execution_count": 1,
1558     "output_type": "execute_result",
1559     "data": {
1560      "text/plain": [
1561       "Symbolic Ring"
1562      ]
1563     },
1565    }
1566   ],
1568  },
1569  {
1570   "source": [
1571    "Here are some issues that can appear with the symbolic ring, when you can represent your numbers in a different way, we strongly encourage you to do so:\n",
1572    "\n",
1574    "-   <https://trac.sagemath.org/wiki/symbolics>\n",
1579    "\n",
1580    "## Conversions, bridges, exactification\n",
1581    "\n",
1582    "**TODO**\n",
1583    "\n",
1584    "Should this be a separate section, or should it appear in the rational/algebraic sections ?\n",
1585    "\n",
1586    "Along a computation, coercion goes towards the less precise numbers:"
1587   ],
1588   "cell_type": "markdown",
1590  },
1591  {
1592   "execution_count": null,
1593   "cell_type": "code",
1594   "source": [
1595    "a = 1/10 ; a"
1596   ],
1597   "outputs": [
1598    {
1599     "execution_count": 1,
1600     "output_type": "execute_result",
1601     "data": {
1602      "text/plain": [
1603       "1/10"
1604      ]
1605     },
1607    }
1608   ],
1610  },
1611  {
1612   "execution_count": null,
1613   "cell_type": "code",
1614   "source": [
1615    "a.parent()"
1616   ],
1617   "outputs": [
1618    {
1619     "execution_count": 1,
1620     "output_type": "execute_result",
1621     "data": {
1622      "text/plain": [
1623       "Rational Field"
1624      ]
1625     },
1627    }
1628   ],
1630  },
1631  {
1632   "execution_count": null,
1633   "cell_type": "code",
1634   "source": [
1635    "b = QQbar(2).sqrt() ; b"
1636   ],
1637   "outputs": [
1638    {
1639     "execution_count": 1,
1640     "output_type": "execute_result",
1641     "data": {
1642      "text/plain": [
1643       "1.414213562373095?"
1644      ]
1645     },
1647    }
1648   ],
1650  },
1651  {
1652   "execution_count": null,
1653   "cell_type": "code",
1654   "source": [
1655    "b.parent()"
1656   ],
1657   "outputs": [
1658    {
1659     "execution_count": 1,
1660     "output_type": "execute_result",
1661     "data": {
1662      "text/plain": [
1663       "Algebraic Field"
1664      ]
1665     },
1667    }
1668   ],
1670  },
1671  {
1672   "execution_count": null,
1673   "cell_type": "code",
1674   "source": [
1675    "c = a+b ; c"
1676   ],
1677   "outputs": [
1678    {
1679     "execution_count": 1,
1680     "output_type": "execute_result",
1681     "data": {
1682      "text/plain": [
1683       "1.514213562373096?"
1684      ]
1685     },
1687    }
1688   ],
1690  },
1691  {
1692   "execution_count": null,
1693   "cell_type": "code",
1694   "source": [
1695    "c.parent()"
1696   ],
1697   "outputs": [
1698    {
1699     "execution_count": 1,
1700     "output_type": "execute_result",
1701     "data": {
1702      "text/plain": [
1703       "Algebraic Field"
1704      ]
1705     },
1707    }
1708   ],
1710  },
1711  {
1712   "execution_count": null,
1713   "cell_type": "code",
1714   "source": [
1715    "d = RDF(0.1) ; d"
1716   ],
1717   "outputs": [
1718    {
1719     "execution_count": 1,
1720     "output_type": "execute_result",
1721     "data": {
1722      "text/plain": [
1723       "0.1"
1724      ]
1725     },
1727    }
1728   ],
1730  },
1731  {
1732   "execution_count": null,
1733   "cell_type": "code",
1734   "source": [
1735    "d.parent()"
1736   ],
1737   "outputs": [
1738    {
1739     "execution_count": 1,
1740     "output_type": "execute_result",
1741     "data": {
1742      "text/plain": [
1743       "Real Double Field"
1744      ]
1745     },
1747    }
1748   ],
1750  },
1751  {
1752   "execution_count": null,
1753   "cell_type": "code",
1754   "source": [
1755    "e = c+d ; e"
1756   ],
1757   "outputs": [
1758    {
1759     "execution_count": 1,
1760     "output_type": "execute_result",
1761     "data": {
1762      "text/plain": [
1763       "1.614213562373095"
1764      ]
1765     },
1767    }
1768   ],
1770  },
1771  {
1772   "execution_count": null,
1773   "cell_type": "code",
1774   "source": [
1775    "e.parent()"
1776   ],
1777   "outputs": [
1778    {
1779     "execution_count": 1,
1780     "output_type": "execute_result",
1781     "data": {
1782      "text/plain": [
1783       "Complex Double Field"
1784      ]
1785     },
1787    }
1788   ],
1790  },
1791  {
1792   "execution_count": null,
1793   "cell_type": "code",
1794   "source": [
1795    "f = pi ; f"
1796   ],
1797   "outputs": [
1798    {
1799     "execution_count": 1,
1800     "output_type": "execute_result",
1801     "data": {
1802      "text/plain": [
1803       "pi"
1804      ]
1805     },
1807    }
1808   ],
1810  },
1811  {
1812   "execution_count": null,
1813   "cell_type": "code",
1814   "source": [
1815    "f.parent()"
1816   ],
1817   "outputs": [
1818    {
1819     "execution_count": 1,
1820     "output_type": "execute_result",
1821     "data": {
1822      "text/plain": [
1823       "Symbolic Ring"
1824      ]
1825     },
1827    }
1828   ],
1830  },
1831  {
1832   "execution_count": null,
1833   "cell_type": "code",
1834   "source": [
1835    "g = e+f ; g"
1836   ],
1837   "outputs": [
1838    {
1839     "execution_count": 1,
1840     "output_type": "execute_result",
1841     "data": {
1842      "text/plain": [
1843       "pi + 1.614213562373095"
1844      ]
1845     },
1847    }
1848   ],
1850  },
1851  {
1852   "execution_count": null,
1853   "cell_type": "code",
1854   "source": [
1855    "g.parent()"
1856   ],
1857   "outputs": [
1858    {
1859     "execution_count": 1,
1860     "output_type": "execute_result",
1861     "data": {
1862      "text/plain": [
1863       "Symbolic Ring"
1864      ]
1865     },
1867    }
1868   ],
1870  },
1871  {
1872   "source": [
1873    "But you might want to swim upstream. Given a floating point approximation of a real number, there are various ways to recover an exact representation it *may* come from.\n",
1874    "\n",
1875    "### From numerical to rational\n",
1876    "\n",
1877    "An element x of RealField(prec) can be transformed into a rational in different ways:\n",
1878    "\n",
1879    "-   x.exact_rational() returns the exact rational representation of the floating-point number x. In particular its denominator is a power of\n",
1880    "    1.  It contains the same information as .sign_mantissa_exponent :"
1881   ],
1882   "cell_type": "markdown",
1884  },
1885  {
1886   "execution_count": null,
1887   "cell_type": "code",
1888   "source": [
1889    "a = RR(pi) ; a"
1890   ],
1891   "outputs": [
1892    {
1893     "execution_count": 1,
1894     "output_type": "execute_result",
1895     "data": {
1896      "text/plain": [
1897       "3.14159265358979"
1898      ]
1899     },
1901    }
1902   ],
1904  },
1905  {
1906   "execution_count": null,
1907   "cell_type": "code",
1908   "source": [
1909    "a.sign_mantissa_exponent()"
1910   ],
1911   "outputs": [
1912    {
1913     "execution_count": 1,
1914     "output_type": "execute_result",
1915     "data": {
1916      "text/plain": [
1917       "(1, 7074237752028440, -51)"
1918      ]
1919     },
1921    }
1922   ],
1924  },
1925  {
1926   "execution_count": null,
1927   "cell_type": "code",
1928   "source": [
1929    "s,m,e = _\n",
1930    "s*m*2^(e)"
1931   ],
1932   "outputs": [
1933    {
1934     "execution_count": 1,
1935     "output_type": "execute_result",
1936     "data": {
1937      "text/plain": [
1938       "884279719003555/281474976710656"
1939      ]
1940     },
1942    }
1943   ],
1945  },
1946  {
1947   "execution_count": null,
1948   "cell_type": "code",
1949   "source": [
1950    "b = a.exact_rational()"
1951   ],
1952   "outputs": [
1953    {
1954     "execution_count": 1,
1955     "output_type": "execute_result",
1956     "data": {
1957      "text/plain": [
1958       "884279719003555/281474976710656"
1959      ]
1960     },
1962    }
1963   ],
1965  },
1966  {
1967   "source": [
1968    "-   x.simplest_rational() returns the simplest rational (i.e. with smallest denominator) that will become equal to x when converted into the parent of x :"
1969   ],
1970   "cell_type": "markdown",
1972  },
1973  {
1974   "execution_count": null,
1975   "cell_type": "code",
1976   "source": [
1977    "a = RR(pi) ; a"
1978   ],
1979   "outputs": [
1980    {
1981     "execution_count": 1,
1982     "output_type": "execute_result",
1983     "data": {
1984      "text/plain": [
1985       "3.14159265358979"
1986      ]
1987     },
1989    }
1990   ],
1992  },
1993  {
1994   "execution_count": null,
1995   "cell_type": "code",
1996   "source": [
1997    "b = a.simplest_rational() ; b"
1998   ],
1999   "outputs": [
2000    {
2001     "execution_count": 1,
2002     "output_type": "execute_result",
2003     "data": {
2004      "text/plain": [
2005       "245850922/78256779"
2006      ]
2007     },
2009    }
2010   ],
2012  },
2013  {
2014   "execution_count": null,
2015   "cell_type": "code",
2016   "source": [
2017    "a == b"
2018   ],
2019   "outputs": [
2020    {
2021     "execution_count": 1,
2022     "output_type": "execute_result",
2023     "data": {
2024      "text/plain": [
2025       "True"
2026      ]
2027     },
2029    }
2030   ],
2032  },
2033  {
2034   "execution_count": null,
2035   "cell_type": "code",
2036   "source": [
2037    "RR(b)"
2038   ],
2039   "outputs": [
2040    {
2041     "execution_count": 1,
2042     "output_type": "execute_result",
2043     "data": {
2044      "text/plain": [
2045       "3.14159265358979"
2046      ]
2047     },
2049    }
2050   ],
2052  },
2053  {
2054   "source": [
2055    "-   nearby_rational(max_error=e) finds the simplest rational that is at distance at most e of x :"
2056   ],
2057   "cell_type": "markdown",
2059  },
2060  {
2061   "execution_count": null,
2062   "cell_type": "code",
2063   "source": [
2064    "RR(pi).nearby_rational(max_error=0.01)"
2065   ],
2066   "outputs": [
2067    {
2068     "execution_count": 1,
2069     "output_type": "execute_result",
2070     "data": {
2071      "text/plain": [
2072       "22/7"
2073      ]
2074     },
2076    }
2077   ],
2079  },
2080  {
2081   "source": [
2082    "-   nearby_rational(max_denominator=d) finds the closest rational to x with denominator at most d :"
2083   ],
2084   "cell_type": "markdown",
2086  },
2087  {
2088   "execution_count": null,
2089   "cell_type": "code",
2090   "source": [
2091    "RR(pi).nearby_rational(max_denominator=120)"
2092   ],
2093   "outputs": [
2094    {
2095     "execution_count": 1,
2096     "output_type": "execute_result",
2097     "data": {
2098      "text/plain": [
2099       "355/113"
2100      ]
2101     },
2103    }
2104   ],
2106  },
2107  {
2108   "source": [
2109    "### From numerical to algebraic\n",
2110    "\n",
2111    "If x is an element of RealField(prec), x.algebraic_dependency(d) returns an integer polynomial P of degree at most d such that P(x) is almost zero:"
2112   ],
2113   "cell_type": "markdown",
2115  },
2116  {
2117   "execution_count": null,
2118   "cell_type": "code",
2119   "source": [
2120    "a = 3.0.nth_root(3)\n",
2121    "a += 5*a.ulp()\n",
2122    "a.algebraic_dependency(3)"
2123   ],
2124   "outputs": [
2125    {
2126     "execution_count": 1,
2127     "output_type": "execute_result",
2128     "data": {
2129      "text/plain": [
2130       "x^3 - 3"
2131      ]
2132     },
2134    }
2135   ],
2137  },
2138  {
2139   "source": [
2140    "### From algebraic field to (smaller, faster) embedded number fields\n",
2141    "\n",
2142    "Algebraic numbers have an as_number_field_element method:"
2143   ],
2144   "cell_type": "markdown",
2146  },
2147  {
2148   "execution_count": null,
2149   "cell_type": "code",
2150   "source": [
2151    "a = QQbar(2).sqrt() + QQbar(3).nth_root(3)\n",
2152    "a.as_number_field_element()"
2153   ],
2154   "outputs": [],
2156  },
2157  {
2158   "execution_count": null,
2159   "cell_type": "code",
2160   "source": [
2161    "(Number Field in a with defining polynomial y^6 - 6*y^4 + 6*y^3 + 12*y^2 + 36*y + 1,"
2162   ],
2163   "outputs": [
2164    {
2165     "execution_count": 1,
2166     "output_type": "execute_result",
2167     "data": {
2168      "text/plain": [
2169       " 96/755*a^5 - 54/755*a^4 - 128/151*a^3 + 936/755*a^2 + 1003/755*a + 2184/755,\n",
2170       " Ring morphism:\n",
2171       "   From: Number Field in a with defining polynomial y^6 - 6*y^4 + 6*y^3 + 12*y^2 + 36*y + 1\n",
2172       "   To:   Algebraic Real Field\n",
2173       "   Defn: a |--> -0.02803600793431334?)\n"
2174      ]
2175     },
2177    }
2178   ],
2180  },
2181  {
2182   "execution_count": null,
2183   "cell_type": "code",
2184   "source": [
2185    "field, number, morphism = _\n",
2186    "morphism(number) == a"
2187   ],
2188   "outputs": [
2189    {
2190     "execution_count": 1,
2191     "output_type": "execute_result",
2192     "data": {
2193      "text/plain": [
2194       "True"
2195      ]
2196     },
2198    }
2199   ],
2201  },
2202  {
2203   "source": [
2204    "For more than one algebraic number, the previous method generalizes with the number_field_elements_from_algebraics function:"
2205   ],
2206   "cell_type": "markdown",
2208  },
2209  {
2210   "execution_count": null,
2211   "cell_type": "code",
2212   "source": [
2213    "from sage.rings.qqbar import number_field_elements_from_algebraics\n",
2214    "sqrt2 = AA(2).sqrt()\n",
2215    "sqrt3 = AA(3).sqrt()\n",
2216    "number_field_elements_from_algebraics([sqrt2, sqrt3])"
2217   ],
2218   "outputs": [
2219    {
2220     "execution_count": 1,
2221     "output_type": "execute_result",
2222     "data": {
2223      "text/plain": [
2224       "(Number Field in a with defining polynomial y^4 - 4*y^2 + 1,\n",
2225       " [-a^3 + 3*a, -a^2 + 2],\n",
2226       " Ring morphism:\n",
2227       "   From: Number Field in a with defining polynomial y^4 - 4*y^2 + 1\n",
2228       "   To:   Algebraic Real Field\n",
2229       "   Defn: a |--> 0.5176380902050415?)"
2230      ]
2231     },
2233    }
2234   ],
2236  },
2237  {
2238   "execution_count": null,
2239   "cell_type": "code",
2240   "source": [
2241    "field, numbers, morphism = _\n",
2242    "morphism(numbers[0]) == sqrt2"
2243   ],
2244   "outputs": [
2245    {
2246     "execution_count": 1,
2247     "output_type": "execute_result",
2248     "data": {
2249      "text/plain": [
2250       "True"
2251      ]
2252     },
2254    }
2255   ],
2257  },
2258  {
2259   "source": [
2260    "### From algebraic to symbolic\n",
2261    "\n",
2262    "Algebraic numbers have a radical_expression method:"
2263   ],
2264   "cell_type": "markdown",
2266  },
2267  {
2268   "execution_count": null,
2269   "cell_type": "code",
2270   "source": [
2271    "x = polygen(QQ,'x')\n",
2272    "P = x^3-x^2-x-1\n",
2273    "P"
2274   ],
2275   "outputs": [
2276    {
2277     "execution_count": 1,
2278     "output_type": "execute_result",
2279     "data": {
2280      "text/plain": [
2281       "x^3 - x^2 - x - 1"
2282      ]
2283     },
2285    }
2286   ],
2288  },
2289  {
2290   "execution_count": null,
2291   "cell_type": "code",
2292   "source": [
2293    "P.roots(QQbar)[0][0]"
2294   ],
2295   "outputs": [
2296    {
2297     "execution_count": 1,
2298     "output_type": "execute_result",
2299     "data": {
2300      "text/plain": [
2301       "1.839286755214161?"
2302      ]
2303     },
2305    }
2306   ],
2308  },
2309  {
2310   "execution_count": null,
2311   "cell_type": "code",
2312   "source": [
2314   ],
2315   "outputs": [
2316    {
2317     "execution_count": 1,
2318     "output_type": "execute_result",
2319     "data": {
2320      "text/plain": [
2321       "(1/9*sqrt(11)*sqrt(3) + 19/27)^(1/3) + 4/9/(1/9*sqrt(11)*sqrt(3) +\n",
2322       "19/27)^(1/3) + 1/3"
2323      ]
2324     },
2326    }
2327   ],
2329  },
2330  {
2331   "source": [
2332    "> **note**\n",
2333    ">\n",
2334    "> Note that the radical_expression method does not implement the whole Galois theory yet, meaning that some algebraic numbers that can be represented by radicals, might not be modified that way:\n",
2335    ">\n",
2336    "> {.python .input}\n",
2337    "> a = QQbar(2).sqrt() + QQbar(3).nth_root(3)\n",
2339    "> \n",
2340    ">\n",
2341    "> {.json .output}\n",
2342    "> [{\"execution_count\": 1, \"output_type\": \"execute_result\", \"data\": {\"text/plain\": \"2.856463132680504?\"}, \"metadata\": {}}]\n",
2343    "> \n",
2344    "\n",
2345    "### Exercise: reverse engineering\n",
2346    "\n",
2347    "Suppose that, with a numerical experiment, you have approached an irrational algebraic number by:"
2348   ],
2349   "cell_type": "markdown",
2351  },
2352  {
2353   "execution_count": null,
2354   "cell_type": "code",
2355   "source": [
2356    "a = 3.14626436994197234"
2357   ],
2358   "outputs": [],
2360  },
2361  {
2362   "source": [
2363    "What is this algebraic number likely to be ?"
2364   ],
2365   "cell_type": "markdown",
2367  },
2368  {
2369   "execution_count": null,
2370   "cell_type": "code",
2371   "source": [
2372    "# edit here"
2373   ],
2374   "outputs": [],
2376  },
2377  {
2378   "source": [
2379    "## Complex numbers\n",
2380    "\n",
2381    "Most real number representation have their complex equivalent.\n",
2382    "\n",
2383    "-   RR becomes CC=ComplexField()\n",
2384    "-   RDF becomes CDF=ComplexDoubleField()\n",
2385    "-   RIF becomes CIF=ComplexIntervalField()\n",
2386    "-   RBF becomes CBF=ComplexBallField()\n",
2387    "-   QQ becomes QQ[I] (Gauss integers)\n",
2388    "-   AA becomes QQbar (the algebraic closure of QQ)\n",
2389    "-   SR remains SR (which contains much more than just complex numbers)\n",
2390    "\n",
2391    "Note that by default, I and i belong to the Symbolic Ring, so you will have to convert it, e.g. QQbar(I), CDF(I), ... or obtain it as QQbar.gen(), QQbar.gen(), ... . If you overwrite it, you can recover it with:"
2392   ],
2393   "cell_type": "markdown",
2395  },
2396  {
2397   "execution_count": null,
2398   "cell_type": "code",
2399   "source": [
2400    "from sage.symbolic.all import i as I"
2401   ],
2402   "outputs": [],
2404  },
2405  {
2406   "source": [
2407    "### Diamonds are not round\n",
2408    "\n",
2409    "We saw previously that contracting makes computations safe. Consider the following example:"
2410   ],
2411   "cell_type": "markdown",
2413  },
2414  {
2415   "execution_count": null,
2416   "cell_type": "code",
2417   "source": [
2418    "a = CBF(1/10+I*4/5)\n",
2419    "b = CBF((1+I)*2/3)\n",
2420    "a.diameter()"
2421   ],
2422   "outputs": [
2423    {
2424     "execution_count": 1,
2425     "output_type": "execute_result",
2426     "data": {
2427      "text/plain": [
2428       "4.4408921e-16"
2429      ]
2430     },
2432    }
2433   ],
2435  },
2436  {
2437   "execution_count": null,
2438   "cell_type": "code",
2439   "source": [
2440    "b.diameter()"
2441   ],
2442   "outputs": [
2443    {
2444     "execution_count": 1,
2445     "output_type": "execute_result",
2446     "data": {
2447      "text/plain": [
2448       "4.4408921e-16"
2449      ]
2450     },
2452    }
2453   ],
2455  },
2456  {
2457   "execution_count": null,
2458   "cell_type": "code",
2459   "source": [
2460    "b.abs() < 1"
2461   ],
2462   "outputs": [
2463    {
2464     "execution_count": 1,
2465     "output_type": "execute_result",
2466     "data": {
2467      "text/plain": [
2468       "True"
2469      ]
2470     },
2472    }
2473   ],
2475  },
2476  {
2477   "execution_count": null,
2478   "cell_type": "code",
2479   "source": [
2480    "(a*b).diameter()"
2481   ],
2482   "outputs": [
2483    {
2484     "execution_count": 1,
2485     "output_type": "execute_result",
2486     "data": {
2487      "text/plain": [
2488       "1.1768364e-15"
2489      ]
2490     },
2492    }
2493   ],
2495  },
2496  {
2497   "source": [
2498    "Explain this loss of precision.\n",
2499    "\n",
2500    "*Hint*: draw a picture !\n",
2501    "\n",
2502    "Complex ball elements are pairs of real ball elements, see <http://arblib.org/acb.html>\n",
2503    "\n",
2504    "## Transitional representations (under the hood)\n",
2505    "\n",
2506    "Some \"transitional representations\" were not included in the big zoo for simplicity, and because they exist mostly as wrappers for internal stuff.\n",
2507    "\n",
2508    "### Real literals: between your fingers and Sage\n",
2509    "\n",
2510    "When the user writes:"
2511   ],
2512   "cell_type": "markdown",
2514  },
2515  {
2516   "execution_count": null,
2517   "cell_type": "code",
2518   "source": [
2519    "a = 0.1"
2520   ],
2521   "outputs": [],
2523  },
2524  {
2525   "source": [
2526    "this defines a floating-point number though the precision is not explicitely set.\n",
2527    "\n",
2528    "Sage defines some default precision (depending on the number of digits that are provided, 53 in the current case):"
2529   ],
2530   "cell_type": "markdown",
2532  },
2533  {
2534   "execution_count": null,
2535   "cell_type": "code",
2536   "source": [
2537    "type(a)"
2538   ],
2539   "outputs": [
2540    {
2541     "execution_count": 1,
2542     "output_type": "execute_result",
2543     "data": {
2544      "text/plain": [
2545       "<type 'sage.rings.real_mpfr.RealLiteral'>"
2546      ]
2547     },
2549    }
2550   ],
2552  },
2553  {
2554   "execution_count": null,
2555   "cell_type": "code",
2556   "source": [
2557    "a.parent()"
2558   ],
2559   "outputs": [
2560    {
2561     "execution_count": 1,
2562     "output_type": "execute_result",
2563     "data": {
2564      "text/plain": [
2565       "Real Field with 53 bits of precision"
2566      ]
2567     },
2569    }
2570   ],
2572  },
2573  {
2574   "source": [
2575    "but it allows casting into higher precision rings, without losing the precision from the default choice:"
2576   ],
2577   "cell_type": "markdown",
2579  },
2580  {
2581   "execution_count": null,
2582   "cell_type": "code",
2583   "source": [
2584    "RealField(200)(0.1)"
2585   ],
2586   "outputs": [
2587    {
2588     "execution_count": 1,
2589     "output_type": "execute_result",
2590     "data": {
2591      "text/plain": [
2592       "0.10000000000000000000000000000000000000000000000000000000000"
2593      ]
2594     },
2596    }
2597   ],
2599  },
2600  {
2601   "execution_count": null,
2602   "cell_type": "code",
2603   "source": [
2604    "RealField(200)(RDF(0.1))"
2605   ],
2606   "outputs": [
2607    {
2608     "execution_count": 1,
2609     "output_type": "execute_result",
2610     "data": {
2611      "text/plain": [
2612       "0.10000000000000000555111512312578270211815834045410156250000"
2613      ]
2614     },
2616    }
2617   ],
2619  },
2620  {
2621   "source": [
2622    "### Real Lazy Field: between exact and numerical representations\n",
2623    "\n",
2624    "Elements of RLF = RealLazyField() behaves like floating-point numbers, This represen\n",
2625    "\n",
2626    "while\n",
2627    "\n",
2628    "RLF = RealLazyField() is lazy, in the sense that it doesn't really do anything but simply sits between exact rings of characteristic 0 and the real numbers. The values are actually computed when they are cast into a field of fixed precision.\n",
2629    "\n",
2630    "RLF = RealLazyField() like RealField() but can wait that the user fixes the precision to get one. <http://www.sagemath.org/doc/reference/rings_numerical/sage/rings/real_lazy.html>"
2631   ],
2632   "cell_type": "markdown",
2634  },
2635  {
2636   "execution_count": null,
2637   "cell_type": "code",
2638   "source": [
2639    "a = RLF(pi) + 2\n",
2640    "a"
2641   ],
2642   "outputs": [
2643    {
2644     "execution_count": 1,
2645     "output_type": "execute_result",
2646     "data": {
2647      "text/plain": [
2648       "5.141592653589794?"
2649      ]
2650     },
2652    }
2653   ],
2655  },
2656  {
2657   "execution_count": null,
2658   "cell_type": "code",
2659   "source": [
2660    "RealField(100)(a)"
2661   ],
2662   "outputs": [
2663    {
2664     "execution_count": 1,
2665     "output_type": "execute_result",
2666     "data": {
2667      "text/plain": [
2668       "5.1415926535897932384626433833"
2669      ]
2670     },
2672    }
2673   ],
2675  },
2676  {
2677   "execution_count": null,
2678   "cell_type": "code",
2679   "source": [
2680    "RealField(150)(a)"
2681   ],
2682   "outputs": [
2683    {
2684     "execution_count": 1,
2685     "output_type": "execute_result",
2686     "data": {
2687      "text/plain": [
2688       "5.1415926535897932384626433832795028841971694"
2689      ]
2690     },
2692    }
2693   ],
2695  },
2696  {
2697   "source": [
2698    "The precision of a can be increased on demand since, internally, a keeps the exact representation it comes from:"
2699   ],
2700   "cell_type": "markdown",
2702  },
2703  {
2704   "execution_count": null,
2705   "cell_type": "code",
2706   "source": [
2707    "a._value"
2708   ],
2709   "outputs": [
2710    {
2711     "execution_count": 1,
2712     "output_type": "execute_result",
2713     "data": {
2714      "text/plain": [
2715       "pi + 2"
2716      ]
2717     },
2719    }
2720   ],
2722  },
2723  {
2724   "source": [
2725    "Its use is sometimes tricky, and shares some weirdnesses of the symbolic ring since it also warps its elements:"
2726   ],
2727   "cell_type": "markdown",
2729  },
2730  {
2731   "execution_count": null,
2732   "cell_type": "code",
2733   "source": [
2734    "RLF(pi) + RLF(QQbar(2).sqrt())"
2735   ],
2736   "outputs": [
2737    {
2738     "execution_count": 1,
2739     "output_type": "execute_result",
2740     "data": {
2741      "text/plain": [
2742       "4.555806215962889?"
2743      ]
2744     },
2746    }
2747   ],
2749  },
2750  {
2751   "execution_count": null,
2752   "cell_type": "code",
2753   "source": [
2754    "_._value"
2755   ],
2756   "outputs": [
2757    {
2758     "execution_count": 1,
2759     "output_type": "execute_result",
2760     "data": {
2761      "text/plain": [
2762       "pi + 1.414213562373095?"
2763      ]
2764     },