From 22b15a3c6de45bb6ee45ea278b7923cb135291f1 Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Tue, 26 Dec 2023 21:17:34 +0200 Subject: [PATCH 1/2] Make prime limit approximator more sane --- src/__tests__/approximation.spec.ts | 6 +++ src/approximation.ts | 63 +++++++++++++++++------------ 2 files changed, 44 insertions(+), 25 deletions(-) diff --git a/src/__tests__/approximation.spec.ts b/src/__tests__/approximation.spec.ts index 23e41f7..1d4669f 100644 --- a/src/__tests__/approximation.spec.ts +++ b/src/__tests__/approximation.spec.ts @@ -62,6 +62,12 @@ describe('Prime limit approximator', () => { expect(calculatedError).toBeLessThanOrEqual(10); }); }); + + it('has somewhat sane default behavior', () => { + expect(() => + approximatePrimeLimit(valueToCents(Math.PI), 8, 2) + ).not.toThrow(); + }); }); describe('Convergent calculator', () => { diff --git a/src/approximation.ts b/src/approximation.ts index 169640e..b37970e 100644 --- a/src/approximation.ts +++ b/src/approximation.ts @@ -168,16 +168,19 @@ export function approximatePrimeLimitWithErrors( limitIndex: number, maxExponent: number, maxError = 600, - maxLength?: number + maxLength = 100 ) { if (maxError > 600) { throw new Error('Maximum search distance is 600 cents'); } const results: [Fraction, number][] = []; - function push(error: number, generator: () => Fraction) { - if (maxLength === undefined || results.length < maxLength) { - results.push([generator(), error]); + function push(error: number, result: Fraction) { + if (!result.n) { + return; + } + if (results.length < maxLength) { + results.push([result, error]); if (results.length === maxLength) { results.sort((a, b) => a[1] - b[1]); } @@ -187,7 +190,7 @@ export function approximatePrimeLimitWithErrors( } for (let index = results.length - 1; ; index--) { if (error <= results[index][1]) { - results.splice(index, 0, [generator(), error]); + results.splice(index, 0, [result, error]); break; } } @@ -200,6 +203,9 @@ export function approximatePrimeLimitWithErrors( approximationCents: number, index: number ) { + if (approximation.n > 10e10 || approximation.d > 10e10) { + return; + } if (index > limitIndex) { // Procedure is the same as in approximateOddLimitWithErrors const remainder = mmod(cents - approximationCents, 1200); @@ -208,34 +214,41 @@ export function approximatePrimeLimitWithErrors( if (error > maxError) { return; } - push(error, () => - approximation.mul( - TWO.pow( - Math.round((cents - approximationCents - remainder) / 1200) - )! - ) + const exponent = Math.round( + (cents - approximationCents - remainder) / 1200 ); + if (Math.abs(exponent) > 52) { + return; + } + const result = approximation.mul(TWO.pow(exponent)!); + push(error, result); } else { const error = 1200 - remainder; if (error > maxError) { return; } - push(error, () => - approximation.mul( - TWO.pow( - Math.round((cents - approximationCents - remainder) / 1200) + 1 - )! - ) - ); + const exponent = + Math.round((cents - approximationCents - remainder) / 1200) + 1; + if (Math.abs(exponent) > 52) { + return; + } + const result = approximation.mul(TWO.pow(exponent)!); + push(error, result); } return; } - approximation = approximation.div(PRIMES[index] ** maxExponent); - approximationCents -= PRIME_CENTS[index] * maxExponent; - for (let i = -maxExponent; i <= maxExponent; ++i) { - accumulate(approximation, approximationCents, index + 1); - approximation = approximation.mul(PRIMES[index]); - approximationCents += PRIME_CENTS[index]; + accumulate(approximation, approximationCents, index + 1); + for (let i = 1; i <= maxExponent; ++i) { + accumulate( + approximation.mul(PRIMES[index] ** i), + approximationCents + PRIME_CENTS[index] * i, + index + 1 + ); + accumulate( + approximation.div(PRIMES[index] ** i), + approximationCents - PRIME_CENTS[index] * i, + index + 1 + ); } } accumulate(new Fraction(1), 0, 1); @@ -259,7 +272,7 @@ export function approximatePrimeLimit( limitIndex: number, maxExponent: number, maxError = 600, - maxLength?: number + maxLength = 100 ) { return approximatePrimeLimitWithErrors( cents, From e6e9e2cdce04a08ada4881d53b634fbaf8fa1f0c Mon Sep 17 00:00:00 2001 From: Lumi Pakkanen Date: Tue, 26 Dec 2023 21:37:48 +0200 Subject: [PATCH 2/2] Make convergent calculator more sane --- src/__tests__/approximation.spec.ts | 6 ++++++ src/approximation.ts | 22 +++++++++++++--------- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/src/__tests__/approximation.spec.ts b/src/__tests__/approximation.spec.ts index 1d4669f..c3e6a4c 100644 --- a/src/__tests__/approximation.spec.ts +++ b/src/__tests__/approximation.spec.ts @@ -120,6 +120,12 @@ describe('Convergent calculator', () => { expect(semiconvergents[3].equals('10/3')).toBeTruthy(); expect(semiconvergents[4].equals('13/4')).toBeTruthy(); }); + + it('calculates convergents for 1\\5', () => { + expect(() => + getConvergents(1.148698354997035, undefined, 256, true, false) + ).not.toThrow(); + }); }); describe('Radical approximator', () => { diff --git a/src/approximation.ts b/src/approximation.ts index b37970e..a9737f4 100644 --- a/src/approximation.ts +++ b/src/approximation.ts @@ -52,17 +52,21 @@ export function getConvergents( if (includeNonMonotonic) { result.push(convergent); } else { - // See https://en.wikipedia.org/wiki/Continued_fraction#Semiconvergents - // for the origin of this half-rule if (2 * i > cfDigit) { result.push(convergent); - } else if ( - convergent - .sub(value_) - .abs() - .compare(result[result.length - 1].sub(value_).abs()) < 0 - ) { - result.push(convergent); + } else { + // See https://en.wikipedia.org/wiki/Continued_fraction#Semiconvergents + // for the origin of this half-rule + try { + const halfRule = + convergent + .sub(value_) + .abs() + .compare(result[result.length - 1].sub(value_).abs()) < 0; + if (halfRule) { + result.push(convergent); + } + } catch {} } } if (result.length >= maxLength!) {