Jaime Villate | 1 Apr 08:31 2010
Picon

Re: CVS server down

On Wed, 2010-03-31 at 16:49 -0600, Robert Dodier wrote:
> No problem. I'm kind of distracted by other things at the moment
> so how about if we postpone the release branch for this weekend.
> I'm guessing CVS will return before then.
> You can send me a message when you get your stuff committed. 
OK, thanks.
Jaime
Dieter Kaiser | 1 Apr 14:58 2010
Picon

Simplification of exp(x*%i*%pi)

I had a look at the implementation of simplifications for expressions
like exp(x*%i*%pi) where x is an arbitrary expression.

With revision 1.88 of simp.lisp the call of %especial has been disabled
for x a float or bigfloat value to avoid curious results. But %especial
can do a lot of more general simplifications. Unfortunately, these have
been disabled too (my error).

Some simplifications of %especial are problematic, e.g.

exp((2+x)^2*%pi*%i)  -> %e^(%i*%pi*(x+2)^2-4*%i*%pi)
exp((2+x)^10*%pi*%i) -> %e^(%i*%pi*(x+2)^10-1024*%i*%pi)

Now the handling of expressions is inconsistent, e.g. the following is
simplified

   exp((2+x)*%i*%pi) -> exp(%i*%pi*x)

but not the expanded form

   exp(2*%i*%pi%+%i*%pi*x) -> e^(%i*%pi*x+2*%i*%pi)

Furthermore, we have inconsistent simplifications of expressions which
contain float or bigfloat values.

The following shows some examples of a testfile to check the
simplification of expressions like exp(x*%pi*%) and the desired results:

********************** Problem 5 ***************
Input:
(Continue reading)

Stavros Macrakis | 1 Apr 16:01 2010
Picon

Re: Simplification of exp(x*%i*%pi)

Dieter,

As a rule, the general simplifier does not make radical changes to the form of expressions.  If the user inputs sinh(x)+cosh(x) or (x^100-1)/(x-1), s/he presumably wants that form for some reason, and the general simplifier does not automatically convert this to the forms %e^x (available through the explicit sequence of exponentialize, ratsimp) or x^99 + x^98 + ... + x + 1 (expand).  It doesn't even convert (x^100-x^99)/(x-1) to x^99 (ratsimp).  This is not an implementation oversight or an over-concern for efficiency (though that is part of it), but an explicit philosophy -- see Moses' paper "Algebraic Simplification: A Guide for the Perplexed" (http://groups.csail.mit.edu/mac/users/gjs/6.945/readings/simplification/moses-simp.pdf).

I have not analyzed or tested your code in detail, but the presence of ratcoef and expand in your code make me think that it does in some cases make radical changes, which may not be ideal.

                -s

On Thu, Apr 1, 2010 at 08:58, Dieter Kaiser <drdieterkaiser <at> web.de> wrote:
I had a look at the implementation of simplifications for expressions
like exp(x*%i*%pi) where x is an arbitrary expression.

With revision 1.88 of simp.lisp the call of %especial has been disabled
for x a float or bigfloat value to avoid curious results. But %especial
can do a lot of more general simplifications. Unfortunately, these have
been disabled too (my error).

Some simplifications of %especial are problematic, e.g.

exp((2+x)^2*%pi*%i)  -> %e^(%i*%pi*(x+2)^2-4*%i*%pi)
exp((2+x)^10*%pi*%i) -> %e^(%i*%pi*(x+2)^10-1024*%i*%pi)

Now the handling of expressions is inconsistent, e.g. the following is
simplified

  exp((2+x)*%i*%pi) -> exp(%i*%pi*x)

but not the expanded form

  exp(2*%i*%pi%+%i*%pi*x) -> e^(%i*%pi*x+2*%i*%pi)

Furthermore, we have inconsistent simplifications of expressions which
contain float or bigfloat values.

The following shows some examples of a testfile to check the
simplification of expressions like exp(x*%pi*%) and the desired results:

********************** Problem 5 ***************
Input:
exp(2*%i*%pi)


Result:
1

... Which was correct.

********************** Problem 6 ***************
Input:
exp((x+2)*%i*%pi)


Result:
%e^(%i*%pi*x)

... Which was correct.

********************** Problem 7 ***************
Input:
exp(x*%i*%pi+2*%i*%pi)


Result:
%e^(%i*%pi*x+2*%i*%pi)

This differed from the expected result:
exp(x*%i*%pi)

********************** Problem 8 ***************
Input:
log(exp((x+2)^2*%i*%pi))


Result:
%i*%pi*(x+2)^2-4*%i*%pi

This differed from the expected result:
(x+2)^2*%i*%pi

********************** Problem 9 ***************
Input:
exp(2.0*%i*%pi)


Result:
%e^(2.0*%i*%pi)

This differed from the expected result:
1.0

********************** Problem 10 ***************
Input:
exp((x+2.0)*%i*%pi)


Result:
%e^(%i*%pi*x)

This differed from the expected result:
exp(1.0*x*%i*%pi)

********************** Problem 11 ***************
Input:
exp(x*%i*%pi+2.0*%i*%pi)


Result:
%e^(%i*%pi*x+2.0*%i*%pi)

This differed from the expected result:
exp(1.0*x*%i*%pi)

********************** Problem 12 ***************
Input:
log(exp((x+2.0)^2*%i*%pi))


Result:
%i*%pi*(x+2.0)^2-4*%i*%pi

This differed from the expected result:
(x+2.0)^2*%i*%pi

********************** Problem 13 ***************
Input:
exp(2.0b0*%i*%pi)


Result:
%e^(2.0b0*%i*%pi)

This differed from the expected result:
1.0b0

********************** Problem 14 ***************
Input:
exp((x+2.0b0)*%i*%pi)


Result:
%e^(%i*%pi*x)

This differed from the expected result:
exp(1.0b0*x*%i*%pi)

********************** Problem 15 ***************
Input:
exp(x*%i*%pi+2.0b0*%i*%pi)


Result:
%e^(%i*%pi*x+2.0b0*%i*%pi)

This differed from the expected result:
exp(1.0b0*x*%i*%pi)

********************** Problem 16 ***************
Input:
log(exp((x+2.0b0)^2*%i*%pi))


Result:
%i*%pi*(x+2.0b0)^2-4*%i*%pi

This differed from the expected result:
(x+2.0b0)^2*%i*%pi

********************** Problem 17 ***************
Input:
exp(3*%pi*%i/2)


Result:
-%i

... Which was correct.

********************** Problem 18 ***************
Input:
exp(1.5*%pi*%i)


Result:
%e^(1.5*%i*%pi)

This differed from the expected result:
-1.0*%i

********************** Problem 19 ***************
Input:
exp(1.5b0*%pi*%i)


Result:
%e^(1.5b0*%i*%pi)

This differed from the expected result:
-1.0b0*%i

********************** Problem 20 ***************
Input:
exp((x+3/2)*%pi*%i)


Result:
-%i*%e^(%i*%pi*x)

... Which was correct.

********************** Problem 21 ***************
Input:
exp((x+1.5)*%pi*%i)


Result:
-%i*%e^(%i*%pi*x)

This differed from the expected result:
-%i*exp(1.0*%i*%pi*x)

********************** Problem 22 ***************
Input:
exp((x+1.5b0)*%pi*%i)


Result:
-%i*%e^(%i*%pi*x)

This differed from the expected result:
-%i*exp(1.0b0*%i*%pi*x)


To get all desired simplifications we have to modify the routine %
especial. This could be the new code:

(defun %especial (e)
 (prog (varlist y k kk j ans $%emode $ratprint genvar)
    (let ()
      (unless (setq y (pip ($ratcoef e '$%i))) (return nil))

      ;; Subtract the term y*%i*%pi from the expression e.
      (setq k ($expand (add e (mul -1 '$%pi '$%i y)) 1))

      ;; This is a workaround to get the type (integer, float, or
bigfloat)
      ;; of the expression. kk must evaluate to 1, 1.0, or 1.0b0.
      ;; Furthermore, if e is nonlinear, kk does not simplify to a
number ONE.
      ;; Because of this we do not simplify something like exp((2+x)^2*
%i*%pi)
      (setq kk (div (sub ($expand e) k) (mul '$%i '$%pi y)))

      ;; Return if kk is not an integer or kk is ONE, but y not an
integer
      ;; or a half integral value.
      (if (not (or (integerp kk)
                   (and (onep1 kk)
                        (integerp (add y y)))))
          (return nil))

      (setq j (trigred y))
      (setq ans (spang1 j t)))

    (cond ((among '%sin ans)
           (cond ((equal y j) (return nil))
                 ((zerop1 k)
                  ;; To preverse the type we add k into the expression.
                  (return (power '$%e (mul %p%i (add k j)))))
                 (t
                   ;; To preserve the type we multiply kk into the
result.
                   (return (power '$%e (add (mul kk k) (mul kk %p%i
j))))))))
    (setq y (spang1 j nil))
    ;; To preserve the type we multiply kk into the expression.
    (return (mul (power '$%e (mul kk k)) (add y (mul '$%i ans))))))

Dieter Kaiser


_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Dieter Kaiser | 1 Apr 16:37 2010
Picon

Re: Simplification of exp(x*%i*%pi)

Hello Stavros,

 

Some more explanations:

 

All these simplifications you can see in the code I have posted are already present in the routine %especial. I have not invented much, but I tried to do it better for integer values, expressions and more consistent for floating point and bigfloat values. Therefore, I have added one calculation and one test to get the desired consistent results. The rest is (very) old code. So, nothing is really new.

 

Furthermore, this type of simplification can be switched off setting the flag %emode to false. The default value is true.

 

With revision 1.88 of simp.lisp %especial is no longer called for float and bigfloat as a factor of %i*%pi, but I have not done the changes complete and consistent.

 

I had again a look at this problem. Now I have tried to do it more complete and consistent. To get it more consistent I had to modify %especial.

 

By the way: The “old” code has a big problem with expression like exp((2+x)^n*%pi*%i). Even for moderate integer values of n %especial needs a very long time to return. One of the modifications I have implemented in %especial is not to try to  simplify such an expression.

 

Thank your very much for the paper. I will study it in detail.

 

Dieter Kaiser

 

Von: Stavros Macrakis [mailto:macrakis <at> gmail.com]
Gesendet: Donnerstag, 1. April 2010 16:01
An: Dieter Kaiser
Cc: Maxima
Betreff: Re: [Maxima] Simplification of exp(x*%i*%pi)

 

Dieter,

As a rule, the general simplifier does not make radical changes to the form of expressions.  If the user inputs sinh(x)+cosh(x) or (x^100-1)/(x-1), s/he presumably wants that form for some reason, and the general simplifier does not automatically convert this to the forms %e^x (available through the explicit sequence of exponentialize, ratsimp) or x^99 + x^98 + ... + x + 1 (expand).  It doesn't even convert (x^100-x^99)/(x-1) to x^99 (ratsimp).  This is not an implementation oversight or an over-concern for efficiency (though that is part of it), but an explicit philosophy -- see Moses' paper "Algebraic Simplification: A Guide for the Perplexed" (http://groups.csail.mit.edu/mac/users/gjs/6.945/readings/simplification/moses-simp.pdf).

I have not analyzed or tested your code in detail, but the presence of ratcoef and expand in your code make me think that it does in some cases make radical changes, which may not be ideal.

                -s

On Thu, Apr 1, 2010 at 08:58, Dieter Kaiser <drdieterkaiser <at> web.de> wrote:

I had a look at the implementation of simplifications for expressions
like exp(x*%i*%pi) where x is an arbitrary expression.

With revision 1.88 of simp.lisp the call of %especial has been disabled
for x a float or bigfloat value to avoid curious results. But %especial
can do a lot of more general simplifications. Unfortunately, these have
been disabled too (my error).

Some simplifications of %especial are problematic, e.g.

exp((2+x)^2*%pi*%i)  -> %e^(%i*%pi*(x+2)^2-4*%i*%pi)
exp((2+x)^10*%pi*%i) -> %e^(%i*%pi*(x+2)^10-1024*%i*%pi)

Now the handling of expressions is inconsistent, e.g. the following is
simplified

  exp((2+x)*%i*%pi) -> exp(%i*%pi*x)

but not the expanded form

  exp(2*%i*%pi%+%i*%pi*x) -> e^(%i*%pi*x+2*%i*%pi)

Furthermore, we have inconsistent simplifications of expressions which
contain float or bigfloat values.

The following shows some examples of a testfile to check the
simplification of expressions like exp(x*%pi*%) and the desired results:

********************** Problem 5 ***************
Input:
exp(2*%i*%pi)


Result:
1

... Which was correct.

********************** Problem 6 ***************
Input:
exp((x+2)*%i*%pi)


Result:
%e^(%i*%pi*x)

... Which was correct.

********************** Problem 7 ***************
Input:
exp(x*%i*%pi+2*%i*%pi)


Result:
%e^(%i*%pi*x+2*%i*%pi)

This differed from the expected result:
exp(x*%i*%pi)

********************** Problem 8 ***************
Input:
log(exp((x+2)^2*%i*%pi))


Result:
%i*%pi*(x+2)^2-4*%i*%pi

This differed from the expected result:
(x+2)^2*%i*%pi

********************** Problem 9 ***************
Input:
exp(2.0*%i*%pi)


Result:
%e^(2.0*%i*%pi)

This differed from the expected result:
1.0

********************** Problem 10 ***************
Input:
exp((x+2.0)*%i*%pi)


Result:
%e^(%i*%pi*x)

This differed from the expected result:
exp(1.0*x*%i*%pi)

********************** Problem 11 ***************
Input:
exp(x*%i*%pi+2.0*%i*%pi)


Result:
%e^(%i*%pi*x+2.0*%i*%pi)

This differed from the expected result:
exp(1.0*x*%i*%pi)

********************** Problem 12 ***************
Input:
log(exp((x+2.0)^2*%i*%pi))


Result:
%i*%pi*(x+2.0)^2-4*%i*%pi

This differed from the expected result:
(x+2.0)^2*%i*%pi

********************** Problem 13 ***************
Input:
exp(2.0b0*%i*%pi)


Result:
%e^(2.0b0*%i*%pi)

This differed from the expected result:
1.0b0

********************** Problem 14 ***************
Input:
exp((x+2.0b0)*%i*%pi)


Result:
%e^(%i*%pi*x)

This differed from the expected result:
exp(1.0b0*x*%i*%pi)

********************** Problem 15 ***************
Input:
exp(x*%i*%pi+2.0b0*%i*%pi)


Result:
%e^(%i*%pi*x+2.0b0*%i*%pi)

This differed from the expected result:
exp(1.0b0*x*%i*%pi)

********************** Problem 16 ***************
Input:
log(exp((x+2.0b0)^2*%i*%pi))


Result:
%i*%pi*(x+2.0b0)^2-4*%i*%pi

This differed from the expected result:
(x+2.0b0)^2*%i*%pi

********************** Problem 17 ***************
Input:
exp(3*%pi*%i/2)


Result:
-%i

... Which was correct.

********************** Problem 18 ***************
Input:
exp(1.5*%pi*%i)


Result:
%e^(1.5*%i*%pi)

This differed from the expected result:
-1.0*%i

********************** Problem 19 ***************
Input:
exp(1.5b0*%pi*%i)


Result:
%e^(1.5b0*%i*%pi)

This differed from the expected result:
-1.0b0*%i

********************** Problem 20 ***************
Input:
exp((x+3/2)*%pi*%i)


Result:
-%i*%e^(%i*%pi*x)

... Which was correct.

********************** Problem 21 ***************
Input:
exp((x+1.5)*%pi*%i)


Result:
-%i*%e^(%i*%pi*x)

This differed from the expected result:
-%i*exp(1.0*%i*%pi*x)

********************** Problem 22 ***************
Input:
exp((x+1.5b0)*%pi*%i)


Result:
-%i*%e^(%i*%pi*x)

This differed from the expected result:
-%i*exp(1.0b0*%i*%pi*x)


To get all desired simplifications we have to modify the routine %
especial. This could be the new code:

(defun %especial (e)
 (prog (varlist y k kk j ans $%emode $ratprint genvar)
    (let ()
      (unless (setq y (pip ($ratcoef e '$%i))) (return nil))

      ;; Subtract the term y*%i*%pi from the expression e.
      (setq k ($expand (add e (mul -1 '$%pi '$%i y)) 1))

      ;; This is a workaround to get the type (integer, float, or
bigfloat)
      ;; of the expression. kk must evaluate to 1, 1.0, or 1.0b0.
      ;; Furthermore, if e is nonlinear, kk does not simplify to a
number ONE.
      ;; Because of this we do not simplify something like exp((2+x)^2*
%i*%pi)
      (setq kk (div (sub ($expand e) k) (mul '$%i '$%pi y)))

      ;; Return if kk is not an integer or kk is ONE, but y not an
integer
      ;; or a half integral value.
      (if (not (or (integerp kk)
                   (and (onep1 kk)
                        (integerp (add y y)))))
          (return nil))

      (setq j (trigred y))
      (setq ans (spang1 j t)))

    (cond ((among '%sin ans)
           (cond ((equal y j) (return nil))
                 ((zerop1 k)
                  ;; To preverse the type we add k into the expression.
                  (return (power '$%e (mul %p%i (add k j)))))
                 (t
                   ;; To preserve the type we multiply kk into the
result.
                   (return (power '$%e (add (mul kk k) (mul kk %p%i
j))))))))
    (setq y (spang1 j nil))
    ;; To preserve the type we multiply kk into the expression.
    (return (mul (power '$%e (mul kk k)) (add y (mul '$%i ans))))))

Dieter Kaiser


_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima

 

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
Dieter Kaiser | 1 Apr 17:44 2010
Picon

Re: Bug report ID:635627 - subst([...] is order-dependent

Am Mittwoch, den 31.03.2010, 14:10 -0400 schrieb Joseph Cusumano:
> 
> Hello:
> 
> I admit up front to not having the perspective (or skills) of a
> developer, so maybe I'm missing something fundamental. However, I do
> not understand why this should be considered a bug. I have a fair
> amount of code in Maple (which has a similar property to its
> substitution command) and now Maxima that depends on (or, at the vary
> least, assumes) this feature property.
> 
> I would say at the very least this should be controllable via an
> external variable, e.g. "subst_sequence" equals "serial" or
> "parallel".

Hello Joseph,

thank you for your comment.

Yes, I think it is no problem to get both the parallel substitution and
the serial substitution we already have.

Because, now the serial substitution is the default and the code of
to_poly_solver depends on this behavior we might first introduce the
parallel substitution as an option to the function $substitute.

Dieter Kaiser
Raymond Toy | 1 Apr 17:59 2010
Picon

Re: Bug report ID:635627 - subst([...] is order-dependent

On 4/1/10 11:44 AM, Dieter Kaiser wrote:
> Am Mittwoch, den 31.03.2010, 14:10 -0400 schrieb Joseph Cusumano:
>   
>> Hello:
>>
>> I admit up front to not having the perspective (or skills) of a
>> developer, so maybe I'm missing something fundamental. However, I do
>> not understand why this should be considered a bug. I have a fair
>> amount of code in Maple (which has a similar property to its
>> substitution command) and now Maxima that depends on (or, at the vary
>> least, assumes) this feature property.
>>
>> I would say at the very least this should be controllable via an
>> external variable, e.g. "subst_sequence" equals "serial" or
>> "parallel".
>>     
> Hello Joseph,
>
> thank you for your comment.
>
> Yes, I think it is no problem to get both the parallel substitution and
> the serial substitution we already have.
>
> Because, now the serial substitution is the default and the code of
> to_poly_solver depends on this behavior we might first introduce the
> parallel substitution as an option to the function $substitute.
>   
I didn't follow this dicussion too closely, but why not just leave subst
to do serial substitution, and add, say, psubst to do parallel
substitution? 

I think it's bad to change subst to do parallel because that changes
existing behavior.  There could be lots of code out there that depends
on this and we should gratuitiously break them.  I know this is the
opposite of what Robert said, but I don't see a good reason to break
other people's code in this case.

Ray
dlakelan | 1 Apr 21:31 2010

Lie Groups, Fluid Mechanics, and a Paper on Use of Maxima in symmetry analysis

I was randomly writing a blog entry today when I looked up the "Law of 
The Wall" on Wikipedia. This is a concept in fluid mechanics of 
turbulent flows. The wikipedia article pointed me to a ten year old 
paper on deriving turbulent flow solutions via Lie symmetry analysis.

The paper by Martin Oberlack is here:

http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=64901

It offers a description of a variety of symmetry group software packages 
developed in the late 1990's and then a description of the use of a 
particular package SIMMGROUP.MAX, a version of which is available here 
from Professor Hereman's page.

http://inside.mines.edu/~whereman/software.html

Perhaps someone finds this interesting, and/or would stick it in a 
useful place on the wiki?
dlakelan | 1 Apr 21:46 2010

Re: Lie Groups, Fluid Mechanics, and a Paper on Use of Maxima in symmetry analysis

On 04/01/2010 12:31 PM, dlakelan wrote:

> It offers a description of a variety of symmetry group software packages
> developed in the late 1990's and then a description of the use of a
> particular package SIMMGROUP.MAX, a version of which is available here
> from Professor Hereman's page.

Whoops, slightly mixing up two different papers, the One by Oberlack 
focuses on the application to fluid mechanics, and the one by Hereman here:

http://inside.mines.edu/~whereman/papers/Hereman-AMS-SIAMLNAP-28-1993.pdf

Focuses on the general symmetry calculation software. That's what I get 
for having too many windows open at once:
Barton Willis | 1 Apr 22:06 2010

Re: Bug report ID:635627 - subst([...] is order-dependent

> I think it's bad to change subst to do parallel because that changes
> existing behavior.  There could be lots of code out there that depends
> on this and we should gratuitiously break them. 

I agree. My enthusiasm for breaking existing code is low. But the user
documentation for subst doesn't specify the substitution order (or
simplification order), so there is an argument for saying that code
that depends on a specific substitution order is buggy.

--Barton
Dieter Kaiser | 1 Apr 22:32 2010
Picon

Re: Bug report ID:635627 - subst([...] is order-dependent

Am Donnerstag, den 01.04.2010, 15:06 -0500 schrieb Barton Willis:
> > I think it's bad to change subst to do parallel because that changes
> > existing behavior.  There could be lots of code out there that depends
> > on this and we should gratuitiously break them. 
> 
> I agree. My enthusiasm for breaking existing code is low. But the user
> documentation for subst doesn't specify the substitution order (or
> simplification order), so there is an argument for saying that code
> that depends on a specific substitution order is buggy.

When I summarize the different opinions I think it is the best not to
change the current implementation of subst, to add some documentation
about the way substitute works, and to close the open bug report as
"won't fix".

Perhaps, we have the following additional options:

1. 
Adding a function like psubstitute for parallel substitution which has a
functionality comparable to the function substitute.

2. 
Adding a flag to the existing function substitute to switch on parallel
substitution.

Remark: I think it is not a good idea to try to extend sublis to get a
functionality like substitute, because we have to double a lot of code.

Dieter Kaiser

Gmane