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

Jump to bottom

Extending Summation and Product Functionality #137

Merged
merged 3 commits into from
Dec 18, 2023
Merged

Extending Summation and Product Functionality #137

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a5e3023
extended summation functionality to parse a sequence in the sub and s…
Dec 16, 2023
80f878a
extended evalSummation and evalProduct
Dec 17, 2023
11e948e
Merge branch 'main' into summation_extension
arnog Dec 18, 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
325 changes: 203 additions & 122 deletions src/compute-engine/library/arithmetic-add.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,13 @@ import { Sum } from '../symbolic/sum';
import { asBignum, asFloat, MAX_SYMBOLIC_TERMS } from '../numerics/numeric';
import { widen } from '../boxed-expression/boxed-domain';
import { sortAdd } from '../boxed-expression/order';
import { canonicalIndexingSet, normalizeIndexingSet } from './utils';
import {
MultiIndexingSet,
SingleIndexingSet,
normalizeIndexingSet,
cartesianProduct,
range,
} from './utils';
import { each, isIndexableCollection } from '../collection-utils';

/** The canonical form of `Add`:
Expand Down Expand Up @@ -114,30 +120,52 @@ export function canonicalSummation(
ce: IComputeEngine,
body: BoxedExpression,
indexingSet: BoxedExpression | undefined
) {
): BoxedExpression | null {
// Sum is a scoped function (to declare the index)
ce.pushScope();

body ??= ce.error('missing');

indexingSet = canonicalIndexingSet(indexingSet);
const result = indexingSet
? ce._fn('Sum', [body.canonical, indexingSet])
: ce._fn('Sum', [body.canonical]);
var result: BoxedExpression | undefined = undefined;

if (
indexingSet &&
indexingSet.ops &&
indexingSet.ops[0]?.head === 'Delimiter'
) {
var multiIndex = MultiIndexingSet(indexingSet);
if (!multiIndex) return null;
var bodyAndIndex = [body.canonical];
multiIndex.forEach((element) => {
bodyAndIndex.push(element);
});
result = ce._fn('Sum', bodyAndIndex);
} else {
var singleIndex = SingleIndexingSet(indexingSet);
result = singleIndex
? ce._fn('Sum', [body.canonical, singleIndex])
: ce._fn('Sum', [body.canonical]);
}

ce.popScope();
return result;
}

export function evalSummation(
ce: IComputeEngine,
expr: BoxedExpression,
indexingSet: BoxedExpression | undefined,
summationEquation: BoxedExpression[],
mode: 'simplify' | 'N' | 'evaluate'
): BoxedExpression | undefined {
let expr = summationEquation[0];
let indexingSet: BoxedExpression[] = [];
if (summationEquation) {
indexingSet = [];
for (let i = 1; i < summationEquation.length; i++) {
indexingSet.push(summationEquation[i]);
}
}
let result: BoxedExpression | undefined | null = null;

if (!indexingSet) {
if (indexingSet?.length === 0) {
// The body is a collection, e.g. Sum({1, 2, 3})
const body =
mode === 'simplify'
Expand Down Expand Up @@ -178,141 +206,194 @@ export function evalSummation(
return result ?? undefined;
}

const [index, lower, upper, isFinite] = normalizeIndexingSet(
indexingSet.evaluate()
);

if (!index) return undefined;
var indexArray: string[] = [];
let lowerArray: number[] = [];
let upperArray: number[] = [];
let isFiniteArray: boolean[] = [];
indexingSet.forEach((indexingSetElement) => {
const [index, lower, upper, isFinite] = normalizeIndexingSet(
indexingSetElement.evaluate()
);
if (!index) return undefined;
indexArray.push(index);
lowerArray.push(lower);
upperArray.push(upper);
isFiniteArray.push(isFinite);
});

const fn = expr;
const savedContext = ce.swapScope(fn.scope);
ce.pushScope();
fn.bind();

if (lower >= upper) return undefined;
for (let i = 0; i < indexArray.length; i++) {
const index = indexArray[i];
const lower = lowerArray[i];
const upper = upperArray[i];
const isFinite = isFiniteArray[i];
if (lower >= upper) return undefined;

if (mode === 'simplify' && upper - lower >= MAX_SYMBOLIC_TERMS)
return undefined;
if (mode === 'simplify' && upper - lower >= MAX_SYMBOLIC_TERMS)
return undefined;

if (mode === 'evaluate' && upper - lower >= MAX_SYMBOLIC_TERMS) mode = 'N';
if (mode === 'evaluate' && upper - lower >= MAX_SYMBOLIC_TERMS) mode = 'N';

const savedContext = ce.swapScope(fn.scope);
ce.pushScope();
fn.bind();
if (mode === 'simplify') {
const terms: BoxedExpression[] = [];
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
terms.push(fn.simplify());
}
result = ce.add(terms).simplify();
}
}

if (mode === 'simplify') {
const terms: BoxedExpression[] = [];
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
terms.push(fn.simplify());
// create cartesian product of ranges
let cartesianArray: number[][] = [];
if (indexArray.length > 1) {
for (let i = 0; i < indexArray.length - 1; i++) {
if (cartesianArray.length === 0) {
cartesianArray = cartesianProduct(
range(lowerArray[i], upperArray[i]),
range(lowerArray[i + 1], upperArray[i + 1])
);
} else {
cartesianArray = cartesianProduct(
cartesianArray.map((x) => x[0]),
range(lowerArray[i + 1], upperArray[i + 1])
);
}
}
result = ce.add(...terms).simplify();
} else {
cartesianArray = range(lowerArray[0], upperArray[0]).map((x) => [x]);
}

if (mode === 'evaluate') {
const terms: BoxedExpression[] = [];
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
for (const element of cartesianArray) {
const index = indexArray.map((x, i) => {
ce.assign(x, element[i]);
return x;
});
terms.push(fn.evaluate());
}
result = ce.add(...terms).evaluate();
}

for (let i = 0; i < indexArray.length; i++) {
// unassign indexes once done because if left assigned to an integer value,
// in double summations the .evaluate will assume the inner index value = upper
// for example in the following code latex: \\sum_{n=0}^{4}\\sum_{m=4}^{8}{n+m}`
// if the indexes aren't unassigned, once the first pass is done, every following pass
// will assume m is 8 for the m=4->8 iterations
ce.assign(indexArray[i], undefined);
}

if (mode === 'N') {
// if (result === null && !fn.scope) {
// //
// // The term is not a function of the index
// //

// const n = fn.N();
// if (!isFinite) {
// if (n.isZero) result = ce._ZERO;
// else if (n.isPositive) result = ce._POSITIVE_INFINITY;
// else result = ce._NEGATIVE_INFINITY;
// }
// if (result === null && fn.isPure)
// result = ce.mul([ce.number(upper - lower + 1), n]);

// // If the term is not a function of the index, but it is not pure,
// // fall through to the general case
// }

//
// Finite series. Evaluate each term and add them up
//
if (result === null && isFinite) {
if (bignumPreferred(ce)) {
let sum = ce.bignum(0);
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
const term = asBignum(fn.N());
if (term === null) {
result = undefined;
break;
}
if (!term.isFinite()) {
sum = term;
break;
}
sum = sum.add(term);
}
if (result === null) result = ce.number(sum);
} else {
// Machine precision
const numericMode = ce.numericMode;
const precision = ce.precision;
ce.numericMode = 'machine';
let sum = 0;
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
const term = asFloat(fn.N());
if (term === null) {
result = undefined;
break;
for (let i = 0; i < indexArray.length; i++) {
const index = indexArray[i];
const lower = lowerArray[i];
const upper = upperArray[i];
const isFinite = isFiniteArray[i];
// if (result === null && !fn.scope) {
// //
// // The term is not a function of the index
// //

// const n = fn.N();
// if (!isFinite) {
// if (n.isZero) result = ce._ZERO;
// else if (n.isPositive) result = ce._POSITIVE_INFINITY;
// else result = ce._NEGATIVE_INFINITY;
// }
// if (result === null && fn.isPure)
// result = ce.mul([ce.number(upper - lower + 1), n]);

// // If the term is not a function of the index, but it is not pure,
// // fall through to the general case
// }

//
// Finite series. Evaluate each term and add them up
//
if (result === null && isFinite) {
if (bignumPreferred(ce)) {
let sum = ce.bignum(0);
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
const term = asBignum(fn.N());
if (term === null) {
result = undefined;
break;
}
if (!term.isFinite()) {
sum = term;
break;
}
sum = sum.add(term);
}
if (!Number.isFinite(term)) {
sum = term;
break;
if (result === null) result = ce.number(sum);
} else {
// Machine precision
const numericMode = ce.numericMode;
const precision = ce.precision;
ce.numericMode = 'machine';
let sum = 0;
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
const term = asFloat(fn.N());
if (term === null) {
result = undefined;
break;
}
if (!Number.isFinite(term)) {
sum = term;
break;
}
sum += term;
}
sum += term;
ce.numericMode = numericMode;
ce.precision = precision;
if (result === null) result = ce.number(sum);
}
ce.numericMode = numericMode;
ce.precision = precision;
if (result === null) result = ce.number(sum);
}
} else if (result === null) {
//
// Infinite series.
//

// First, check for divergence
ce.assign(index, 1000);
const nMax = fn.N();
ce.assign(index, 999);
const nMaxMinusOne = fn.N();

const ratio = asFloat(ce.div(nMax, nMaxMinusOne).N());
if (ratio !== null && Number.isFinite(ratio) && Math.abs(ratio) > 1) {
result = ce.PositiveInfinity;
} else {
// Potentially converging series.
// Evaluate as a machine number (it's an approximation to infinity, so
// no point in calculating with high precision), and check for convergence
let sum = 0;
const numericMode = ce.numericMode;
const precision = ce.precision;
ce.numericMode = 'machine';
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
const term = asFloat(fn.N());
if (term === null) {
result = undefined;
break;
} else if (result === null) {
//
// Infinite series.
//

// First, check for divergence
ce.assign(index, 1000);
const nMax = fn.N();
ce.assign(index, 999);
const nMaxMinusOne = fn.N();

const ratio = asFloat(ce.div(nMax, nMaxMinusOne).N());
if (ratio !== null && Number.isFinite(ratio) && Math.abs(ratio) > 1) {
result = ce.PositiveInfinity;
} else {
// Potentially converging series.
// Evaluate as a machine number (it's an approximation to infinity, so
// no point in calculating with high precision), and check for convergence
let sum = 0;
const numericMode = ce.numericMode;
const precision = ce.precision;
ce.numericMode = 'machine';
for (let i = lower; i <= upper; i++) {
ce.assign(index, i);
const term = asFloat(fn.N());
if (term === null) {
result = undefined;
break;
}
// Converged (or diverged), early exit
if (Math.abs(term) < Number.EPSILON || !Number.isFinite(term))
break;
sum += term;
}
// Converged (or diverged), early exit
if (Math.abs(term) < Number.EPSILON || !Number.isFinite(term)) break;
sum += term;
ce.numericMode = numericMode;
ce.precision = precision;
if (result === null) result = ce.number(sum);
}
ce.numericMode = numericMode;
ce.precision = precision;
if (result === null) result = ce.number(sum);
}
}
}
Expand Down
Loading