[Axiom-developer] Re: [Axiom-mail] comparison operators and %pi

[Martin Rubey]

Dear Alasdair, Francois,

(didn't get through, did it?)

"Alasdair McAndrew" <address@hidden> writes:

(I move this to developer)

> I also notice that pi is defined as an operator, so that pi() returns the
> result %pi.  I wonder how hard it would be to include useful constants in
> Axiom, and to be able to enter (for example)
> 
> solve(integrate(1/x,x=1..t)=1,t)
> 
> and obtain the answer %e.

(6) -> solve(integrate(1/x,x=1..t, "noPole")::EXPR INT = 1, t)

               +---+     +---+
               |  2      |  2
   (6)  [t= - \|%e  ,t= \|%e  ]
                                       Type: List Equation Expression Integer

Yes, we should teach Axiom that %e is positive.  In fact, I don't think that's
too difficult in a hackish way.  In expr.spad, you find

  CF  ==> CombinatorialFunction(R, %)

      x:% ** y:%                == x **$CF y

In our case, x is %e^2 and y is 1/2.  The case of integer y is dealt with
seperately in expr.spad, although maybe this should be changed.  (Note that
algorithms in CF need to hold in greater generality, namely for all
"FunctionSpace"s.  But at the moment I do not see whether the special
definitions can be moved to CF or not, see below for them.

In CF (combfunc.spad) you find

    x ** y               == 
      -- Do some basic simplifications
      is?(x,POWER) =>
        args : List F := argument first kernels x
        not(#args = 2) => error "Too many arguments to **"
        number?(first args) and number?(y) =>
          oppow [first(args)**y, second args]
        oppow [first args, (second args)* y]
      -- Generic case
      exp : Union(Record(val:F,exponent:Z),"failed") := isPower x
      exp case Record(val:F,exponent:Z) =>
        expr := exp::Record(val:F,exponent:Z)
        oppow [expr.val, (expr.exponent)*y]
      oppow [x, y]

Of course, you cannot simply check whether x is positive, since the ordering is
not the mathematical one.  So we can only deal with special cases for the
moment (unless we get provisos).  I think there is only one special case that
will somewhat work: check whether x is numerical, convert it to float and check
whether this float is positive.

Note that x**y might be of the form sqrt(-1)^2.  Thus, even if x is a square,
it is not necessarily positive.

Francois, you have already some expertise.  Maybe you can help here.  Although,
in reality, we need a redesign.

Martin


Here are the special cases from EXPR.  

      algkernels l     == select_!(has?(operator #1, ALGOP), l)

      simplifyPower(x:%,n:Integer):% ==
        k : List K := kernels x
        is?(x,POWER) =>
          -- Look for a power of a number in case we can do a simplification
          args : List % := argument first k
          not(#args = 2) => error "Too many arguments to **"
          number?(args.1) =>
             reduc((args.1) **$Rep n, algkernels kernels (args.1))**(args.2)
          (first args)**(n*second(args))
        reduc(x **$Rep n, algkernels k)

Hm, no idea what reduc does.

      reduc(x, l) ==
        for k in l repeat
          p := minPoly k
          x := evl(numer x, k, p) /$Rep evl(denom x, k, p)
        x

      evl0(p, k) ==
        numer univariate(p::Fraction(MP),
                     k)$PolynomialCategoryQuotientFunctions(IndexedExponents K,
                                                            K,R,MP,Fraction MP)

      -- uses some operations from Rep instead of % in order not to
      -- reduce recursively during those operations.
      evl(p, k, m) ==
        degree(p, k) < degree m => p::Fraction(MP)
        (((evl0(p, k) pretend SparseUnivariatePolynomial($)) rem m)
           pretend SparseUnivariatePolynomial Fraction MP)
      -- (k::MP::Fraction(MP))

In any case, we cannot simply copy-paste it to CF, since we are using special
features from EXPR here...  Some documentation would be good, though.

      x:% ** n:NonNegativeInteger ==
        n = 0 => 1$%
        n = 1 => x
        simplifyPower(numerator x,n pretend Integer) /
        simplifyPower(denominator x,n pretend Integer)

      x:% ** n:Integer ==
        n = 0 => 1$%
        n = 1 => x
        n = -1 => 1/x
        simplifyPower(numerator x,n) / simplifyPower(denominator x,n)

      x:% ** n:PositiveInteger == 
        n = 1 => x
        simplifyPower(numerator x,n pretend Integer) /
        simplifyPower(denominator x,n pretend Integer)