simplify, piecewise
Walter Roberson
> simplify(`if`(M<=N,2^(M+1) - M - 2+ (N-j)*(2^M-1) ,Q));
Error, (in tools/map) invalid input: if expects 3 arguments, but received 1
The problem starts as early as
> simplify(`if`(M<=N,2^M,Q));
M
if(M <= N, 2 , Q)
> simplify(`if`(M<=N,2*2^M,Q));
Error, (in tools/map) invalid input: if expects 3 arguments, but received 1
That is, as soon as the simplification has to recurse [due to the
exponent formula M+1 derived from 2*2^M], there are problems.
These problems can, to some extent, be postponed by using layers of
evalb() around the condition; for example, using this technique,
(N-j)*(2^M-1) can be made to pass in the centre, and even slightly
more complicated than that.
However, with the expression above, the recursion needs go beyond
any reasonable number of evalb() levels I was willing to type in.
simplify() is able to handle complex expressions if you convert to
piecewise() notation.
--
Okay, buzzwords only. Two syllables, tops. -- Laurie Anderson
Walter Roberson
> (N, M) -> piecewise(M <= N, P, hypergeom([1, -1 - N], [1 + M - N], -1)):
> codegen[optimize](%);
proc(N, M)
local t7;
option operator, arrow;
t7 := hypergeom([1, -1 - N], [1 + M - N], -1);
piecewise(M <= N, P, t7)
end proc
Notice that the possibly highly time-consuming expression
hypergeom([1, -1 - N], [1 + M - N], -1) is calculated no matter
what the result of the condition is going to turn out to be.
This is surely not an optimized procedure!!
Maple 7, Maple 10:
> codegen[optimize](M->piecewise(M=0,P,1/M+Q/M));
proc(M)
local t2;
option operator, arrow;
t2 := 1/M; piecewise(M = 0, P, t2 + Q*t2)
end proc
What happens if M is 0 ? "Bad Maple, no doggy biscuit!"
I speculate that there are similar problems in waiting for
Dirac or Heaviside.
--
Prototypes are supertypes of their clones. -- maplesoft
Walter Roberson
is unfair criticism:
Maple 7, Maple 10:
> convert(piecewise(M=0,P,1/M+Q/M),Heaviside);
Dirac(M) Dirac(M) Q
P Dirac(M) + 1/M + Q/M - -------- - ----------
M M
> eval(%,M=0);
Error, numeric exception: division by zero
> piecewise(M=0,P,1/M+Q/M);
{ P M = 0
{
{ 1/M + Q/M otherwise
> eval(%,M=0);
P
We see that in the conversion to Heaviside, no precautions were
taken to ensure that the terms would all be calculable outside of the
domain of the restrictions implied by the piecewise() conditions.
Similarily,
> convert(piecewise(M=Pi/2,P,tan(M)),Heaviside);
Pi Pi
P Dirac(- ---- + M) + tan(M) - tan(M) Dirac(- ---- + M)
2 2
> subs(M=Pi/2,%);
Pi Pi
P Dirac(0) + tan(----) - tan(----) Dirac(0)
2 2
> simplify(%);
Error, (in tan) numeric exception: division by zero
This would obviously be harder to catch than problems with elementary
operations such as division by 0 -- but isn't a substantial
point of piecewise() that you can use it to stay within boundary
conditions ? The conversion to Heaviside somehow seems incomplete
to me if the boundary conditions are not respected.
--
There are some ideas so wrong that only a very intelligent person
could believe in them. -- George Orwell
Axel Vogt
:-((
At least it works with the option 'tryhard', it seems
to respect that (ie: does not change it):
tmp:= M -> piecewise(M=0,P,1/M+Q/M);
tmp(M);
tmp:=makeproc(%);
tmp:=optimize(tmp,tryhard);
showstat(tmp);
tmp:=(N, M) -> piecewise(M <= N, P, hypergeom([1, -1 - N], [1 + M - N], -1)):
tmp(N, M);
codegen[makeproc](%);
tmp:=optimize(tmp,tryhard);
showstat(tmp);
For the 'convert' (in your subsequent post) I would like
to know, whether M states that the expressions are really
equivalent (or something like using simplify with option
symbolic acn happen).
Walter Roberson
>At least it works with the option 'tryhard', it seems
>to respect that (ie: does not change it):
>tmp:= M -> piecewise(M=0,P,1/M+Q/M);
>tmp(M);
>tmp:=makeproc(%);
>tmp:=optimize(tmp,tryhard);
>showstat(tmp);
With the tryhard option, 10.04 codegen[optimize] just leaves the
piecewise expression unoptimized. For example, it doesn't even
remove the common subexpression for
piecewise(M=0,M^2+5,1/(M^2+5)+Q/(M^2+5))
--
I was very young in those days, but I was also rather dim.
-- Christopher Priest
Walter Roberson
>to know, whether M states that the expressions are really
>equivalent (or something like using simplify with option
>symbolic acn happen).
That refers to a convert() of a piecewise to Heaviside in which
the boundary conditions were not checked.
According to ?convert,Heaviside the output might differ from
the original at finitely many points -- so failing to check
for M=0 in convert(piecewise(M=0,P,1/M),Heaviside) would fall
within the 'finitely many points' limitation.
Consider, though:
> convert(piecewise(sin(Pi*M)=0,P,1/sin(Pi*M)),Heaviside);
1 Dirac(sin(Pi M))
P Dirac(sin(Pi M)) + --------- - ----------------
sin(Pi M) sin(Pi M)
This differs from the piecewise result at infinitely many points
(i.e., all integers).
--
"It is important to remember that when it comes to law, computers
never make copies, only human beings make copies. Computers are given
commands, not permission. Only people can be given permission."
-- Brad Templeton
Walter Roberson
Walter Roberson <robe...@ibd.nrc-cnrc.gc.ca> wrote:
>Maple 10:
Here's a case where codegen[optimize] 'tryhard' produces provably
wrong code:
> F1 := unapply(sum(k*f(x-k),k=M..N)/sum(f(x-k),k=M..N),x,M,N):
> F2 := codegen[optimize](F1,'tryhard'):
> F1(x,M,N);
N
-----
\
) k f(x - k)
/
-----
k = M
----------------
N
-----
\
) f(x - k)
/
-----
k = M
> F2(x,M,N);
2 2
(1/2 f(x - k) (N + 1) - 1/2 (N + 1) f(x - k) - 1/2 f(x - k) M
+ 1/2 M f(x - k))/((N - M + 1) f(x - k))
> F1(x,3,5);
3 f(x - 3) + 4 f(x - 4) + 5 f(x - 5)
------------------------------------
f(x - 3) + f(x - 4) + f(x - 5)
> F2(x,3,5);
4
> eval(F1(x,3,5),{x=7,f=proc(t) t^2 end proc});
104
---
29
> op(F2);
proc(x, M, N)
local result, t96;
global k;
t96 := f(x - k);
result := (
1/2*t96*(1 + N)^2 - 1/2*(1 + N)*t96 - 1/2*t96*M^2 + 1/2*M*t96)/(
(N - M + 1)*t96)
end proc
Notice that the summations have disappeared, and that k has become
a global variable with a free value -- whereas in the original it
is the dummy variable of a definite summation. Even if one
uses a modified version of F1 in which k is explicitly declared
as local, the only change to the optimized version is that k will
be local in the revised F2 -- but will never be given a value.
--
"law -- it's a commodity"
-- Andrew Ryan (The Globe and Mail, 2005/11/26)
Walter Roberson
Walter Roberson <robe...@ibd.nrc-cnrc.gc.ca> wrote:
(some problems with codegen[optimize])
Maple 10:
> simplify(`if`(a>b,a,b));
if(b < a, a, b)
> codegen[optimize](%);
t2 = if(b < a, a, b)
> codegen[optimize](%,tryhard);
[Ummm, it's still running 330 seconds later, on a fast machine!]
--
If you lie to the compiler, it will get its revenge. -- Henry Spencer
Walter Roberson
Walter Roberson <robe...@ibd.nrc-cnrc.gc.ca> wrote:
>In article <eeerp8$d1b$1...@canopus.cc.umanitoba.ca>,
>Walter Roberson <robe...@ibd.nrc-cnrc.gc.ca> wrote:
>Maple 10:
>> simplify(`if`(a>b,a,b));
> if(b < a, a, b)
>
>> codegen[optimize](%);
> t2 = if(b < a, a, b)
>
>> codegen[optimize](%,tryhard);
>[Ummm, it's still running 330 seconds later, on a fast machine!]
I interrupted and quit maple and re-entered and tried the above
sequence, and it gave an answer in a fraction of a second. That
suggests that something else I'd done in the session triggered the
loop behaviour. Unfortunately those previous things have scrolled
off, and I can't seem to reproduce this easily.
--
"No one has the right to destroy another person's belief by
demanding empirical evidence." -- Ann Landers
Axel Vogt
I think that error comes through sum, not through simplification,
as it works with Sum, but I have not located the error's source.
The evaluated form has a behaviour I really dislike. From the help:
"The call sum(f, k=m..n) computes the definite sum of f(k) over
the given range m..n. That is, it computes f(m) + f(m+1) + ... +
f(n). If m = n+1, the value returned is 0. If m > n+1, the value
returned is -sum(f, k=n+1..m-1)."
For me a sum is over index in a set, ordered or not, and we do
it in a commutative group - where the 'extra rule' can not hold.
Walter Roberson
>Walter Roberson wrote:
>> > F1 := unapply(sum(k*f(x-k),k=M..N)/sum(f(x-k),k=M..N),x,M,N):
>> > F2 := codegen[optimize](F1,'tryhard'):
>> > F2(x,M,N);
>> 2 2
>> (1/2 f(x - k) (N + 1) - 1/2 (N + 1) f(x - k) - 1/2 f(x - k) M
>>
>> + 1/2 M f(x - k))/((N - M + 1) f(x - k))
>I think that error comes through sum, not through simplification,
>as it works with Sum, but I have not located the error's source.
`codegen/optimize/tryhard/optimize` near line 76,
is extracting the f(x-k) as a constant, and moving it outside
the summation, which of course produces bogus answers.
[t3 = f(x-k), n[1] = sum(k*t3,k = M .. N)/sum(t3,k = M .. N)]
This will probably occur for any constant-looking expression
that actually depends upon an implicit sequencing, such as in
sum(), prod(), seq(), or `$` .
Walter Roberson
Walter Roberson <robe...@ibd.nrc-cnrc.gc.ca> wrote:
>Maple 10:
>The problem starts as early as
>> simplify(`if`(M<=N,2*2^M,Q));
>Error, (in tools/map) invalid input: if expects 3 arguments, but received 1
>That is, as soon as the simplification has to recurse [due to the
>exponent formula M+1 derived from 2*2^M], there are problems.
I have traced down this problem, along with the related one that
shows up in PDEtools/eval/2 when simplifying other expressions
inside of `if`. The two bugs have the same root cause.
In both cases, the code is doing roughly
op(0,zz)(op(zz))
where zz is the unevaluated `if`(condition,truecase,falsecase)
Because `if` has special evaluation semantics, the op(zz)
is not evaluated, and `if`() is passed the unevaluated 'op(zz)'.
It then complains because there is only one argument instead of 3.
Code akin to this will work:
subs(ARGS=op(zz),FUNC=op(0,zz),FUNC(ARGS))
though possibly a layer of eval() around that is called for.
The above fix might perhaps break on other functions with special
evaluation semantics, such as for a user defined function that declares
some arguments as ::uneval or ::evaln
If anyone is interested in patching their own copies, then in
10.04 the lines are `tools/map` line 24 and `PDEtools/eval/2` line 23 .