Skip to content

Commit 90b426a

Browse files
arnogarnog
authored
Merge pull request #137 from emily0622/summation_extension
Extending Summation and Product Functionality
2 parents 19d4d5f + 11e948e commit 90b426a

File tree

6 files changed

+617
-325
lines changed

6 files changed

+617
-325
lines changed

Diff for: ‎src/compute-engine/library/arithmetic-add.ts

+203-122
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,13 @@ import { Sum } from '../symbolic/sum';
77
import { asBignum, asFloat, MAX_SYMBOLIC_TERMS } from '../numerics/numeric';
88
import { widen } from '../boxed-expression/boxed-domain';
99
import { sortAdd } from '../boxed-expression/order';
10-
import { canonicalIndexingSet, normalizeIndexingSet } from './utils';
10+
import {
11+
MultiIndexingSet,
12+
SingleIndexingSet,
13+
normalizeIndexingSet,
14+
cartesianProduct,
15+
range,
16+
} from './utils';
1117
import { each, isIndexableCollection } from '../collection-utils';
1218

1319
/** The canonical form of `Add`:
@@ -114,30 +120,52 @@ export function canonicalSummation(
114120
ce: IComputeEngine,
115121
body: BoxedExpression,
116122
indexingSet: BoxedExpression | undefined
117-
) {
123+
): BoxedExpression | null {
118124
// Sum is a scoped function (to declare the index)
119125
ce.pushScope();
120126

121127
body ??= ce.error('missing');
122-
123-
indexingSet = canonicalIndexingSet(indexingSet);
124-
const result = indexingSet
125-
? ce._fn('Sum', [body.canonical, indexingSet])
126-
: ce._fn('Sum', [body.canonical]);
128+
var result: BoxedExpression | undefined = undefined;
129+
130+
if (
131+
indexingSet &&
132+
indexingSet.ops &&
133+
indexingSet.ops[0]?.head === 'Delimiter'
134+
) {
135+
var multiIndex = MultiIndexingSet(indexingSet);
136+
if (!multiIndex) return null;
137+
var bodyAndIndex = [body.canonical];
138+
multiIndex.forEach((element) => {
139+
bodyAndIndex.push(element);
140+
});
141+
result = ce._fn('Sum', bodyAndIndex);
142+
} else {
143+
var singleIndex = SingleIndexingSet(indexingSet);
144+
result = singleIndex
145+
? ce._fn('Sum', [body.canonical, singleIndex])
146+
: ce._fn('Sum', [body.canonical]);
147+
}
127148

128149
ce.popScope();
129150
return result;
130151
}
131152

132153
export function evalSummation(
133154
ce: IComputeEngine,
134-
expr: BoxedExpression,
135-
indexingSet: BoxedExpression | undefined,
155+
summationEquation: BoxedExpression[],
136156
mode: 'simplify' | 'N' | 'evaluate'
137157
): BoxedExpression | undefined {
158+
let expr = summationEquation[0];
159+
let indexingSet: BoxedExpression[] = [];
160+
if (summationEquation) {
161+
indexingSet = [];
162+
for (let i = 1; i < summationEquation.length; i++) {
163+
indexingSet.push(summationEquation[i]);
164+
}
165+
}
138166
let result: BoxedExpression | undefined | null = null;
139167

140-
if (!indexingSet) {
168+
if (indexingSet?.length === 0) {
141169
// The body is a collection, e.g. Sum({1, 2, 3})
142170
const body =
143171
mode === 'simplify'
@@ -178,141 +206,194 @@ export function evalSummation(
178206
return result ?? undefined;
179207
}
180208

181-
const [index, lower, upper, isFinite] = normalizeIndexingSet(
182-
indexingSet.evaluate()
183-
);
184-
185-
if (!index) return undefined;
209+
var indexArray: string[] = [];
210+
let lowerArray: number[] = [];
211+
let upperArray: number[] = [];
212+
let isFiniteArray: boolean[] = [];
213+
indexingSet.forEach((indexingSetElement) => {
214+
const [index, lower, upper, isFinite] = normalizeIndexingSet(
215+
indexingSetElement.evaluate()
216+
);
217+
if (!index) return undefined;
218+
indexArray.push(index);
219+
lowerArray.push(lower);
220+
upperArray.push(upper);
221+
isFiniteArray.push(isFinite);
222+
});
186223

187224
const fn = expr;
225+
const savedContext = ce.swapScope(fn.scope);
226+
ce.pushScope();
227+
fn.bind();
188228

189-
if (lower >= upper) return undefined;
229+
for (let i = 0; i < indexArray.length; i++) {
230+
const index = indexArray[i];
231+
const lower = lowerArray[i];
232+
const upper = upperArray[i];
233+
const isFinite = isFiniteArray[i];
234+
if (lower >= upper) return undefined;
190235

191-
if (mode === 'simplify' && upper - lower >= MAX_SYMBOLIC_TERMS)
192-
return undefined;
236+
if (mode === 'simplify' && upper - lower >= MAX_SYMBOLIC_TERMS)
237+
return undefined;
193238

194-
if (mode === 'evaluate' && upper - lower >= MAX_SYMBOLIC_TERMS) mode = 'N';
239+
if (mode === 'evaluate' && upper - lower >= MAX_SYMBOLIC_TERMS) mode = 'N';
195240

196-
const savedContext = ce.swapScope(fn.scope);
197-
ce.pushScope();
198-
fn.bind();
241+
if (mode === 'simplify') {
242+
const terms: BoxedExpression[] = [];
243+
for (let i = lower; i <= upper; i++) {
244+
ce.assign(index, i);
245+
terms.push(fn.simplify());
246+
}
247+
result = ce.add(terms).simplify();
248+
}
249+
}
199250

200-
if (mode === 'simplify') {
201-
const terms: BoxedExpression[] = [];
202-
for (let i = lower; i <= upper; i++) {
203-
ce.assign(index, i);
204-
terms.push(fn.simplify());
251+
// create cartesian product of ranges
252+
let cartesianArray: number[][] = [];
253+
if (indexArray.length > 1) {
254+
for (let i = 0; i < indexArray.length - 1; i++) {
255+
if (cartesianArray.length === 0) {
256+
cartesianArray = cartesianProduct(
257+
range(lowerArray[i], upperArray[i]),
258+
range(lowerArray[i + 1], upperArray[i + 1])
259+
);
260+
} else {
261+
cartesianArray = cartesianProduct(
262+
cartesianArray.map((x) => x[0]),
263+
range(lowerArray[i + 1], upperArray[i + 1])
264+
);
265+
}
205266
}
206-
result = ce.add(...terms).simplify();
267+
} else {
268+
cartesianArray = range(lowerArray[0], upperArray[0]).map((x) => [x]);
207269
}
208270

209271
if (mode === 'evaluate') {
210272
const terms: BoxedExpression[] = [];
211-
for (let i = lower; i <= upper; i++) {
212-
ce.assign(index, i);
273+
for (const element of cartesianArray) {
274+
const index = indexArray.map((x, i) => {
275+
ce.assign(x, element[i]);
276+
return x;
277+
});
213278
terms.push(fn.evaluate());
214279
}
215280
result = ce.add(...terms).evaluate();
216281
}
217282

283+
for (let i = 0; i < indexArray.length; i++) {
284+
// unassign indexes once done because if left assigned to an integer value,
285+
// in double summations the .evaluate will assume the inner index value = upper
286+
// for example in the following code latex: \\sum_{n=0}^{4}\\sum_{m=4}^{8}{n+m}`
287+
// if the indexes aren't unassigned, once the first pass is done, every following pass
288+
// will assume m is 8 for the m=4->8 iterations
289+
ce.assign(indexArray[i], undefined);
290+
}
291+
218292
if (mode === 'N') {
219-
// if (result === null && !fn.scope) {
220-
// //
221-
// // The term is not a function of the index
222-
// //
223-
224-
// const n = fn.N();
225-
// if (!isFinite) {
226-
// if (n.isZero) result = ce._ZERO;
227-
// else if (n.isPositive) result = ce._POSITIVE_INFINITY;
228-
// else result = ce._NEGATIVE_INFINITY;
229-
// }
230-
// if (result === null && fn.isPure)
231-
// result = ce.mul([ce.number(upper - lower + 1), n]);
232-
233-
// // If the term is not a function of the index, but it is not pure,
234-
// // fall through to the general case
235-
// }
236-
237-
//
238-
// Finite series. Evaluate each term and add them up
239-
//
240-
if (result === null && isFinite) {
241-
if (bignumPreferred(ce)) {
242-
let sum = ce.bignum(0);
243-
for (let i = lower; i <= upper; i++) {
244-
ce.assign(index, i);
245-
const term = asBignum(fn.N());
246-
if (term === null) {
247-
result = undefined;
248-
break;
249-
}
250-
if (!term.isFinite()) {
251-
sum = term;
252-
break;
253-
}
254-
sum = sum.add(term);
255-
}
256-
if (result === null) result = ce.number(sum);
257-
} else {
258-
// Machine precision
259-
const numericMode = ce.numericMode;
260-
const precision = ce.precision;
261-
ce.numericMode = 'machine';
262-
let sum = 0;
263-
for (let i = lower; i <= upper; i++) {
264-
ce.assign(index, i);
265-
const term = asFloat(fn.N());
266-
if (term === null) {
267-
result = undefined;
268-
break;
293+
for (let i = 0; i < indexArray.length; i++) {
294+
const index = indexArray[i];
295+
const lower = lowerArray[i];
296+
const upper = upperArray[i];
297+
const isFinite = isFiniteArray[i];
298+
// if (result === null && !fn.scope) {
299+
// //
300+
// // The term is not a function of the index
301+
// //
302+
303+
// const n = fn.N();
304+
// if (!isFinite) {
305+
// if (n.isZero) result = ce._ZERO;
306+
// else if (n.isPositive) result = ce._POSITIVE_INFINITY;
307+
// else result = ce._NEGATIVE_INFINITY;
308+
// }
309+
// if (result === null && fn.isPure)
310+
// result = ce.mul([ce.number(upper - lower + 1), n]);
311+
312+
// // If the term is not a function of the index, but it is not pure,
313+
// // fall through to the general case
314+
// }
315+
316+
//
317+
// Finite series. Evaluate each term and add them up
318+
//
319+
if (result === null && isFinite) {
320+
if (bignumPreferred(ce)) {
321+
let sum = ce.bignum(0);
322+
for (let i = lower; i <= upper; i++) {
323+
ce.assign(index, i);
324+
const term = asBignum(fn.N());
325+
if (term === null) {
326+
result = undefined;
327+
break;
328+
}
329+
if (!term.isFinite()) {
330+
sum = term;
331+
break;
332+
}
333+
sum = sum.add(term);
269334
}
270-
if (!Number.isFinite(term)) {
271-
sum = term;
272-
break;
335+
if (result === null) result = ce.number(sum);
336+
} else {
337+
// Machine precision
338+
const numericMode = ce.numericMode;
339+
const precision = ce.precision;
340+
ce.numericMode = 'machine';
341+
let sum = 0;
342+
for (let i = lower; i <= upper; i++) {
343+
ce.assign(index, i);
344+
const term = asFloat(fn.N());
345+
if (term === null) {
346+
result = undefined;
347+
break;
348+
}
349+
if (!Number.isFinite(term)) {
350+
sum = term;
351+
break;
352+
}
353+
sum += term;
273354
}
274-
sum += term;
355+
ce.numericMode = numericMode;
356+
ce.precision = precision;
357+
if (result === null) result = ce.number(sum);
275358
}
276-
ce.numericMode = numericMode;
277-
ce.precision = precision;
278-
if (result === null) result = ce.number(sum);
279-
}
280-
} else if (result === null) {
281-
//
282-
// Infinite series.
283-
//
284-
285-
// First, check for divergence
286-
ce.assign(index, 1000);
287-
const nMax = fn.N();
288-
ce.assign(index, 999);
289-
const nMaxMinusOne = fn.N();
290-
291-
const ratio = asFloat(ce.div(nMax, nMaxMinusOne).N());
292-
if (ratio !== null && Number.isFinite(ratio) && Math.abs(ratio) > 1) {
293-
result = ce.PositiveInfinity;
294-
} else {
295-
// Potentially converging series.
296-
// Evaluate as a machine number (it's an approximation to infinity, so
297-
// no point in calculating with high precision), and check for convergence
298-
let sum = 0;
299-
const numericMode = ce.numericMode;
300-
const precision = ce.precision;
301-
ce.numericMode = 'machine';
302-
for (let i = lower; i <= upper; i++) {
303-
ce.assign(index, i);
304-
const term = asFloat(fn.N());
305-
if (term === null) {
306-
result = undefined;
307-
break;
359+
} else if (result === null) {
360+
//
361+
// Infinite series.
362+
//
363+
364+
// First, check for divergence
365+
ce.assign(index, 1000);
366+
const nMax = fn.N();
367+
ce.assign(index, 999);
368+
const nMaxMinusOne = fn.N();
369+
370+
const ratio = asFloat(ce.div(nMax, nMaxMinusOne).N());
371+
if (ratio !== null && Number.isFinite(ratio) && Math.abs(ratio) > 1) {
372+
result = ce.PositiveInfinity;
373+
} else {
374+
// Potentially converging series.
375+
// Evaluate as a machine number (it's an approximation to infinity, so
376+
// no point in calculating with high precision), and check for convergence
377+
let sum = 0;
378+
const numericMode = ce.numericMode;
379+
const precision = ce.precision;
380+
ce.numericMode = 'machine';
381+
for (let i = lower; i <= upper; i++) {
382+
ce.assign(index, i);
383+
const term = asFloat(fn.N());
384+
if (term === null) {
385+
result = undefined;
386+
break;
387+
}
388+
// Converged (or diverged), early exit
389+
if (Math.abs(term) < Number.EPSILON || !Number.isFinite(term))
390+
break;
391+
sum += term;
308392
}
309-
// Converged (or diverged), early exit
310-
if (Math.abs(term) < Number.EPSILON || !Number.isFinite(term)) break;
311-
sum += term;
393+
ce.numericMode = numericMode;
394+
ce.precision = precision;
395+
if (result === null) result = ce.number(sum);
312396
}
313-
ce.numericMode = numericMode;
314-
ce.precision = precision;
315-
if (result === null) result = ce.number(sum);
316397
}
317398
}
318399
}

0 commit comments

Comments
 (0)