Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

More sane approximation #12

Merged
merged 2 commits into from
Dec 27, 2023
Merged

More sane approximation #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
22b15a3
Make prime limit approximator more sane
frostburn Dec 26, 2023
e6e9e2c
Make convergent calculator more sane
frostburn Dec 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions src/__tests__/approximation.spec.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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', () => {
Expand Down Expand Up @@ -114,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', () => {
Expand Down
85 changes: 51 additions & 34 deletions src/approximation.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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!) {
Expand Down Expand Up @@ -168,16 +172,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]);
}
Expand All @@ -187,7 +194,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;
}
}
Expand All @@ -200,6 +207,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);
Expand All @@ -208,34 +218,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);
Expand All @@ -259,7 +276,7 @@ export function approximatePrimeLimit(
limitIndex: number,
maxExponent: number,
maxError = 600,
maxLength?: number
maxLength = 100
) {
return approximatePrimeLimitWithErrors(
cents,
Expand Down