#734 term-substitution in function sum is too late

closed
nobody
5
2005-11-24
2005-04-30
Anonymous
No

(%i1) "Lower Riemann-Sum, Intervall [0,1], 10 parts"\$
(%i2) /*Function*/ f(x):=x^2\$
(%i3) /*area width*/ b: 1/10\$
(%i4) /*area depending on k=0,...,9*/ Ak: b*f(k*b)\$
(%i5) Ak;
2
k
(%o5) ----
1000
(%i6) /* wrong */ sum(Ak,k,0,9);
2
k
(%o6) ---
100
obviously, we first have the summation and then the
symbol-substitution
(%i7) sum(''Ak,k,0,9);
57
(%o7) ---
200

finally, this is ok, but the usage of the double-quote
is hard to explain to beginners.

it gets even worse, if you use a function
A(k):=k*f(k*b)\$ for a sum with n partitions:
(%i8) "Lower Riemann-Sum, Intervall [0,1], n parts"\$
(%i9) /*area width*/ b:1/n\$
(%i10) /*area depending on k=0,...,n-1*/ A(k):=k*f(k*b)\$
(%i11) A(k);
3
k
(%o11) --
2
n
(%i12) sum(A(k),k,0,n-1);
n - 1
====
\ (%o12) > A(k)
/
====
k = 0

in that case you need ''(A(k)), a bracket and the
double-quote.

is it possible to patch the sum-function in that way,
that we first have the term- or symbol-substitution and
second the summation?

Discussion

• Barton Willis - 2005-05-01

Logged In: YES
user_id=895922

We've discussed this before; see

Maybe a simple workaround such as

mysum(f,v,lo,hi):=block([acc:0],
if integerp(lo) and integerp(hi) then for i from lo thru
hi do acc:acc+substitute(i,v,f)
else acc:funmake('mysum,[f,v,lo,hi]),
acc)

or

mysum(f,lo,hi):= block([acc:0],
if integerp(lo) and integerp(hi) then
for i : lo thru hi do acc : acc + apply(f,[i])
else acc : funmake('mysum,[f,lo,hi]),
acc);

might work for you,

Barton

• James Amundson - 2005-05-30
• labels: 460521 --> Lisp Core
• milestone: 482214 -->

• Robert Dodier - 2005-06-13

Logged In: YES
user_id=501686

OK, for the record, here is a definition of \$sum which
implements one of the ideas that has been proposed, namely
this policy: evaluate the summand after binding the
summation variable to itself.

(defmspec \$sum (l) (setq l (cdr l))
(if (= (length l) 4)
(wna-err '\$sum)))

This version yields the results expected by the person who
originated this bug report, and yields the results predicted
by S Macrakis in the email referenced below
(http://www.math.utexas.edu/pipermail/maxima/2003/004869.html),
and run_testsuite() runs to completion with no errors.

I'm in favor of making this change, but I'm just recording
this code snippet here for future reference; no changes
planned at the moment.

For comparison here is the current version of \$sum:

(defmspec \$sum (l) (setq l (cdr l))
(if (= (length l) 4)
(wna-err '\$sum)))

• Barton Willis - 2005-06-13

Logged In: YES
user_id=895922

If we decide to use Robert's \$sum function, the translate
property for sum will need to be changed. Consider:

(1) Redefine \$sum as Robert suggested.

(2) Try this:

(%i1) ak : k^2\$
(%i2) g(a,n) := sum(a,k,1,n)\$
(%i3) g(ak,5);
(%o3) 55 <--- correct for Robert's \$sum function

(%i4) translate(g);
(%o4) [g]
(%i5) g(ak,5);
(%o5) 5*k^2 <---- yeech

There might be other things that need fixing:

(%i8) properties(sum);
(%o8) [Special Evaluation Form,OUTATIVE,NOUN,RULE]

Barton

• Robert Dodier - 2005-06-14

Logged In: YES
user_id=501686

Actually, after looking at the code some more (specifically
DOSUM in src/asum.lisp) it looks to me like the intent is
indeed to evaluate the summand after binding the summation
index to itself. In addition, there is code to assume the
index is between its lower and upper limits (although
assume/forget is handled incorrectly, see 851765). I believe
that the observed behavior is a bug, in the narrow sense
that the observed behavior is different from what the code

Here is a different patch. I like this one better than the
previous one.

In DOSUM:

;; ORIGINAL: MEVAL ONCE
;; (mbinding (lind l*i) (meval exp))
(mbinding (lind l*i) (meval (meval exp)))

In MEVALSUMARG:

;; ORIGINAL: MEVALATOMS
;; (resimplify (mevalatoms (if (and (not (atom exp))
(resimplify (meval (if (and (not (atom exp))

With these changes, the results are as expected for the
examples cited by the original poster, and also the example
about the translated function produces the correct result,
and some examples adapted from another bug report (740134)
yield correct results.

ex 1:

f(x):=x^2\$ b: 1/10\$ Ak: b*f(k*b)\$

sum(Ak,k,0,9); => 57/200

ex 2:

A(k):=k*f(k*b)\$ b:1/n\$

A(k) => k^3/n^2

sum(A(k),k,0,n-1) => ('sum(k^3,k,0,n-1))/n^2

sum(A(k),k,0,n-1), simpsum => ((n-1)^4 + 2*(n-1)^3 +
(n-1)^2) / (4*n^2)

ex 3:

ak : k^2\$ g(a,n) := sum(a,k,1,n)\$

g(ak,5) => 55

translate (g)\$

g(ak,5) => 55

some other examples, adapted from bug # 740134:

ex 4:

sum (print (i), i, 1, 3) => prints 1, 2, 3 then returns 6

ex 5:

sum (integrate (x^i ,x), i, 0, 2) => x^3/3 + x^2/2 + x

ex 6:

sum (integrate (1/(x^i + 1), x), i, 0, 1) => log(x+1)
+ x/2

ex 7:

f[i](x):=x^i\$ g[i](x):=x^i\$ h[i](x):=x^i\$

/* reference f[i] and g[i] -- see 740134 for the effect
this has on previous defn of sum */

f[i]\$ g[i](t)\$

sum (f[i](x), i, 0, n) => 'sum (x^i, i, 0, n)
sum (g[i](x), i, 0, n) => 'sum (x^i, i, 0, n)
sum (h[i](x), i, 0, n) => 'sum (x^i, i, 0, n)

ex 8:

sum (integrate (x^i, x), i, 0, n) => 'sum (x^(i+1) /
(i+1), i, 0, n)

• Robert Dodier - 2005-06-14

Logged In: YES
user_id=501686

OK, here is another attempt. I think I now understand this
comment by Stavros in bug report # 740134 -- "First evaluate
foo with i bound to itself, and check if that result is free
of i. If so, return the product. If not, *substitute* (don't
evaluate) i=lowerlimit, i=lowerlimit+1, etc."

;; ORIGINAL: MEVAL ONCE
;; (mbinding (lind l*i) (meval exp))
;; TRIED THIS: MEVAL TWICE
;; (mbinding (lind l*i) (meval (meval exp)))
;; THIRD TIME IS A CHARM ??
(\$substitute *i ind (mbinding (lind l*i)
(meval exp)))

This version now yields this result --

ex 9:

f(x) := sum (x, i, 1, 3);

f(-x) => -3 x

as expected, as well as the same results for ex 1 through ex 8.

• Robert Dodier - 2005-06-15

Logged In: YES
user_id=501686

I'm attaching a list of test cases for sum, which can be
executed by batch ("rtestsum.mac", test) . Doubtless there
will be some discussion as to what the expected output of
each test should be.

All tests pass after making these 2 changes:

;; ORIGINAL: MEVAL ONCE
;; (mbinding (lind l*i) (meval exp))
;; TRIED THIS: MEVAL TWICE
;; (mbinding (lind l*i) (meval (meval exp)))
;; THIRD TIME IS A CHARM ??
(\$substitute *i ind (mbinding (lind l*i) (meval exp)))

;; ORIGINAL: MEVALATOMS
;; (resimplify (mevalatoms (if (and (not (atom exp))
(resimplify (meval (if (and (not (atom exp))

• Volker van Nek - 2005-06-18

Logged In: YES
user_id=1269745

Dear Maxima-developers,
I tried the last fix from Robert Dodier and it works.
Volker van Nek
original (anonymous) poster

• Robert Dodier - 2005-06-19

Logged In: YES
user_id=501686

I'm removing rtestsum.mac, in anticipation of attaching a
new version in a minute.

• Robert Dodier - 2005-06-19

Logged In: YES
user_id=501686

Attaching a longer list of test cases. batch
("rtestsum.mac", test); executes the test. As always the
"correct" results (2nd line in each pair) are subject to review.

• Robert Dodier - 2005-06-19

Logged In: YES
user_id=501686

I'm attaching my latest and greatest revision
(tmp-sum5.lisp) of DOSUM, DO%SUM, and MEVALSUMARGS.

I find that

batch ("rtestsum.mac", test);

where rtestsum.mac is the currently attached version, yields
all tests passed for Maxima 5.9.1cvs built with clisp 2.33.2
and with gcl 2.6.6.

This revision uses a gensym index for cases handled by
MEVALSUMARGS (i.e., cases for which (upper limit) minus
(lower limit) is not an integer); the gensym index is the
argument for assume/forget. When upper minus lower is an
integer, the index is bound to lower, lower+1, lower+2,
etc., the summand is evaluated, and then lower, lower+1,
lower+2, etc. is substituted for any remaining instances of
index.

• Robert Dodier - 2005-06-20

Logged In: YES
user_id=501686

Removing rtestsum.mac, will upload a new version momentarily.

• Robert Dodier - 2005-06-20

Logged In: YES
user_id=501686

New version of rtestsum.mac. Two new test cases, sum (lambda
([i], i^2, 1, 3) and sum (lambda ([i], i^2), i, 1, n).

• Robert Dodier - 2005-06-20

Logged In: YES
user_id=501686

Yet another revision (attached file tmp-sum6.lisp) of DOSUM
and friends. This version passes run_testsuite() and batch
("rtestsum.mac", test) using the currently attached version
of rtestsum.mac.

This version avoids substituting for the index in
expressions which are free of the index.

This version also has an alternate definition of the
function FREE. The simplification code (and other code)
calls FREE instead of FREEOF. FREEOF looks for dummy
variables, but FREE does not, so with the original defn of
FREE, sum (lambda([i], i^2), i, 1, n) does not simplify to n
lambda ([i], i^2), i, 1, n). I'm not very happy about
changing FREE since it is called from all over the place,
but I'm also not very happy about the presence of FREE,
either; why not call FREEOF ?

I don't think changing FREE at this point is such a good
idea; it seems to me the two choices are (1) inspect the
sum/product simplification code and change FREE to FREEOF
where possible. That would be a lot of work. (2) Leave FREE
as it is, and just accept that expressions with dummy
variables can't be factored out of summations.

At this point I think I'd rather leave FREE alone, and file
a separate bug report about the inability to factor out
expressions with dummy variables, and just leave it for
another day.

• Barton Willis - 2005-06-20

Logged In: YES
user_id=895922

Does tmp-sum6.lisp need a ratdisrep somewhere?

(%i2) sum(k,k,1,3);
(%o2) 6
(%i3) sum(rat(k),k,1,3);
Maxima encountered a Lisp error:
Error in COND [or a callee]: Caught fatal error [memory may
be damaged]
Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.
(%i4)

Barton

• Robert Dodier - 2005-06-20

Logged In: YES
user_id=501686

makes it happy again. Oh well. I think I'll just hang onto
the current rev for a day or two to see what other changes
accumulate, before posting a new version of tmp-sum*.lisp.

A related bit -- sum(rat(k), k, 1, n) returns 'sum (k, k, 1,
n), which isn't exactly wrong, but it does make me wonder
what happened to the rat... it turns out this is what happened:

subst (foo, k, rat(k)) => foo

Not rat(foo) as I would have expected. Dunno what this is about.

• Robert Dodier - 2005-06-21

Logged In: YES
user_id=501686

Delete rtestsum.mac, going to upload a new one momentarily.

• Robert Dodier - 2005-06-21

Logged In: YES
user_id=501686

New rtestsum.mac. A couple of new tests for sum, and also
copied all sum tests and substituted product for sum (and

• Robert Dodier - 2005-06-21

Logged In: YES
user_id=501686

New DOSUM, etc, in tmp-sum8.lisp. A few new sum tests have
been added to rtestsum.mac, and all sum tests copied and
turned into product tests. This revision passes all tests in
current rtestsum.mac except for test 38:

(assume(i < 0, j < 0, n > 0), product(product(abs(j) +
abs(i), i, 1, j), j, 1, n));

=> Lower bound to product is > upper bound.

This emanates from SIMPROD1 (src/combin.lisp) which is
called after the assumptions about i and j (actually,
assumptions about the gensym stand-ins for i and j) have
been forgotten.

This looks like a bug in combin.lisp -- there is no

• Robert Dodier - 2005-06-27

Logged In: YES
user_id=501686

attempting to remove old version of rtestsum.mac and upload
new version in the same update... let's see if it works.

• Robert Dodier - 2005-08-20

zip file containing sum reimplementation

• Robert Dodier - 2005-08-20

Logged In: YES
user_id=501686

attaching a zip file containing unksum-6.lisp,
tmp-sum13.lisp (two different sum reimplementations),
rtestsum.mac, and mand.lisp (a reimplementation of
mcond-related stuff -- needed for the one test which makes
use of an unevaluated conditional).

• Robert Dodier - 2005-08-20

Logged In: YES
user_id=501686

removing old versions of files

• Robert Dodier - 2005-11-24
• status: open --> closed

• Robert Dodier - 2005-11-24

Logged In: YES
user_id=501686

The reported bug is not present in the current cvs version of
Maxima.

Thank you for your report. If you see this bug in a later version
of Maxima, please submit a new bug report.