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,