2f:T4cf,# Real numbers with `BigInt`, and the Arithmetic-Geometric Mean (AGM) Following Alan Turing, we define a (computable) real number as an algorithm: given a (rational) *precision* ${tex`\epsilon > 0`}, it outputs two *rational* numbers, the lower and upper approximations (or bounds) of the "number" that are less than ${tex`\epsilon`} apart. Two algorithms would define the same number if they "converge" to the same limit as ${tex`\epsilon\to 0`}; or rather, two algorithms are *different* if, for some precision, the upper bound of one is less than the lower bound of the other. Much of the subject of *calculus* (of differentials and integrals) can be regarded as the art and science of deciding if two algorithms give the same number. Since your browser now supports `BigInt`, why not implement rational numbers as pairs of `BigInt`s and calculate real numbers to any precision? For example, to compute the sum of two real numbers, for a given precision ${tex`\epsilon`}, is to compute each number with precision ${tex`\epsilon/2`}, and add the results. Our goal is to compute ${tex`\pi`} using the arithmetic-geometric mean (AGM) algorithm. Let's proceed step by step. We need first to define a `class` for rational numbers.30:T55c,Different types of real numbers would naturally have different types of algorithms, some of which are given below. Note that we are not concerned with the complexity of the algorithm, nor how "good" (with smaller denominator) the outputs are, so long as it does the job. - algebraic numbers, solutions to polynomial equations, e.g., ${tex`\sqrt 2`} (Newton's method works great). - infinite series, which supposedly cover all numbers (by definition): ${tex.block`e=1+\frac{1}{1!}+\frac{1}{2!} + \frac{1}{3!} + \cdots `} ${tex.block`\sin 1 = 1 - \frac{1}{3!} + \frac{1}{5!} -\cdots`} (Easy to implement, though you may need to take extra care of the "error estimates".) - periods, i.e., integrals of algebraic functions. ${tex.block`\frac{\pi}{2}=\int_0^1\frac{dx}{\sqrt{1-x^2}}`} (The upper and lower Riemann sums are natural fit for this.) It should be noted that rational numbers also fit here: the algorithm simply outputs that number (for both bounds), no matter the precision. That is, one may simply "store" the two integers in the algorithm. (If the integers are too big, such as ${tex`100!`}, it may be shorter to write an algorithm to compute them. Let's not worry about numbers too big even for `BigInt`.) To get a feel for this, here's one algorithm for ${tex`\sqrt d`}. For easier reading, it does not actually use `QQ` (so we get decimals for the bounds).31:T417,## Arithmetic-geometric mean Start with two positive numbers ${tex`a_0 < b_0`}, we can calculate iteratively the arithmetic and geometric means of the previous pair: ${tex.block` a_{n+1} := \sqrt{a_{n} b_{n}} \qquad b_{n+1} := \frac{a_{n}+b_{n}}{2}`} The two sequences ${tex`\{a_n\}`} and ${tex`\{b_n\}`} converge to the same number, known as the **arithmetic-geometric mean** of ${tex`a_0`} and ${tex`b_0`}, denoted ${tex`M(a_0, b_0)`}. For example, Gauss found that ${tex.block`M(1,\sqrt 2) = \frac{\pi}{\varpi}`} where ${tex`\varpi`} is the Lemniscate constant. The rate of convergence is quadratic, meaning that you double the number of correct digits with each iteration. In 1975, a formula for ${tex`\pi`} was discovered that "launched the modern computer era of the computation of ${tex`\pi`}" ${tex.block`\pi=\frac{M(1/\sqrt2, 1)^2}{\frac{1}{4}-\sum_{n=1}^\infty 2^{n-1} c_n^2}`} where ${tex`c_n := (b_{n-1}-a_{n-1})/2 = b_{n-1}-b_n`} from the iterations. This gives a pretty good algorithm for computing ${tex`\pi`}. Let's implement it.32:T45c,class QQ { constructor(p, q=1, reduced=false) { if (q==0) {p=0; q=1} if (q==Infinity || p==Infinity) {p=0; q=1} if (q==Infinity) {p=0; q=1} if (typeof(p)=="number" && typeof(q)=="number") { while (!Number.isInteger(p) || !Number.isInteger(q)) {p*=10; q*=10} } this.p=BigInt(p) this.q=BigInt(q) if (q<0) {this.p*=-1n; this.q*=-1n} if (!reduced) this.reduce() this.value=Number(this.p)/Number(this.q) this.tex=(this.p>=0n)?"":"-" this.tex+=(this.q==1n)?ZZ.abs(this.p):"\\dfrac{"+ZZ.abs(this.p)+"}{"+this.q+"}" } reduce() { let gcd=ZZ.gcd(this.p, this.q) if (gcd>1n) {this.p/=gcd; this.q/=gcd} } multiply(x) {return new QQ(this.p*x.p, this.q*x.q)} inverse() {return new QQ(this.q, this.p, false)} divide(x) {return new QQ(this.p*x.q, this.q*x.p)} add(x) {return new QQ(this.p*x.q+this.q*x.p, this.q*x.q)} subtract(x) {return new QQ(this.p*x.q-this.q*x.p, this.q*x.q)} less_than(x) {return this.p*x.qe) {\n a=(d+a*b)/(a+b)\n b=(b+d/b)/2\n n++\n }\n return [a,b,n]\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":444,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":534,"value":"sqrt_split(0.001, 2, 2.0012)","pinned":true,"hidden":false,"mode":"ojs"},{"id":770,"value":"{\n let e=new QQ(1n,10n**30n);\n let two=new QQ(2n);\n return sqrt_split(e, two, two.add(e));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":499,"value":"sqrt = function(d, e) {\n e = (typeof(e)==\"number\")?new QQ(e):e\n // d1 and d2 are now QQ\n d=(typeof(d)==\"number\")?new QQ(d):d\n let b=new QQ(0)\n while (b.pow(2).less_than(d2)) b=b.add(new QQ(1))\n let a=b.subtract(new QQ(1))\n let steps=0\n while (e.less_than(b.subtract(a))) {\n a=(d1.add(a.multiply(b))).divide(a.add(b))\n b=b.add(d2.divide(b)).divide(new QQ(2))\n steps++\n if (b.less_than(a)) break;\n }\n return [a,b,steps]\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":497,"value":"AGM = function(a, b, e) {\n e = (typeof(e)==\"number\")?new QQ(e):e\n let [a1,a2,steps]=sqrt_split(e, new QQ(1,2));\n let b1=new QQ(1), b2=new QQ(1);\n let s1=new QQ(1,4), s2=new QQ(1,4);\n let b1_=b1, b2_=b2, steps_=0;\n let r1=a1.pow(2).divide(s2), r2=b2.pow(2).divide(s1);\n while (steps<10 && e.less_than(r2.subtract(r1))) {\n b1_=a1.add(b1).divide(new QQ(2));\n b2_=a2.add(b2).divide(new QQ(2));\n [a1,a2,steps_]=sqrt_split(e, a1.multiply(b1), a2.multiply(b2))\n s1=s1.subtract(b2.subtract(b1_).pow(2))\n s2=s2.subtract(b1.subtract(b2_).pow(2))\n b1=b1_, b2=b2_\n steps+=steps_\n r1=a1.pow(2).divide(s2), r2=b2.pow(2).divide(s1)\n if (r2.less_than(r1)) return [NaN, NaN, steps];\n }\n return [r1, r2, steps]\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":421,"value":"## Code","pinned":false,"hidden":false,"mode":"md"},{"id":11,"value":"ZZ = {\n // class for operations on BigInts\n const abs=(n)=>(n>=0)?n:-n\n const sign=(n)=>(n>=0)?\"+\":\"-\"\n const min=(a,b)=>(a(a>b)?a:b\n \n const ht=(list)=>{\n let m=list[0], M=list[0], n=0\n for (let i=1; iabs(M)) M=list[i]\n }\n return [m, M]\n }\n \n const gcd=(a, b)=>{\n if (abs(a){\n n=BigInt(n)\n if (n<0n) return Infinity\n if (n==0n) return 1n\n return factorial(n-1n)*n\n }\n \n return {abs, sign, gcd, min, max, ht, factorial}\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}