[Axiom-math] Axiom's integration non-deterministic?

[Igor Khavkine]

Start up Axiom, and type in the following integral

integrate(sqrt(1+x^(-2/3)),x)

Then type in the same integral again. Maybe try the same integral a
few more times. Do you get the same answer each time?

I don't. I see the same behavior from axiom-20050901-4.1 (Debian sid)
and axiom-20050201-1 (Debian sarge).

I get two possible answers, one is (x^(3/2)+1)*sqrt(x^(3/2)+1) and the
other one is a much more complicated radical expression. Ironically,
the longer expression takes less time to compute. Differentiating both
answers gets me back to the original integrand. Also, their curves
coincide when plotted. It seems that each time the answer is correct.
However, the shorter one is obviously more attractive because it looks
simpler.

Clearly, Axiom takes two different paths through the integration
algorithm, even when give identical input. What is the cause of the
branch? Is there a non-deterministic step somewhere in the algorithm?

Incidentally, is there a canonical form for radical expressions in
which the two forms of the answer can be compared and directly shown
to be the same?

Thanks in advance.

Igor

[root]

fascinating.

notice that the return types are different also.

the short version is Union(Expression Integer,...)
the long version is Expression Integer

t

[Igor Khavkine]

Post by Igor Khavkine
Clearly, Axiom takes two different paths through the integration
algorithm, even when give identical input. What is the cause of
the branch? Is there a non-deterministic step somewhere in the
algorithm?

I don't know the cause of this behaviour but I get the same effect
on the Windows version of Axiom.

So there must be something strage going on in the integration code.
I've been trying to trace the algorithm's execution to see what it
does exactly in each case. Unfortunately, I'm having little luck. I
played with a few variations of )trace command, but couldn't find the
right one.

I've identified the routine 'integrate' in the FSINT package. Knowing
that, is there a way for me to step through its code line by line like
in a usual debugger?

Thanks.

Igor

[root]

you can start axiom and type:

)trace INTEF )math

which will trace external functions but

the INTEF domain compiles into lisp.
the lisp code lives in int/algebra/INTEF.NRLIB/code.lsp
you can look at that file and it contains DEFUNs
you can trace those functions. For example,

in int/algebra/INTEF.NRLIB/code.lsp the function prim? is defined as:

(DEFUN |INTEF;prim?| (|k| |x| $)
(COND
((SPADCALL |k| (SPADCALL "log" (QREFELT $ 10)) (QREFELT $ 13))
(QUOTE T))
((QUOTE T)
(SPADCALL (SPADCALL |k| (QREFELT $ 15)) "prim" (QREFELT $ 16)))))

)lisp (trace |INTEF;prim?|)

which will give you the lisp level trace of the function.

t

[root]

Igor,

This is a specific example of tracing applied to INTEF.
The symbols used for the lisp trace commands came from
int/algebra/INTEF.NRLIB/code.lsp

)trace INTEF )math
)lisp (trace |INTEF;prim?|)
)lisp (trace |INTEF;tanint|)
)lisp (trace |INTEF;tanint!5|)
)lisp (trace |INTEF;tanint!4|)
)lisp (trace |INTEF;tanint!3|)
)lisp (trace |INTEF;tanint!2|)
)lisp (trace |INTEF;tanint!1|)
)lisp (trace |INTEF;tanint!0|)
)lisp (trace |INTEF;unknownint|)
)lisp (trace |INTEF;droponex|)
)lisp (trace |INTEF;unklimint|)
)lisp (trace |INTEF;unkextint|)
)lisp (trace |INTEF;isx?|)
)lisp (trace |INTEF;alglfint|)
)lisp (trace |INTEF;alglfint!0|)
)lisp (trace |INTEF;alglfextint|)
)lisp (trace |INTEF;alglflimint|)
)lisp (trace |INTEF;alglflimint!0|)
)lisp (trace |INTEF;lfintegrate;FSIr;11|)
)lisp (trace |INTEF;lfintegrate;FSIr;12|)
)lisp (trace |INTEF;lfintegrate;FSIr;13|)
)lisp (trace |INTEF;lfintegrate;FSIr;14|)
)lisp (trace |INTEF;lfintegrate0|)
)lisp (trace |INTEF;lfintegrate0!0|)
)lisp (trace |INTEF;addx|)
)lisp (trace |INTEF;tryChangeVar|)
)lisp (trace |INTEF;tryChangeVar!0|)
)lisp (trace |INTEF;algexpint|)
)lisp (trace |INTEF;algexpint!1|)
)lisp (trace |INTEF;algexpint!0|)
)lisp (trace |INTEF;algprimint|)
)lisp (trace |INTEF;algprimint!1|)
)lisp (trace |INTEF;algprimint!0|)
)lisp (trace |INTEF;lfextendedint;FSFU;20|)
)lisp (trace |INTEF;lfextendedint;FSFU;20!0|)
)lisp (trace |INTEF;lflimitedint;FSLU;21|)
)lisp (trace |INTEF;lflimitedint;FSLU;21!0|)
)lisp (trace |INTEF;lfinfieldint;FSU;22|)
)lisp (trace |INTEF;primextint|)
)lisp (trace |INTEF;primextint!2|)
)lisp (trace |INTEF;primextint!1|)
)lisp (trace |INTEF;primextint!0|)
)lisp (trace |INTEF;expextint|)
)lisp (trace |INTEF;expextint!4|)
)lisp (trace |INTEF;expextint!3|)
)lisp (trace |INTEF;expextint!2|)
)lisp (trace |INTEF;expextint!1|)
)lisp (trace |INTEF;expextint!0|)
)lisp (trace |INTEF;primint|)
)lisp (trace |INTEF;primint!3|)
)lisp (trace |INTEF;primint!2|)
)lisp (trace |INTEF;primint!1|)
)lisp (trace |INTEF;primint!0|)
)lisp (trace |INTEF;lfextlimint;FSKLU;26|)
)lisp (trace |INTEF;cfind|)
)lisp (trace |INTEF;expint|)
)lisp (trace |INTEF;expint!5|)
)lisp (trace |INTEF;expint!4|)
)lisp (trace |INTEF;expint!3|)
)lisp (trace |INTEF;expint!2|)
)lisp (trace |INTEF;expint!1|)
)lisp (trace |INTEF;expint!0|)
)lisp (trace |INTEF;primlimint|)
)lisp (trace |INTEF;primlimint!2|)
)lisp (trace |INTEF;primlimint!1|)
)lisp (trace |INTEF;primlimint!0|)
)lisp (trace |INTEF;explimint|)
)lisp (trace |INTEF;explimint!4|)
)lisp (trace |INTEF;explimint!3|)
)lisp (trace |INTEF;explimint!2|)
)lisp (trace |INTEF;explimint!1|)
)lisp (trace |INTEF;explimint!0|)
)lisp (trace |ElementaryIntegration|)
)lisp (trace |ElementaryIntegration;|)


t

[William Sit]

Post by Igor Khavkine
Start up Axiom, and type in the following integral
integrate(sqrt(1+x^(-2/3)),x)
Then type in the same integral again. Maybe try the same integral a
few more times. Do you get the same answer each time?
I don't. I see the same behavior from axiom-20050901-4.1 (Debian sid)
and axiom-20050201-1 (Debian sarge).
I get two possible answers, one is (x^(3/2)+1)*sqrt(x^(3/2)+1) and the
other one is a much more complicated radical expression. Ironically,
the longer expression takes less time to compute. Differentiating both
answers gets me back to the original integrand. Also, their curves
coincide when plotted. It seems that each time the answer is correct.
However, the shorter one is obviously more attractive because it looks
simpler.
Clearly, Axiom takes two different paths through the integration
algorithm, even when give identical input. What is the cause of the
branch? Is there a non-deterministic step somewhere in the algorithm?
Incidentally, is there a canonical form for radical expressions in
which the two forms of the answer can be compared and directly shown
to be the same?
Thanks in advance.
Igor

I vaguely remembered that Axiom does caching. So the first time the command
above is done, it gives an unsimplified answer and the second time it may be
doing the simplification using the previous answer. Simplifying radical
expresssion is expensive (I think it requires constructing towers of algebraic
extensions). The two answers are identical (on my Windows version)

AXIOM Computer Algebra System
Version of Tuesday November 30, 2004 at 21:11:14
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave AXIOM and return to shell.
-----------------------------------------------------------------------------

(1) -> a:=integrate(sqrt(1+x^(-2/3)),x)

+---------+
3+-+2 3+-+ |3+-+2 3+-+2 3+-+ 2
(4x \|x + 3\|x + 7x)\|\|x + 1 + 6\|x + 9x\|x + 4x + 1
(1) ---------------------------------------------------------------
+---------+
3+-+2 |3+-+2 3+-+
(4\|x + 1)\|\|x + 1 + 3\|x + 4x
Type: Union(Expression Integer,...)
(2) -> b:=integrate(sqrt(1+x^(-2/3)),x)

+---------+
3+-+2 |3+-+2
(2) (\|x + 1)\|\|x + 1
Type: Union(Expression Integer,...)
(3) -> c:=integrate(sqrt(1+x^(-2/3)),x)

+---------+
3+-+2 |3+-+2
(3) (\|x + 1)\|\|x + 1
Type: Union(Expression Integer,...)
(4) -> t:Boolean:=(a=b)

(4) false
Type: Boolean
(5) -> a-b

(5) 0
Type: Expression Integer
(6) -> a-c

(6) 0
Type: Expression Integer

Two comments: (2) and (3) took longer than (1). (4) is not compatible with (5)
or (6) (the reason may be because the domain is a Union, and equality there may
be by data representation. In general, equality and zero testing are different.)

William

[Igor Khavkine]

Post by William Sit
I vaguely remembered that Axiom does caching. So the first time the command
above is done, it gives an unsimplified answer and the second time it may be
doing the simplification using the previous answer. Simplifying radical
expresssion is expensive (I think it requires constructing towers of algebraic
extensions). The two answers are identical (on my Windows version)

The first time I tried this expression, it was the simplified answer
that axiom spat out first, the more complicated answer second, which
then alternated (but without any kind of regular pattern I could
recognize). So it doesn't look like any kind of caching is
responsible, at least in version axiom-20050201. In the 20050901
version, the behavior seems closer to what William found.

Strange in deed.

Igor

[William Sit]

Here is more detail, including loading and timing:

AXIOM Computer Algebra System
Version of Tuesday November 30, 2004 at 21:11:14
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave AXIOM and return to shell.
-----------------------------------------------------------------------------

(1) -> )set mess time on
(1) -> a:=integrate(sqrt(1+x^(-2/3)),x)
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/EXPR.o for domain
Expression
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/KERNEL.o for
domain Kernel
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/BOP.o for domain
BasicOperator
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/SCACHE.o for
package SortedCache
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/VOID.o for domain
Void
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/UPMP.o for package
UnivariatePolynomialMultiplicationPackage
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/AN.o for domain
AlgebraicNumber
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IAN.o for domain
InnerAlgebraicNumber
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ACFS-.o for domain
AlgebraicallyClosedFunctionSpace&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FS-.o for domain
FunctionSpace&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/COMMONOP.o for
package CommonOperators
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/KDAGG-.o for
domain KeyedDictionary&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/DIAGG-.o for
domain Dictionary&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/BOP1.o for package
BasicOperatorFunctions1
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ACF-.o for domain
AlgebraicallyClosedField&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ES-.o for domain
ExpressionSpace&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FRETRCT-.o for
domain FullyRetractableTo&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/EVALAB-.o for
domain Evalable&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/TRANFUN-.o for
domain TranscendentalFunctionCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/BASTYPE-.o for
domain BasicType&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IEVALAB-.o for
domain InnerEvalable&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/TRIGCAT-.o for
domain TrigonometricFunctionCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ATRIG-.o for
domain ArcTrigonometricFunctionCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/HYPCAT-.o for
domain HyperbolicFunctionCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ELEMFUN-.o for
domain ElementaryFunctionCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/RADCAT-.o for
domain RadicalCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/AF.o for package
AlgebraicFunction
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/POLYCATQ.o for
package PolynomialCategoryQuotientFunctions
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/POLYROOT.o for
package PolynomialRoots
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/PGCD.o for package
PolynomialGcdPackage
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FSAGG-.o for
domain FiniteSetAggregate&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FLASORT.o for
package FiniteLinearAggregateSort
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FSINT.o for
package FunctionSpaceIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/TRIGMNIP.o for
package TrigonometricManipulations
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTFACT.o for
package IntegerFactorizationPackage
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IROOT.o for
package IntegerRoots
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FACTFUNC.o for
package FactoredFunctions
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/EF.o for package
ElementaryFunction
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IDPOAMS.o for
domain IndexedDirectProductOrderedAbelianMonoidSup
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IDPOAM.o for
domain IndexedDirectProductOrderedAbelianMonoid
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/SIGNEF.o for
package ElementaryFunctionSign
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MULTSQFR.o for
package MultivariateSquareFree
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/UPSQFREE.o for
package UnivariatePolynomialSquareFree
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/HEUGCD.o for
package HeuGcd
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INMODGCD.o for
package InnerModularGcd
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/EMR.o for domain
EuclideanModularRing
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MDDFACT.o for
package ModularDistinctDegreeFactorizer
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MODRING.o for
domain ModularRing
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/TOOLSIGN.o for
package ToolsForSign
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/EFSTRUC.o for
package ElementaryFunctionStructurePackage
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/POLYLIFT.o for
package PolynomialCategoryLifting
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTTOOLS.o for
package IntegrationTools
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ALGMANIP.o for
package AlgebraicManipulations
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/LIST2MAP.o for
package ListToMap
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTEF.o for
package ElementaryIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FSPRMELT.o for
package FunctionSpacePrimitiveElement
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/PRIMELT.o for
package PrimitiveElement
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/PRS.o for package
PseudoRemainderSequence
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FACUTIL.o for
package FactoringUtilities
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTPAF.o for
package PureAlgebraicIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTG0.o for
package GenusZeroIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/UPCDEN.o for
package UnivariatePolynomialCommonDenominator
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/CDEN.o for package
CommonDenominator
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTRAT.o for
package RationalIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTTR.o for
package TranscendentalIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTHERTR.o for
package TranscendentalHermiteIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MONOTOOL.o for
package MonomialExtensionTools
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IR.o for domain
IntegrationResult
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IR2.o for package
IntegrationResultFunctions2
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IR2F.o for package
IntegrationResultToFunction
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/COMPLEX.o for
domain Complex

+---------+
3+-+2 3+-+ |3+-+2 3+-+2 3+-+ 2
(4x \|x + 3\|x + 7x)\|\|x + 1 + 6\|x + 9x\|x + 4x + 1
(1) ---------------------------------------------------------------
+---------+
3+-+2 |3+-+2 3+-+
(4\|x + 1)\|\|x + 1 + 3\|x + 4x
Type: Union(Expression Integer,...)
Time: 0.02 (IN) + 0.17 (EV) + 0.42 (OT) + 0.12 (GC) = 0.72 sec
(2) ->
(2) -> b:=integrate(sqrt(1+x^(-2/3)),x)
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/CHVAR.o for
package ChangeOfVariable
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ALGFF.o for domain
AlgebraicFunctionField
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MATRIX.o for
domain Matrix
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IIARRAY2.o for
domain InnerIndexedTwoDimensionalArray
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MATCAT-.o for
domain MatrixCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ARR2CAT-.o for
domain TwoDimensionalArrayCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTALG.o for
package AlgebraicIntegrate
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/SAE.o for domain
SimpleAlgebraicExtension
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/INTHERAL.o for
package AlgebraicHermiteIntegration
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FFINTBAS.o for
package FunctionFieldIntegralBasis
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MONOGEN-.o for
domain MonogenicAlgebra&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MATLIN.o for
package MatrixLinearAlgebraFunctions
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MHROWRED.o for
package ModularHermitianRowReduction
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/TRIMAT.o for
package TriangularMatrixOperations
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IBATOOL.o for
package IntegralBasisTools
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/VECTCAT-.o for
domain VectorCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FRAMALG-.o for
domain FramedAlgebra&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FFCAT-.o for
domain FunctionFieldCategory&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MATCAT2.o for
package MatrixCategoryFunctions2
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/IMATLIN.o for
package InnerMatrixLinearAlgebraFunctions
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FINRALG-.o for
domain FiniteRankAlgebra&
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/MCDEN.o for
package MatrixCommonDenominator
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/LIST2.o for
package ListFunctions2
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/FLAGG2.o for
package FiniteLinearAggregateFunctions2
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/ICDEN.o for
package InnerCommonDenominator
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/VECTOR2.o for
package VectorFunctions2
Loading j:/OpenAxiom/axiom014/mnt/windows/algebra/LSMP.o for package
LinearSystemMatrixPackage

+---------+
3+-+2 |3+-+2
(2) (\|x + 1)\|\|x + 1
Type: Union(Expression Integer,...)
Time: 0.02 (IN) + 6.18 (EV) + 0.12 (OT) + 0.78 (GC) = 7.10 sec
(3) ->
(3) -> c:=integrate(sqrt(1+x^(-2/3)),x)

+---------+
3+-+2 |3+-+2
(3) (\|x + 1)\|\|x + 1
Type: Union(Expression Integer,...)
Time: 6.40 (EV) + 0.45 (GC) = 6.85 sec
(4) -> a-b

(4) 0
Type: Expression Integer
Time: 0 sec
(5) -> interpret(a::InputForm) - interpret (b::InputForm)

(5) 0
Type: Expression Integer
Time: 0.03 (EV) + 0.03 (OT) = 0.07 sec
(6) -> t:Boolean:=(a::EXPR INT)=(b::EXPR INT)

(6) false
Type: Boolean
Time: 0 sec
(7) -> d:=integrate(sqrt(1+y^(-2/3)),y)

+---------+
3+-+2 |3+-+2
(7) (\|y + 1)\|\|y + 1
Type: Union(Expression Integer,...)
Time: 0.02 (IN) + 6.98 (EV) + 0.48 (GC) = 7.48 sec

I am no longer sure why (3) took almost as long as (2) (in fact longer EV time),
and why (7) took even longer EV time. However, the longer time compared to (1)
is explained by the extra domains that are loaded. So the system behavior (or at
least the integration algorithm) is affected by what domains are already loaded
(that explains why (7) is comparable to (2) and (3), not (1)). (6) shows that
equality in EXPR INT is not equivalent to zero testing of the difference. Note
that the domain CHVAR was the first loaded in (2), which enabled simplification
of radical expressions. So to answer Igor's question, we way try tracing what
caused CHVAR to be loaded.

William

[William Sit]

More intrique. But I think I am getting close to the reason (which is NOT what I
said before -- those were wrong guesses). But, first, some more "strange"
behavior that led me to the chase. Please bear with me.

The integrate(sqrt(1+x^(-2/3)),x) command is run by selecting integrate from
FSINT(INT, EXPR INT), which when applied to the given function, I believe is
equivalent to:

integrate(f, x) ==
not real? f => complexIntegrate(f, x)
-- [... ] irrelevant for this example
rec := rischNormalize(realElementary(f, x), x)
g := rootSimp(rec.func) -- same as g:=f for this example
--- [... ] more irrelevant code
i:IR
if (comp := goComplex?(rtg, tg, ltg)) then -- false for this example
--- [... ] irrelevant code
else i := lfintegrate(g, x) -- this is what should be executed
--- [... ]
(u := rinteg(i, f, x, el and ht, comp)) case F => --- true, u := i
postSubst(u::F, rec.vals, rec.kers, comp, ltg, x) --- output u
--- [... ] irrelevant code

When the system first started, as we know:

(1) -> a:=integrate(sqrt(1+x^(-2/3)),x)
+---------+
3+-+2 3+-+ |3+-+2 3+-+2 3+-+ 2
(4x \|x + 3\|x + 7x)\|\|x + 1 + 6\|x + 9x\|x + 4x + 1
(1) ---------------------------------------------------------------
+---------+
3+-+2 |3+-+2 3+-+
(4\|x + 1)\|\|x + 1 + 3\|x + 4x
Type: Union(Expression Integer,...)

comparing this with (new session):

(1) -> complexIntegrate(sqrt(1+x^(-2/3)),x)

+---------+
3+-+2 3+-+ |3+-+2 3+-+2 3+-+ 2
(4x \|x + 3\|x + 7x)\|\|x + 1 + 6\|x + 9x\|x + 4x + 1
(1) ---------------------------------------------------------------
+---------+
3+-+2 |3+-+2 3+-+
(4\|x + 1)\|\|x + 1 + 3\|x + 4x
Type: Expression Integer
(2) -> complexIntegrate(sqrt(1+x^(-2/3)),x)

+---------+
3+-+2 |3+-+2
(2) (\|x + 1)\|\|x + 1
Type: Expression Integer

On the other hand (new session):
(1) -> real? sqrt(1+x^(-2/3))

(1) true
Type: Boolean

so the first branch in 'integrate' should simply fall through! It is not clear
whether 'true' is correct since sqrt(1+x^(-2/3)) can take on complex values. But
that does not seem to be the reason either.

The simplified integral can be done by (new session):

(1) -> q:=sqrt(1+x^(-2/3));
(2) -> p:=lfintegrate(q,x::Symbol)$ElementaryIntegration(INT, EXPR INT)

+---------+
|3+-+2
3+-+2 |\|x + 1
(x \|x + x) |---------
| 3+-+2
\| \|x
(2) -------------------------
3+-+2
\|x
Type: IntegrationResult Expression Integer
(3) -> rootSimp p

+---------+
3+-+2 |3+-+2
(3) (\|x + 1)\|\|x + 1
Type: Expression Integer
(4) -> a:=integrate(sqrt(1+x^(-2/3)),x)
+---------+
3+-+2 |3+-+2
(4) (\|x + 1)\|\|x + 1
Type: Union(Expression Integer,...)

Now (4) gave the simplifed answer (which by tracing integrate from FSINT(INT,
EXPR INT), as far as I can tell, follows (1), (2). So does
complexIntegrate(sqrt(1+x^(-2/3)),x). In fact, even tracing complexIntegrate,
which leads to internalIntegrate and internalIntegrate0 will lead back to
lfintegrate for this integrand. It is no longer clear whether loading CHVAR has
anything to do with the simplification. (I used ')set expose add constructor
CHVAR' at the beginning and it does not seem to help.) More digging.

Further investigation with )trace shows that the intermediate steps involve
different integrands (calls to lfintegrate).

AXIOM Computer Algebra System
Version of Tuesday November 30, 2004 at 21:11:14
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave AXIOM and return to shell.
-----------------------------------------------------------------------------

(1) -> )set mess autoload off
(1) -> )trace lfintegrate
lfintegrate is not a function

Nothing is traced now.

(1) -> )trace lfintegrate$INTEF(INT, EXPR INT)

AXIOM does not understand the use of $elt(INTEF INT (EXPR INT))
lfintegrate here.
(1) -> )trace INTEF )math

Parameterized constructors traced:
INTEF
(1) -> q:=sqrt(1+x^(-2/3))

+---------+
|3+-+2
|\|x + 1
(1) |---------
| 3+-+2
\| \|x
Type: Expression Integer
(2) -> complexIntegrate(q,x)
1<enter ElementaryIntegration.lfintegrate,58 :
+---------+
|3+-+2
\|\|x + 1
arg1= ------------
3+-+
\|x
arg2= x
2<enter ElementaryIntegration.lfintegrate,58 :
4 2
- %F + 3%F + 8x %F - 4
arg1= ------------------------
4 2
%F - 3%F - 8x %F + 2
arg2= x
2>exit ElementaryIntegration.lfintegrate,58 :
4 3
3%F + 4x %F + 1
-----------------
3
4%F
1>exit ElementaryIntegration.lfintegrate,58 :
+---------+
3+-+2 3+-+ |3+-+2 3+-+2 3+-+ 2
(4x \|x + 3\|x + 7x)\|\|x + 1 + 6\|x + 9x\|x + 4x + 1
---------------------------------------------------------------
+---------+
3+-+2 |3+-+2 3+-+
(4\|x + 1)\|\|x + 1 + 3\|x + 4x

+---------+
3+-+2 3+-+ |3+-+2 3+-+2 3+-+ 2
(4x \|x + 3\|x + 7x)\|\|x + 1 + 6\|x + 9x\|x + 4x + 1
(2) ---------------------------------------------------------------
+---------+
3+-+2 |3+-+2 3+-+
(4\|x + 1)\|\|x + 1 + 3\|x + 4x
Type: Expression Integer
(3) -> complexIntegrate(q,x)
1<enter ElementaryIntegration.lfintegrate,58 :
+---------+
|3+-+2
\|\|x + 1
arg1= ------------
3+-+
\|x
arg2= x
2<enter ElementaryIntegration.lfintegrate,58 :
5 4 3 2 2
- 3%O - 32x %O + 63%O + 840x %O + (4352x - 486)%O - 4032x
arg1= -----------------------------------------------------------------
5 4 3 2 2
9%O + 96x %O - 189%O - 2520x %O + (- 3840x + 972)%O + 12096x
arg2= x
2>exit ElementaryIntegration.lfintegrate,58 :
5 4 2 3 2
- 1455%O - 15520x %O + (49664x + 27936)%O + 407400x %O
+
2 3 2
(769792x - 164997)%O - 1390592x + 3454976x - 1882188x - 182196
/
2
1787904x - 94284
1>exit ElementaryIntegration.lfintegrate,58 :
+---------+
3+-+2 |3+-+2
(873\|x + 873)\|\|x + 1 + 1687
-----------------------------------
873

+---------+
3+-+2 |3+-+2
(3) (\|x + 1)\|\|x + 1
Type: Expression Integer
(4) -> complexIntegrate(q,x)
1<enter ElementaryIntegration.lfintegrate,58 :
+---------+
|3+-+2
\|\|x + 1
arg1= ------------
3+-+
\|x
arg2= x
2<enter ElementaryIntegration.lfintegrate,58 :
5 4 3 2 2
224%W + 441x %W - 2304%W - 10360x %W + (- 18081x + 9728)%W + 19536
x
arg1= ------------------------------------------------------------------------
-
5 4 3 2 2
448%W + 882x %W - 4608%W - 20720x %W + (- 17640x + 11264)%W + 39072
x
arg2= x
2>exit ElementaryIntegration.lfintegrate,58 :
5 4 2 3 2
495840%W + 976185x %W + (- 1366659x - 4495616)%W - 22932600x %W
+
2 3 2
(- 23623677x + 14280192)%W + 17766567x - 95304951x + 35386448x + 42151936

/
2
19133226x - 8462336
1>exit ElementaryIntegration.lfintegrate,58 :
+---------+
3+-+2 |3+-+2
(2066\|x + 2066)\|\|x + 1 - 10291
--------------------------------------
2066

+---------+
3+-+2 |3+-+2
(4) (\|x + 1)\|\|x + 1
Type: Expression Integer
(5) -> integrate(q,x)
1<enter ElementaryIntegration.lfintegrate,58 :
+---------+
|3+-+2
\|\|x + 1
arg1= ------------
3+-+
\|x
arg2= x
2<enter ElementaryIntegration.lfintegrate,58 :
5 4 3 2 2
- 3%BE - 32x %BE + 63%BE + 840x %BE + (4352x - 486)%BE - 4032x
arg1= ----------------------------------------------------------------------
5 4 3 2 2
9%BE + 96x %BE - 189%BE - 2520x %BE + (- 3840x + 972)%BE + 12096x
arg2= x
2>exit ElementaryIntegration.lfintegrate,58 :
5 4 2 3 2
- 1455%BE - 15520x %BE + (49664x + 27936)%BE + 407400x %BE
+
2 3 2
(769792x - 164997)%BE - 1390592x + 3454976x - 1882188x - 182196
/
2
1787904x - 94284
1>exit ElementaryIntegration.lfintegrate,58 :
+---------+
3+-+2 |3+-+2
(873\|x + 873)\|\|x + 1 + 1687
-----------------------------------
873

+---------+
3+-+2 |3+-+2
(5) (\|x + 1)\|\|x + 1
Type: Union(Expression Integer,...)

So it seems that because the integrand is algebraic over Q(x), Axiom tries to
find a primitive element of the algebraic extension which contains the
integrand. But a primitive elment is not unique, and there might be some
randomness in selecting this primitive element (which involving building a tower
of algebraic extensions to successively remove the radicals appearing in the
integrand, and so there could be some choice of order that is not deterministic
or some factoring algorithm involved in getting to the minimal polynomials at
each stage, with a different integral also at each stage). Tracing
primitiveElement:

)trace FSPRIMELT )math

in addition to

)trace INTEF )math

in repeated calling

complexIntegrate(q,x) did show this to be observed (output omitted). I am
running out of time right now to trace the primitiveElement algorithm, but I
believe that is the reason (still could be wrong of course). Here are at least
two such primitive elements, each with different result. It seems the ground
field used here is Q (rational) or maybe C (complex field).

1<enter FunctionSpacePrimitiveElement.primitiveElement,93 :
3+-+
arg1= \|x
+---------+
|3+-+2
arg2= \|\|x + 1
1>exit FunctionSpacePrimitiveElement.primitiveElement,93 :
+---------+
|3+-+2 3+-+
[primelt= 3\|\|x + 1 + \|x ,

pol1 = ... ,

pol2 = ... ,

6 4 3 2 2
prim= ? - 27? - 56x ? + 243? + 432x ? - 512x - 729]

and

1<enter FunctionSpacePrimitiveElement.primitiveElement,93 :
3+-+
arg1= \|x
+---------+
|3+-+2
arg2= \|\|x + 1
1>exit FunctionSpacePrimitiveElement.primitiveElement,93 :
+---------+
|3+-+2 3+-+
[primelt= - 4\|\|x + 1 + \|x ,

pol1 = ... ,

pol2 = ... ,
6 4 3 2 2
prim= ? - 48? - 98x ? + 768? + 1440x ? - 3375x - 4096]


The other branch of integrate, using lfintegrate, also leads to computing
primitive elements (but less complicated, and apparently CONSISTENT, because the
ground field is Q(x)):

(26) -> lfintegrate(q,x::Symbol)$ElementaryIntegration(INT, EXPR INT)
1<enter ElementaryIntegration.lfintegrate,58 :
+---------+
|3+-+2
|\|x + 1
arg1= |---------
| 3+-+2
\| \|x
arg2= x
1<enter FunctionSpacePrimitiveElement.primitiveElement,93 :
3+-+
arg1= \|x
+---------+
|3+-+2
|\|x + 1
arg2= |---------
| 3+-+2
\| \|x
1>exit FunctionSpacePrimitiveElement.primitiveElement,93 :
+---------+
|3+-+2 2
|\|x + 1 2 6 4 2 - x - 1
[primelt= |--------- ,pol1= x ? - x,pol2= ?,prim= ? - 3? + 3? + --------]

| 3+-+2 2
\| \|x x
2<enter ElementaryIntegration.lfintegrate,58 :
arg1= %CJ
arg2= x
2>exit ElementaryIntegration.lfintegrate,58 :
3
x %CJ
1>exit ElementaryIntegration.lfintegrate,58 :
+---------+
|3+-+2
3+-+2 |\|x + 1
(x \|x + x) |---------
| 3+-+2
\| \|x
-------------------------
3+-+2
\|x

+---------+
|3+-+2
3+-+2 |\|x + 1
(x \|x + x) |---------
| 3+-+2
\| \|x
(26) -------------------------
3+-+2
\|x
Type: IntegrationResult Expression Integer

Hope this give you sufficient leads to go on.


William

[Igor Khavkine]

Post by William Sit
So it seems that because the integrand is algebraic over Q(x), Axiom
tries to find a primitive element of the algebraic extension which
contains the integrand. But a primitive elment is not unique, and
there might be some randomness in selecting this primitive element
(which involving building a tower of algebraic extensions to
successively remove the radicals appearing in the integrand, and so
there could be some choice of order that is not deterministic or some
factoring algorithm involved in getting to the minimal polynomials at
each stage, with a different integral also at each stage). Tracing
)trace FSPRMELT )math

If you also add

)trace PRIMELT )math

you'll see that that William is basically correct. The primitiveElement
routing from PRIMELT gets called with a 4-argument signature and returns
different answers with the same input. Looking at the source
(algebra/primelt.spad.pamphlet) there is in deed a call to
random()$Integer in that routine. If you run the following command a few
times, it will return different results.

primitiveElement(x^(1/3),sqrt(1+x^(2/3)))

BTW, what exactly does this routine do? The documentation is rather
terse and the only mention of "primitive element" in my copy of Geddes
et al. doesn't seem relevant.

So, I guess, the answer to my original question is Yes, Axiom's
integration routine for algebraic functions is non-deterministic. That's
not really a problem, since it seems to return the correct answer but in
different forms. But I am surprized that there doesn't seem to be a
canonical form for algebraic expressions. It may be that I was expecting
too much. However, looking at the last chapter of Geddes et al., there
is a discussion of a canonical form for algebraic functions as they are
fed into the Risch integrator. It doesn't look like any of the
"normalize" routines that I found in Axiom actually perform this trick.
Is this step actually taken or bypassed in the integration code?

Igor

[Page, Bill]

Post by Igor Khavkine
...
Start up Axiom, and type in the following integral
integrate(sqrt(1+x^(-2/3)),x)
Then type in the same integral again
...
Clearly, Axiom takes two different paths through the integration
algorithm, even when give identical input. What is the cause of
the branch? Is there a non-deterministic step somewhere in the
algorithm?

I don't know the cause of this behaviour but I get the same effect
on the Windows version of Axiom.

Post by Igor Khavkine
Incidentally, is there a canonical form for radical expressions
in which the two forms of the answer can be compared and directly
shown to be the same?

It's not based on a canonical form, but this seems to work:

A1 := integrate(sqrt(1+x^(-2/3)),x)
A2 := integrate(sqrt(1+x^(-2/3)),x)
A1-A2

Surprised? :)

Regards,
Bill Page.