31:T1064,class QQ { constructor(p, q = 1n, reduced = false) { if (q === 0n) throw new Error(`Denominator cannot be zero: ${p}/${q}`); this.p = q < 0n ? -BigInt(p) : BigInt(p); this.q = q < 0n ? -BigInt(q) : BigInt(q); if (!reduced) this.simplify(); this.value = Number(this.p) / Number(this.q); this.cf = null; this.cfrac = this.continuedFraction(); } static toQQ(value) { if (value instanceof QQ) return value; if (typeof value === 'bigint') return new QQ(value, 1n); if (typeof value === 'number' && Number.isInteger(value)) { return new QQ(BigInt(value), 1n); } throw new Error(`Invalid input type (${typeof value}). Expected QQ, bigint, or integer number.`); } simplify() { const gcd = (a, b) => { while (b !== 0n) [a, b] = [b, a % b]; return a; }; const divisor = gcd(this.p < 0n ? -this.p : this.p, this.q); this.p /= divisor; this.q /= divisor; } sign() { return this.p < 0n ? -1n : 1n } add(other) { other = QQ.toQQ(other); return new QQ(this.p * other.q + other.p * this.q, this.q * other.q); } subtract(other) { other = QQ.toQQ(other); return new QQ(this.p * other.q - other.p * this.q, this.q * other.q); } multiply(other) { other = QQ.toQQ(other); return new QQ(this.p * other.p, this.q * other.q); } divide(other) { return this.multiply(QQ.toQQ(other).inverse()); } inverse() { return new QQ(this.q, this.p, true); } power(k) { if (k === 0n) return new QQ(1n, 1n); // Any number to the power of 0 is 1 if (k > 0n) return new QQ(this.p ** k, this.q ** k, true); // For negative k, invert the fraction and use the positive exponent return new QQ(this.q ** -k, this.p ** -k, true); } lessThan(other) { other = QQ.toQQ(other); return this.p * other.q < other.p * this.q; } equalTo(other) { other = QQ.toQQ(other); return this.p * other.q === other.p * this.q; } greaterThan(other) { other = QQ.toQQ(other); return this.p * other.q > other.p * this.q; } continuedFraction() { let p = this.p < 0n ? this.p % this.q + this.q : this.p; let q = this.q; const cfrac = []; while (q !== 0n) { cfrac.push(p / q); [p, q] = [q, p % q]; } if (this.p < 0n) cfrac[0] = this.p / this.q - 1n; this.cf = `[${cfrac[0]};${cfrac.slice(1).join(',')}]`; this.cfrac = cfrac; return cfrac; } static from_cfrac(cfrac) { if (cfrac.length === 0) throw new Error('Continued fraction list cannot be empty'); let p = 1n, q = 0n; // Initialize p and q for backward calculation for (let i = cfrac.length - 1; i >= 0; i--) { [p, q] = [BigInt(cfrac[i]) * p + q, p]; // Correct order of swapping } return new QQ(p, q); // Return the rational number QQ(p, q) } approximant(max, direction = null) { if (this.q <= max) { return new QQ(this.p, this.q, true); } const a = this.continuedFraction(); const c = []; let pM2 = 0n, pM1 = 1n, qM2 = 1n, qM1 = 0n; for (const ai of a) { const steps = ai < 0n ? -ai : ai; for (let m = 1n; m <= steps; m++) { const p = m * pM1 + pM2, q = m * qM1 + qM2; if (q > max) break; c.push([p, q]); } const p0 = ai * pM1 + pM2, q0 = ai * qM1 + qM2; pM2 = pM1; qM2 = qM1; pM1 = p0; qM1 = q0; } if (!c.length) { return new QQ(0n, 1n); } let [bestN, bestD] = c[0], abs = (x) => x < 0n ? -x : x; let bestErr = abs(this.p * bestD - this.q * bestN); for (let i = 1; i < c.length; i++) { const [n, d] = c[i], diff = this.p * d - this.q * n; if (direction === 'down' && diff < 0n) continue; if (direction === 'up' && diff > 0n) continue; const err = abs(diff); if (err < bestErr || (err === bestErr && d < bestD)) { bestErr = err; bestN = n; bestD = d; } } return new QQ(bestN, bestD, true); } relax(max) { const bound = max < 0n ? -max : max; if (bound >= this.q) return this; const a = (this.p * bound) / this.q; if (max * this.p < 0n) return new QQ(a, bound); return new QQ(a + (this.p > 0n ? 1n : -1n), bound); } }32:T1295,// RR Class: Computable Real Numbers as Approximation Algorithms class RR { constructor(approx) { this.value = null; this.cf = ''; this.cfrac = []; this.approx = (N) => { if (N < this.N_min) N = this.N_min; let [a, b] = approx(N*2n); const delta = new QQ(1n, N).subtract(b.subtract(a)); if (delta.p > 0n) { const max = (delta.q * 2n) / delta.p; a = a.relax(-max); // Relax a slightly lower b = b.relax(max); // Relax b slightly higher } let cfrac = [], i = 0; while (i < a.cfrac.length && i < b.cfrac.length) { if (a.cfrac[i] != b.cfrac[i]) break; cfrac.push(a.cfrac[i]); i++; } if (cfrac.length > this.cfrac.length) { this.cf = `[${cfrac[0]};${cfrac.slice(1).join(',')},...]`; this.cfrac = cfrac; } return [a, b]; }; this.N_min = 1n; this.sign = this.sign(); } // Static method to create an RR instance from a QQ instance static toRR(q) { q = QQ.toQQ(q); return new RR((N) => [q, q]); } sign() { let N = 1n; while (true) { const [a, b] = this.approx(N); this.N_min = N; if (a.p > 0) return 1n; if (b.p < 0) return -1n; N *= 2n; } } add(other) { return new RR((N) => { if (other instanceof QQ || typeof other === 'bigint') { // Approximate this with precision N and add the QQ directly const [a1, b1] = this.approx(N); return [a1.add(other), b1.add(other)]; } else if (other instanceof RR) { // Use 2N to approximate both numbers and add them const [a1, b1] = this.approx(N * 2n); const [a2, b2] = other.approx(N * 2n); return [a1.add(a2), b1.add(b2)]; } else { throw new TypeError("Unsupported type for addition"); } }); } subtract(other) { return new RR((N) => { if (other instanceof QQ || typeof other === 'bigint') { // Approximate this with precision N and add the QQ directly const [a1, b1] = this.approx(N); return [a1.subtract(other), b1.subtract(other)]; } else if (other instanceof RR) { // Use 2N to approximate both numbers and add them const [a1, b1] = this.approx(N * 2n); const [a2, b2] = other.approx(N * 2n); return [a1.subtract(b2), b1.subtract(a2)]; } else { throw new TypeError("Unsupported type for addition"); } }); } multiply(other) { return new RR((N) => { if (other instanceof QQ || typeof other === 'bigint') { // Multiply with a QQ directly other = QQ.toQQ(other); if (other.p >= 0n) { const [a, b] = this.approx(N * other.p / other.q + 1n); return [a.multiply(other), b.multiply(other)]; } else { const [a, b] = this.approx(N * (-other.p) / other.q + 1n); return [b.multiply(other), a.multiply(other)]; } } else if (other instanceof RR) { const eps = new QQ(1n, N); while (true) { const [a1, b1] = this.approx(N); const [a2, b2] = other.approx(N); let [a, b] = [a1.multiply(a2), b1.multiply(b2)]; if (this.sign * other.sign === -1n) [a, b] = [b, a]; if (b.subtract(a).lessThan(eps)) return [a, b]; N *= 2n; } } else { throw new TypeError("Unsupported type for multiplication"); } }); } divide(other) { if (other instanceof QQ || typeof other === 'bigint') { return this.multiply(QQ.toQQ(other).inverse()); } else if (other instanceof RR) { return this.multiply(other.inverse()); } else { throw new TypeError(`Unsupported type for division (${typeof other})`); } } power(k) { k = BigInt(k); if (k === 0n) return new QQ(1n); return new RR((N) => { const eps = new QQ(1n, N); const est = this.approx(1n)[1]; while (true) { const est_k = est.power(k-1n) let N_prime = N * (k * est_k.p) / est_k.q; if (N_prime < 0n) N_prime = -N_prime; let [a, b] = this.approx(N_prime); let [a_k, b_k] = [a.power(k), b.power(k)]; a_k = a_k.relax(-N*4n); b_k = b_k.relax(N*4n); if (a_k.greaterThan(b_k)) [a_k, b_k] = [b_k, a_k]; if (b_k.subtract(a_k).lessThan(eps)) return [a_k, b_k]; N *= 2n; } }); } inverse() { return new RR((N) => { const eps = new QQ(1n, N); const sgn = this.sign(); while (true) { let [a, b] = this.approx(N); if (a.p * b.p > 0n) { [a, b] = [b.inverse(), a.inverse()]; if (b.subtract(a).lessThan(eps)) return [a, b]; } N *= 2n; } }); } }33:T482,class Alternating extends RR { constructor(terms, start = 0n) { // Approximation function for alternating series // t_0 - t_1 + t_2 - t_3 + ... // where the terms are given by terms(), starting from 0 by default super((N) => { const max = (x, y) => (x > y ? x : y); let eps = new QQ(1n, N), K = 100n; while (terms(K).greaterThan(eps)) K = K * 2n; let k = BigInt(start); // Current term index let b = terms(k); // Start with the first term let a = b.subtract(terms(k+1n)); // Initial lower bound k = k+2n; let sign = 1n; let cache = [terms(k), terms(k+1n), terms(k+2n)]; while (b.subtract(a).greaterThan(eps)) { cache.push(terms(k+3n)); const delta = cache[2].subtract(cache[3]); const bound = max(delta.q / delta.p, K*K); if (sign < 0n) { a = b.subtract(cache[0]).relax(-bound); // Update the lower bound } else { b = a.add(cache[0]).relax(bound); // Update the upper bound } cache.shift(); k++; sign = -sign; // Flip sign for next iteration } return [a, b]; }); } }34:[] 2f:{"title":"Rational, Real and ppp-adic numbers","cells":[{"id":148,"value":"# Rational, Real and ${tex`p`}-adic numbers\n\n* Rational numbers as pairs of `BigInt`s\n* (Computable) real numbers as algorithms that compute the lower and upper *rational* bounds up to a given precision\n* TODO *p*-adic numbers","pinned":false,"hidden":false,"mode":"md"},{"id":1915,"value":"## Rational numbers","pinned":false,"hidden":false,"mode":"md"},{"id":29,"value":"x = new QQ(355n, 113n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1693,"value":"x.power(2n).divide(6n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":5,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":1622,"value":"new QQ(99n,200n).relax(-100n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1917,"value":"## Real numbers","pinned":false,"hidden":false,"mode":"md"},{"id":34,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":2138,"value":"sqrt(4n).power(5n).approx(100n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":357,"value":"function sqrt(d) {\n if (d instanceof RR) {\n if (d.sign == -1n) throw new Error(`Can not take sqrt of negative number: ${d.cf}`);\n return new RR((N)=> {\n const eps = new QQ(1n, N);\n while (true) {\n let [a, b] = d.approx(N);\n if (a.p > 0n) {\n [a, b] = [sqrt(a).approx(N)[0], sqrt(b).approx(N)[1]];\n if (b.subtract(a).lessThan(eps)) return [a, b];\n }\n N *= 2n;\n }\n });\n } else {\n d = QQ.toQQ(d);\n if (d.p < 0n) throw new Error(`Can not take sqrt of negative number: ${d.p}/${d.q}`);;\n return new RR((N) => {\n const eps = new QQ(1n, N);\n let a = new QQ(0n), b = d.add(1n).divide(2n);\n while (b.subtract(a).greaterThan(eps)) {\n a = d.add(a.multiply(b)).divide(a.add(b));\n b = b.add(d.divide(b)).divide(2n);\n a = a.relax(-N*4n); // Relax a slightly lower\n b = b.relax(N*4n); // Relax b slightly higher\n }\n return [a, b];\n });\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":361,"value":"sqrt(6n).approx(1000n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2033,"value":"RR.toRR(2n).multiply(sqrt(3n)).approx(100n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1389,"value":"sqrt(5n).add(1n).divide(2n).approx(10000n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":2046,"value":"sqrt(sqrt(2n)).approx(1000n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1919,"value":"### Alternating series","pinned":false,"hidden":false,"mode":"md"},{"id":1415,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":1417,"value":"function arctan(x) {\n x = QQ.toQQ(x);\n if (x.p < 0n) return arctan(x.multiply(-1n)).multiply(-1n);\n if (x.greaterThan(1n)) {\n const piOver2 = arctan(new QQ(1n,5n)).multiply(4n).subtract(arctan(new QQ(1n,239n))).multiply(2n)\n return piOver2.subtract(arctan(x.inverse()));\n }\n return new Alternating((k) => x.power(2n*k+1n).divide(2n*k+1n));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1574,"value":"pi = arctan(new QQ(1n,5n)).multiply(4n).subtract(arctan(new QQ(1n,239n))).multiply(4n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1578,"value":"pi.approx(1000000000n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1474,"value":"new RR((N) => {\n let sum = new QQ(0n);\n for (let k=1; k x.power(2n).add(1n).inverse(), 0n, -1n);","pinned":true,"hidden":false,"mode":"ojs"},{"id":1566,"value":"arctan_integral.multiply(4n).approx(100n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":95,"value":"// Subclass for computing the integral of a monotonic rational function\nclass MonoIntegral extends RR {\n constructor(f, x1, x2) {\n super((N) => {\n x1 = QQ.toQQ(x1), x2 = QQ.toQQ(x2);\n const eps = new QQ(1n, N);\n let k = 1n, dx = x2.subtract(x1);\n let sum1 = f(x1), sum2 = f(x2);\n if (sum1.multiply(dx).greaterThan(sum2.multiply(dx))) {\n [sum1, sum2] = [sum2, sum1];\n }\n while (true) {\n let xi = x1.add(dx);\n for (let i = 1n; i < k; i += 2n) {\n sum1 = sum1.add(f(xi)).relax(-N*N);\n sum2 = sum2.add(f(xi)).relax(N*N);\n xi = xi.add(dx).add(dx);\n }\n const a = sum1.multiply(dx).relax(-N*4n);\n const b = sum2.multiply(dx).relax(N*4n);\n if (b.subtract(a).lessThan(eps)) return [a, b];\n k *= 2n;\n dx = dx.divide(2n);\n }\n });\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2105,"value":"# ${tex`p`}-adic numbers\n\nOne might think that we simply replace the condition `b.subtract(a).lessThan(eps)` with the ${tex`p`}-adic norm on ${tex`\\mathbb{Q}`}, we'd get the ${tex`p`}-adic version of all the familiar real numbers.","pinned":false,"hidden":false,"mode":"md"}],"files":"$W34","stdlib":"1","titling":"auto"}