Errors-To: maple-test-request@old.daisy.uwaterloo.ca
Message-Id:   <9407301906.AA16315@old.daisy.uwaterloo.ca>
Newsgroups:   bit.listserv.maple-l
Date:         Sat, 30 Jul 1994 15:06:09 -0400
Reply-To: "MAPLE-L Discussion on MAPLE software. Local distribution only."              <MAPLE-L@IRLEARN.BITNET>
Sender: "MAPLE-L Discussion on MAPLE software. Local distribution only."              <MAPLE-L@IRLEARN.BITNET>
From: Maple User Group <maple_gr@OLD.DAISY.UWATERLOO.CA>
Subject:      Re: Error in Assume
To: Lmail Lmail <Lmail.Lmail@CHEMIE.UNI-REGENSBURG.DE>
Status: R

     From: Laurent.Bernardin@inf.ethz.ch
     Subject: Re: Error in Assume

-> From: Maple User Group <maple_gr@old.daisy.uwaterloo.ca>
-> Subject: Re: Error in Assume
->
->
->      From: Joe Riel <jsr@sparc.SanDiegoCA.NCR.COM>
->      Subject: Re: Error in Assume
->
-> The error is not with the assume facility, but rather with simplify, and
-> more specifically with `simplify/sqrt/fraction`.  The following shows
-> your example without using assume.
->
->   tmp := sqrt(a/b)*sqrt(b)/sqrt(a):
->
->   simplify(tmp, symbolic);
->
->                   sqrt(a)
->
->   `simplify/sqrt/fraction`(tmp, symbolic);
->
->                   sqrt(a)
->
-> Joe Riel
->
-> >
-> >
-> >      From: "Charles Warren" <WARREN@chem1.chem.dal.ca>
-> >      Subject:       Error in Assume
-> >
-> > While doing some integations of products of harmonic oscillator
-> > wavefunctions I obtained incorrect results. I traced the
-> > errors to incorrect behaviour of the assume facility. This behaviour
-> > is illustrated below.
-> >
-> > Maple 5 Release 3 gives
-> > an error when assume is used for the following expression
-> >
-> > assume(alpha,positive);
-> > sqrt(alpha/Pi)*sqrt(Pi)/sqrt(alpha);
-> > simplify(");
-> >
-> >
-> > The answer given is sqrt(alpha*Pi)/sqrt(Pi) instead of 1.
-> >
-> > If the above expression is rearranged as
-> >
-> > assume(alpha,positive);
-> > sqrt(alpha/Pi)*sqrt(Pi/alpha);
-> > simplify(");
-> >
-> > then the correct answer 1 is given.
-> >
->

This problem has been fixed. Below you'll find an updated copy of `simplify/sqrt/fraction`.
One way to use it is to save it into a file (e.g. sqrtfraction) and then start maple and
do:
> readlib(`simplify/sqrt`):
> read sqrtfraction;


Sorry for any inconvenience this bug may have caused.

     |\^/|      ---------------------------------------------------------
 ._|\|   |/|_. |       Laurent Bernardin          | bernardin@inf.ethz.ch|
  \  MAPLE  /  |Institute for Scientific Computing| Tel: +41 1 632 74 75 |
  <____V____>  |          ETH Zuerich             | Fax: +41 1 632 11 72 |
       |        ---------------------------------------------------------


`simplify/sqrt/fraction`:=proc(f) local i,nu,de,n,d,nr,dr,g,g1,j,ex,a,gtmp,rmvd,tmp;

option `Copyright 1993 by the University of Waterloo`;
n:=traperror(numer(f));
if n=lasterror then RETURN(f) fi;
d:=denom(f);
if d=1 then RETURN(f) fi;

if type(n,`*`) then n:=[op(n)];
elif type(n,`^`) then n:=[n];
else RETURN(f) fi;

if type(d,`*`) then d:=[op(d)];
elif type(d,`^`) then d:=[d];
else RETURN(f) fi;

nr:=n; dr:=d;
n:=map(proc(x)  if type(x,`^`)
           and type(op(2,x),fraction)
           and op(2,op(2,x))=2
        then x^2 else NULL fi end,n);
d:=map(proc(x)  if type(x,`^`)
           and type(op(2,x),fraction)
           and op(2,op(2,x))=2
        then x^2 else NULL fi end,d);

if nops(d)=0 or nops(n)=0 then RETURN(f) fi;

n:=convert(n,`*`); d:=convert(d,`*`);
g:=traperror(frontend(`simplify/sqrt/GCD`,[n,d]));
if g=lasterror then g:=1 else g:=op(1,g) fi;

if nargs=1 then
    if g=1 then RETURN(f) fi;
    if signum(g)=1 then g;
    elif signum(g)=-1 then g;
    else
        g:=factor(g);
        if not type(g,`*`) then g:=1;
        else
            g:=[op(g)];
            tmp:=map(proc(x) if signum(x)=1 then x
                     elif signum(x)=-1 then x
                     else NULL fi end,g);
            g:=convert(tmp,`*`);
        fi;
    fi;
fi;

if g=1 then RETURN(f) fi;
n:=[];
nu:=[];
for i in nr do
    if type(i,`^`) then
        ex:=op(2,i);
        if type(ex,fraction) and op(2,ex)=2 then
            for j from 1 to op(1,ex) do
                n:=[op(n), op(1,i)];
            od;
        else nu:=[op(nu),i];
        fi;
    else nu:=[op(nu),i];
    fi;
od;
d:=[];
de:=[];
for i in dr do
    if type(i,`^`) then
        ex:=op(2,i);
        if type(ex,fraction) and op(2,ex)=2 then
            for j from 1 to op(1,ex) do
                d:=[op(d), op(1,i)];
            od;
        else de:=[op(de),i];
        fi;
    else de:=[op(de),i];
    fi;
od;

g1:=g;
rmvd:=1;

for j from 1 to nops(n) do
    i:=op(j,n);
    gtmp:=traperror(frontend(`simplify/sqrt/GCD`,[g,i]));
    if gtmp=lasterror then gtmp:=[1]
    else g:=gtmp[2]; a:=gtmp[3];
    fi;
    if gtmp[1]<>1 then
        if signum(gtmp[1])=-1 then
            a:=-a;
        fi;
        n:=subsop(j=a,n);
        rmvd:=rmvd*gtmp[1];
    fi;
od;

g:=rmvd;

for j from 1 to nops(d) do
    i:=op(j,d);
    gtmp:=traperror(frontend(`simplify/sqrt/GCD`,[g,i]));
    if gtmp=lasterror then gtmp:=[1]
    else g:=gtmp[2]; a:=gtmp[3];
    fi;
    if gtmp[1]<>1 then
        if signum(gtmp[1])=-1 then
            a:=-a;
        fi;
        d:=subsop(j=a,d);
    fi;
od;

n:=map(proc(x) sqrt(x) end, n);
d:=map(proc(x) sqrt(x) end, d);

n:=[op(n),op(nu)];
d:=[op(d),op(de)];

convert(n,`*`)/convert(d,`*`);
end:
.
