open-axiom-test 1.4.1+svn~2626-2ubuntu2 / usr / lib / open-axiom / input / fixed.input

--Copyright The Numerical Algorithms Group Limited 1996.
--fixed.input

-------------------------------------------------------------------
-- daly/10/06/92
-------------------------------------------------------------------
)clear all
-- something is wrong with the integration of serices, when the
-- coeffiennts are not constants. while the series is mathematically
-- correct, this is not what one expects as the taylor series, which
-- does not exist. 
series(x/(x+log(x)))
-- here integrate is treating log(x) as a constant, which is incorrect.
integrate(%)

-------------------------------------------------------------------
-- bmt/10/19/92 integrate problem -- serious, wrong answer
-------------------------------------------------------------------
)clear all
2*sin(t)*sqrt(1+cos(t))
--correct answer is 2*sin(t)/sqrt(1+cos(t))
integrate(%,t)

-------------------------------------------------------------------
-- bmt/10/20/92 multiple complex in type tower
-------------------------------------------------------------------
)clear all
n:Complex ?
-- an invalid type. can't have 2 complex constructors in a tower
n:=x/y+%i

-------------------------------------------------------------------------
-- grabm/09/28/92 coercion bug
-- comment: this coercion must carefully check whether there is a symbol
--          among the kernels which is equal to an indeterminate of the
--          polynomial ring, and not simply considering every expression
--          as a constant
-------------------------------------------------------------------------
)clear all
f:=(a-b-c-d)**2::EXPR INT
f::DMP(['a,'b],EXPR INT)
degree %

-------------------------------------------------------------------
-- bmt/10/26/92 wrong answer
-- comment:I believe this problem simplifies to
-- lfintegrate(sqrt(u**3+u**2),u) which returns the
-- wrong answer due to some confusion in prootintegrate in INTPAF.
-- I think the confusion happens with the use of radPoly and rootPoly.
-- The answer is computed with respect to the result returned by rootPoly
-- but the kernel which is substituted for "y" corresponds to radPoly.
-------------------------------------------------------------------
)clear all
-- should be 2*sin(t)/sqrt(1+cos(t))
integrate(sqrt(1+cos(x)),x)

-------------------------------------------------------------------
-- bmt/09/28/92 bug in ITRIGMNP
-- gives:
--    Cannot convert kernel to gaussian function
-------------------------------------------------------------------
)clear all
integrate(exp(x**2),x)

-------------------------------------------------------------------
-- themos/11/06/92 integration bug
--comment:
--I suggest changing line 11 of limitedLogs in INTTR from:
--      rows := max(degree den, 1 + "max"/[degree(u.contrib) for u in l1])
--to:
--      rows := max(degree den,
--                  1 + reduce(max,[degree(u.contrib) for u in l1],0)$List(N))
--Does you agree, thus max/ empty list would return 0.
-------------------------------------------------------------------
)clear all
f:=log(1-(b*x/(a+c*x**2)))/x
--               2
--            c x  - b x + a
--        log(--------------)
--                  2
--               c x  + a
--   (1)  -------------------
--                 x
integrate(%,x)
--                   2
--                 %J c - %J b + a
--             log(---------------)
--           x          2
--         ++         %J c + a
--   (2)   |   -------------------- d%J
--        ++            %J
expand f
--                 2               2
--        - log(c x  + a) + log(c x  - b x + a)
--   (3)  -------------------------------------
--                          x
integrate(%,x)
--   >> Error detected within library code:
--   No identity element for reduce of empty list using operation
--   max

-------------------------------------------------------------------
-- bmt/11/17/92 interpreter resolve problem
-------------------------------------------------------------------
-- says cannot find an appropriate *
)clear all
%i/m*sqrt(m) 

-------------------------------------------------------------------------
-- tpd/09/22/92 this generates extra brackets
-------------------------------------------------------------------------
)clear all
foo n == matrix[[i for i in 1..n] for j in 1..n]
foo

-------------------------------------------------------------------------
--- tpd/09/24/92 matrix multiplication appears broken in the product.
-- fails with
--   >> System error:
--   The function |deleteAssoc| is undefined.
--   You are being returned to the top level of the interpreter.
-- comment: (barry)
-- This bug is due to the new autoloading scheme for the compiler, the
-- function deleteAssoc needs to be moved into the in core system.
-- As a work around, 
-- if you compile a non existent file, the apropriate files will be loaded,
-- e.g.:
-- )co foo
--  comment: this works in version 2
-------------------------------------------------------------------------
msq := Matrix SquareMatrix(2,POLY INT)
m : msq := zero(2,2)
m(1,1) := matrix([[1,2],[a,b]])
m(1,2) := matrix([[a,b],[2,b]])
m(2,2) := matrix([[1,2],[2,b]])
m
m*m
m**2
m**3
(m*m)*m
mm:=m*m
mm*m

-------------------------------------------------------------------------
-- tpd/09/25/92 new equation.spad from johannes grabmeier
-------------------------------------------------------------------------


--Furthermore, we have function to put 0 or 1 on one side
--and factor the left hand side, after the right hand side
--is 0 and we have an IntegralDomain.
--Johannes

eq1 := (-6*x**3+13*x**2+4)=(-x**4+12*x)
--
--
--            3      2         4
--   (1)  - 6x  + 13x  + 4= - x  + 12x
--                                            Type: Equation Polynomial Integer
--                                       Time: 1.61 (IN) + 1.47 (OT) = 3.08 sec
eq2 := x**4+13*x**2-12*x = 6*x**3-4
--
--
--         4      2          3
--   (2)  x  + 13x  - 12x= 6x  - 4
--                                            Type: Equation Polynomial Integer
--                                       Time: 0.16 (IN) + 0.11 (OT) = 0.27 sec
eq := eq1*y**2+eq2
--
--
--             3      2      2    4      2            4        2     3
--   (3)  (- 6x  + 13x  + 4)y  + x  + 13x  - 12x= (- x  + 12x)y  + 6x  - 4
--                                            Type: Equation Polynomial Integer
--               Time: 0.26 (IN) + 0.01 (EV) + 0.01 (OT) + 1.54 (GC) = 1.82 sec
swap %
--
--
--            4        2     3           3      2      2    4      2
--   (4)  (- x  + 12x)y  + 6x  - 4= (- 6x  + 13x  + 4)y  + x  + 13x  - 12x
--                                            Type: Equation Polynomial Integer
--                                                   Time: 0.07 (OT) = 0.07 sec
% + 4
--
--
--            4        2     3       3      2      2    4      2
--   (5)  (- x  + 12x)y  + 6x = (- 6x  + 13x  + 4)y  + x  + 13x  - 12x + 4
--                                            Type: Equation Polynomial Integer
--                                       Time: 0.69 (IN) + 0.01 (OT) = 0.70 sec
%-6*x**3
--
--
--            4        2       3      2      2    4     3      2
--   (6)  (- x  + 12x)y = (- 6x  + 13x  + 4)y  + x  - 6x  + 13x  - 12x + 4
--                                            Type: Equation Polynomial Integer
--                                       Time: 0.19 (IN) + 0.01 (OT) = 0.20 sec
leftZero %
--
--
--             4     3      2            2    4     3      2
--   (7)  0= (x  - 6x  + 13x  - 12x + 4)y  + x  - 6x  + 13x  - 12x + 4
--                                            Type: Equation Polynomial Integer
--                                       Time: 0.01 (IN) + 0.05 (OT) = 0.06 sec
swap %
--
--
--          4     3      2            2    4     3      2
--   (8)  (x  - 6x  + 13x  - 12x + 4)y  + x  - 6x  + 13x  - 12x + 4= 0
--                                            Type: Equation Polynomial Integer
--                                                   Time: 0.01 (IN) = 0.01 sec
factor lhs %
--
--
--               2       2  2
--   (9)  (x - 2) (x - 1) (y  + 1)
--                                            Type: Factored Polynomial Integer
--                           Time: 0.50 (IN) + 0.09 (EV) + 0.30 (OT) = 0.89 sec
factorAndSplit eq
--
--
--                             2
--   (10)  [x - 2= 0,x - 1= 0,y  + 1= 0]
--                                       Type: List Equation Polynomial Integer
--                                       Time: 0.09 (EV) + 0.21 (OT) = 0.30 sec
inv (eq :: EQ FRAC POLY INT)
--
--
--                            1                                1
--   (11)  - ------------------------------------= - ----------------------
--              3      2      2    4      2            4        2     3
--           (6x  - 13x  - 4)y  - x  - 13x  + 12x    (x  - 12x)y  - 6x  + 4
--                                   Type: Equation Fraction Polynomial Integer
--                                                   Time: 0.03 (IN) = 0.03 sec
- %
--
--
--                          1                              1
--   (12)  ------------------------------------= ----------------------
--            3      2      2    4      2          4        2     3
--         (6x  - 13x  - 4)y  - x  - 13x  + 12x  (x  - 12x)y  - 6x  + 4
--                                   Type: Equation Fraction Polynomial Integer
--                                       Time: 0.01 (IN) + 0.01 (OT) = 0.02 sec


-------------------------------------------------------------------
-- bmt/09/29/92 coercion bug
-- fails with:
-- failed cannot be coerced to mode 
-- (Record (: coef1 (Integer)) (: coef2 (Integer)))
-- fixed on 09/29/92
-------------------------------------------------------------------
)clear all
(p,q):UP(x,INT)
p:=3*x**4+11*x**2-4
q:=9*x**4+9*x**2-4
myNextPrime: (INT,NNI) -> INT
myNextPrime(x,n)==nextPrime(x)$PRIMES(INT)
--  runs forever due to algebra bug in handling leading coefficients
modularGcd([p,q])$InnerModularGcd(INT,UP(x,INT),67108859,myNextPrime)

-------------------------------------------------------------------
-- dewar/10/02/92
-- actually, this was never wrong
-------------------------------------------------------------------
numeric(%e ** %pi)
--gives
   (61)
--  23.1406926327 7926900563 5493084578 4861294789 0532254433 2912200665 
-- 6744971913 2534837468 7163222983 2730035

-------------------------------------------------------------------
-- themos/10/07/92
-------------------------------------------------------------------
)clear all
y:=operator 'y
deqx:=differentiate(y(x),x,2)+differentiate(y(x),x)+y(x)
solve(deqx,y,x)
solve(deqx,y,x=0,[1])
deqt:=differentiate(y(t),t,2)+differentiate(y(t),t)+y(t)
solve(deqt,y,t)
solve(deqt,y,t=0,[1]) -- bug
deqz:=differentiate(y(z),z,2)+differentiate(y(z),z)+y(z)
solve(deqz,y,z)
solve(deqz,y,z=0,[1]) --bug

-------------------------------------------------------------------
-- themos/10/07/92 equation
-------------------------------------------------------------------
)clear all
y:=operator 'y
deq:=D(y(x),x)+x**2=(y x)/x-(y x)**2


-------------------------------------------------------------------
-- bmt/10/08/92 laplace
-------------------------------------------------------------------
)clear all
laplace(exp(-x**3)*x**7,x,s) 

-------------------------------------------------------------------
-- bmt/09/30/92 nonlinear ODE 
-------------------------------------------------------------------
-- Any hope of solving the following which appeared in Barry
-- Simon's computer algebra test suite:
-- (Apparently it has Bernoulli form)
y:=operator 'y
x**2 * D(y x, x) + 2*x*(y x) - (y x)**3 = 0
solve(%,y,x)


-------------------------------------------------------------------
-- sit/10/08/92 factor bug
-------------------------------------------------------------------
p:POLY FRAC INT := 3*(x+1)
-- is wrong (missing factor of 3).
factor(p)**2 

-------------------------------------------------------------------
-- henderson/10/08/92 
-------------------------------------------------------------------
)clear all
--TextFile cliams that the operator readLineIfCan! returns "failed" when
--the end of the file is reached.  In fact, I get a system error, which is
--not very useful:
--(16) ->readLineIfCan!(fp)
--   >> System error:
--   Unexpected end of #<input stream "/home/mcd/work/work">.
--   You are being returned to the top level of the interpreter.
fout:TextFile:=open("/tmp/foo","output")
write!(fout,"foo")
close!(fout)
fin:TextFile:=open("/tmp/foo","input")
readLineIfCan!(fin)
readLineIfCan!(fin)
close!(fin)
)lis (system "rm /tmp/foo")

-------------------------------------------------------------------
-- bmt/10/08/92 factoring over SAEs
-------------------------------------------------------------------
)clear all
a | a**2+1
(x+a)*(x+a+1)
factor %

-------------------------------------------------------------------
-- miller/10/09/92 local vars in types
-------------------------------------------------------------------
)clear all
--Fixed (Victor's) bug with using local variables in type expressions
--assigned to local variables.  It now gives a message and correctly uses
--interpret-only mode.
--The function changed was getAndEvalConstructorArgument from i-spec1.boot.

-------------------------------------------------------------------
-- themos/10/19/92
-------------------------------------------------------------------
)clear all
-- Do this in a virgin system
)set expose add constructor SquareMatrix
S2:= SquareMatrix(2,FRAC POLY INT)
V2: S2 := matrix([[v,-v],[-v,v]])
I2: S2 := 1
m:=5
l: List(S2) := append(cons(V2+h*I2,_
            [(V2+2*h*I2) for i in 2 .. (m-1)]),_
             [V2+h*I2])
A: SquareMatrix(m, S2) := diagonalMatrix(l)
-- did you get an error of deleteAssoc not defined ??

-- load definition of deleteAssoc
)li (load "/spad/obj/rios/interp/c-util")

A: SquareMatrix(m, S2) := diagonalMatrix(l)
--now , can print.

-------------------------------------------------------------------
-- bmt/10/19/92 squareFree bug
-------------------------------------------------------------------
-- has an extra factor of 2.
squareFree((2*x*y+1)*(x*y+1)**2)

-------------------------------------------------------------------
-- bmt/10/19/92 division by zero error
-------------------------------------------------------------------
limit(atan(1/sin(x)),x=0)

-------------------------------------------------------------------
-- bmt/10/19/92 limit problem
-------------------------------------------------------------------
)clear all
-- the left and right hand limits don't agree
limit(atan(-sin(x)/(cos(x)+e)),x=acos(-e))
-- the left and right hand limits don't agree
limit(atan(1/(cos(x)+e)),x=acos(-e))

-------------------------------------------------------------------
-- grabmeier/10/19/92 unable to compute the determinant of an 8x8
-------------------------------------------------------------------
)clear all
D := MATRIX FRAC(POLY INT)


d : (INT, Boolean) -> POLY INT
d(i,ss)  ==
  ex := i rem 4
  if ex < 0 then ex := ex+4
  ss =>
    ex = 0 => 's
    ex = 1 => 'sd
    ex = 2 => 'sdd
    'sddd
  ex = 0 => 1
  ex = 1 => 'd
  ex = 2 => 'dd
  'ddd



-- 1,d,dd,ddd,s,sd,sdd,sddd

mTV4 : D := new(8,8,0)

for i in 1..8 repeat
  for j in 1..8 repeat
    mTV4(i,j) :=
      i <= 4 =>
        j <= 4 => d(i+j-2, false)
        d(-i+j,true)
      j <= 4 => d(i+j-2,true)
      d(-i+j,false)
mTV4
gdd4 := determinant mTV4

-------------------------------------------------------------------
-- sutor/09/28/92 operator bug?
-------------------------------------------------------------------
)clear all
L n ==
  n = 0 => 1
  n = 1 => x
  (2*n-1)/n * x * L(n-1) - (n-1)/n * L(n-2)
dx:=operator("D")::OP(POLY FRAC INT)
evaluate(dx,p+-> differentiate(p,'x))
E n == (1-x**2)*dx**2-2*x*dx+n*(n+1)
L 15
-- used to fail on this line
E 15

-------------------------------------------------------------------
-- bmt/10/12/92 EFSTRUC recursion problem
-------------------------------------------------------------------
)clear all
bug:=(1+x**(1/4))**(1/3)/(x**(1/2))
-- gives a value stack overflow
integrate(bug,x) -- seems to go into infinite recursion in rischNormalize

-------------------------------------------------------------------
-- james/10/28/92 coerce bug
-------------------------------------------------------------------
)clear all
-- improper character in Roman numeral: G
g::ROMAN

-------------------------------------------------------------------
-- grabm/10/28/92 runs forever
-------------------------------------------------------------------
)clear all
-- this is quick
factor 1068303355883998767544567663620885466990173600
-- but this runs forever
sqrt 1068303355883998767544567663620885466990173600


-------------------------------------------------------------------
-- themos/11/05/92 fortran output bug
-------------------------------------------------------------------
)clear all
--      REAL T7,T6,T5,T4,T3,T2,T1
--      T1=x*x
--      T2=2.*y**5
--      T3=4.*T1
--      T4=y**3
--      T5=2.*T7
--      T6=x**3
--      T7=x**4
--      R46=((T2+(T3+8.*x+8.)*T4+(T5+8.*T6+(-40.*T1))*y)*SIN(x)+(-T2+(-T3+
--     &16.*x)*T4+(-T5+16.*T6)*y)*COS(x))/(y**8+4.*T1*y**6+6.*T7*y**4+4.*x
--     &**6*y*y+x**8)
-- T7 is referenced before it is defined
a1:=sin(x)/(x**2+y**2);a2:=D(a1,[x,y]);a2+D(a2,x)

-------------------------------------------------------------------
-- grabm/11/05/92 
-- comment:
--I believe the "bug" is that the interpreter is converting gb which is
--a list of dmp style polynomials to a list of POLY type polynomials.
--Unfortunately the result of such a conversion is no longer a groebner basis.
--I don't see any reasonable way to prevent this, and would therefore have to
--classify this as "user error", when working with grobner bases, the user
--needs to force all his polynomials to be the appropriate type (thus having
--the correct ordering).
-------------------------------------------------------------------
-- expressing a symmetric functions in terms of elementary symm. functions
-- working in POLY INT is o.k.
R := FRAC INT
--   (1)  Fraction Integer
lpx : List POLY R := [x1+x2+x3, x1*x2+x1*x3+x2*x3,x1*x2*x3]
--   (2)  [x3 + x2 + x1,(x2 + x1)x3 + x1 x2,x1 x2 x3]
lip := [lpx.1-e3, lpx.2-e2, lpx.3-e1]
--   (3)  [x3 + x2 + x1 - e3,(x2 + x1)x3 + x1 x2 - e2,x1 x2 x3 - e1]
gbp := groebner lip
--   (4)
--                         2                   2
--   [x3 + x2 + x1 - e3, x2  + (x1 - e3)x2 + x1  - e3 x1 + e2,
--      3        2
--    x1  - e3 x1  + e2 x1 - e1]
normalForm(x1**2+x2**2+x3**2,gbp)
--          2
--   (5)  e3  - 2e2
-- working in DMP is o.k.
dmp := DMP([x1,x2,x3,e1,e2,e3],INT)
--   (6)  DistributedMultivariatePolynomial([x1,x2,x3,e1,e2,e3],Integer)
li : List dmp := [lpx.1-e1, lpx.2-e2, lpx.3-e3]
--   (7)  [x1 + x2 + x3 - e1,x1 x2 + x1 x3 + x2 x3 - e2,x1 x2 x3 - e3]
gb := groebner li
--   (8)
--                         2                     2
--   [x1 + x2 + x3 - e1, x2  + x2 x3 - x2 e1 + x3  - x3 e1 + e2,
--      3     2
--    x3  - x3 e1 + x3 e2 - e3]
p:dmp:=(x1**2+x2**2+x3**2)
--          2     2     2
--   (9)  x1  + x2  + x3
normalForm(p,gb)
--           2
--   (10)  e1  - 2e2
-- but we have some problems (coercion? ordering?) here:
normalForm(x1**2+x2**2+x3**2,gb)
--            2                      2              2
--   (11)  2x2  + (2x1 - 2e1)x2 + 2x1  - 2e1 x1 + e1




-------------------------------------------------------------------
-- bmt/11/17/92 Ifintegrate
-------------------------------------------------------------------
--It seems that one cannot simply parametrize reducible curves.
--Any attempt at parametrization can only follow one branch.
--For curves of the form y**2=a**2*x**2+b*x+c, these are only reducible
--if a**2*x**2+b*x+c is a perfect square. Basically I think we need to
--choose 1 branch in this case. Under some circumstances we could
--try to proceed down both branches and then unify the result, but for
--the time being, perhaps we should just perform the substitution
--y=sqrt(a**2*x**2+b*x+c) whenever that is a perfect square.
--(then we don't even need to do any back subsitution or change to dx).
--
--For another example, perhaps due to the same essential difficulty try
integrate(normalize(sqrt(1+cos(x)),x),x)
-- this gets a non-invertible error

--I have fixed integrate(sqrt(1+cos(x)),x) and
--integrate(normalize(sqrt(1+cos(x)),x),x)
--fix was to function prootintegrate in INTPAF in intaf.spad


-------------------------------------------------------------------
-- james/02/11/93 integrate
-- used to give error:
--   No identity element for reduce of empty list using operation
--   max
-------------------------------------------------------------------
integrate(((-x-1)*log((x**2+x))**2+2*log(x))/(x+1),x)

-------------------------------------------------------------------
-- dewar/02/15/93
-------------------------------------------------------------------
)clear all
--i::Polynomial(Integer)
--%::Union(Symbol,Polynomial Integer)
--list %
--
--this is a bug in old style union handling, this style of unions
--will become obsolete, so the suggested fix is to use tagged unions instead:
i::POLY INT
%::Union(s:Symbol, p:POLY INT)
list % -- works

-------------------------------------------------------------------
-- grabm/03/04/93
-- parser bug? not for system commands
-------------------------------------------------------------------
)clear all
-- gives error: Sorry, but no option begins with upd;
)load /spad/mnt/rios/algebra/EQ )upd; a:=1
-- works
)load /spad/mnt/rios/algebra/EQ )upd ; a:=1

-------------------------------------------------------------------
-- bronstein/03/08/93
-------------------------------------------------------------------
)clear all
I:=operator 'I
J:=operator 'J
eq := mu * D(I x,x) = - (K + S) * I(x) + S*J(x)
solve(eq,I,x)
--> wrong particular solution

-------------------------------------------------------------------
-- dewar/03/16/93
-- bug writing backslash to files
--comment:
--This is not a bug, what is written to the file is the lisp representation
--of the string "\\test" (since '\' is lisp escape character).
--The file is read and written using lisp primitives and thus the escape
--characters need to be doubled. This is consistent, since if you read the
--item back in from the file, the correct string is reconstituted.
--
--If you want to create a text file which simply contains lines of text
--instead of lines of lisp strings, use:
--ofile : TextFile := open("test","output")  ...
--then you will get a file of literal strings with no string delimiters and
--no escape characters.
-------------------------------------------------------------------
)clear all
ofile: File String := open("/tmp/test","output")
-- this writes "\\\\test" but should write "\\test"
write!(ofile,"\\test"::String)
close! ofile

-------------------------------------------------------------------
-- themos/04/07/93
-- this sometimes fails
-- comment:
--thanks for tracing down the problem. Basically the lines for monicizing
--the result need to be done before we try to 
--retract back from the algebraic extension.
-------------------------------------------------------------------
)clear all
pol:DMP([x,y,z],PF(2)):=x**2*y**2+x**2*y*z+x**2*z**2+x*y*z**2+y**3*z+y*z**3
factor pol

-------------------------------------------------------------------
-- williamson/04/21/93
-- comment:
--I ran into some problems when factoring polynomials over number fields.
--Evidently the interpretor doesn't choose the correct function 'factor'.
--It uses a function from MULTFACT, when it should be using a function from
--SAEFACT. 
-------------------------------------------------------------------
)clear all
up := UP('w,FRAC INT)
p : up := w**4 + w**3 + w**2 + w + 1
sae := SAE(FRAC INT,up,p)
q : UP('x,sae) := x**5 - 1
factor q -- used to report: x**5-1 
saefact := SAEFACT(up,sae,UP('x,sae))
factor(q)$saefact

-------------------------------------------------------------------
-- grabm/04/30/93
-- comment:
--This bug has been fixed at yorktown, in radix.spad, change line 163, the
--third line of intgroup to be:
--      empty? rest li => intToExpr first(li)
--and then recompile RADIX
-------------------------------------------------------------------
)clear all
-- results in 10 and not A
10::RadixExpansion(16)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch (manuel bronstein)
-------------------------------------------------------------------
)clear all
r:=rule 'x == 1
-- used to return x, now returns 1
r x

-------------------------------------------------------------------
-- ADK@scri.fsu.edu (tony kennedy)
-------------------------------------------------------------------
)clear all
factor(-12)
%**2

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch (manuel bronstein)
-- >> Error: cannot retract nonconstant polynomial
-------------------------------------------------------------------
)clear all
complexNumeric(log(sqrt(-3)))

-------------------------------------------------------------------
-- quitte@knuth.univ-poitiers.fr/8/15/93 (Claude Quitte)
-------------------------------------------------------------------
)clear all
--)abbreviation package TEST Test
--
--Test() : with
--
--   leftProductBy : Integer -> (Integer -> Integer)
--   rightProductBy : Integer -> Mapping(Integer, Integer)
--
-- == add
--
--   leftProductBy(n) == n * #1
--
--   rightProductBy(n : Integer) : Mapping(Integer, Integer) == #1 * n
--------------------------------------------------------------------------
--
--Why is it impossible to specify the target/cible types for the
--functions leftProductBy and rightProductBy in the body of the package ???
--
--The code for leftProductBy is correct but the code for rightProductBy
--troubles the compiler !!!! It seems to compile a LOCAL function with
--ANOTHER signature.
--
--------------------------------------------------------------------------
--   initializing NRLIB TEST for Test
--   compiling into NRLIB TEST
--   compiling exported leftProductBy : Integer -> Integer -> Integer
--Time: 0.35 SEC.
--
--   compiling local rightProductBy : Integer -> Integer
--****** comp fails at level 2 with expression: ******
--error in function rightProductBy
--
--****** level 2  ******
--$x:= #1
--$m:= $EmptyMode
--$f:=
--((((|n| # #) (|rightProductBy| # # #) (|n| #) (|rightProductBy| # # #) ...)))
--
--   >> Apparent user error:
--   no mode found for
--   #1
--   You are being returned to the top level of the interpreter.
------------------------------------------------------------------------
-- answer (barry):
--
--Your problem in TEST about specifying target types which are Mapping's
--was a bug in the compiler and has been fixed for OpenAxiom release 2.0.
--

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/10/6/93 (Manuel Bronstein)
-------------------------------------------------------------------
)clear all
--Here is algfunc.spad with an old bug fixed: eventh-roots of negative floats
--coerced to EXPR FLOAT (e.g. sqrt(-1::EXPR FLOAT)) used to produced an error
--message since the sqrt from FLOAT was called. 
--I have now fixed iroot from AF(R,)
--to call the sqrt from R only if either R is algebraically closed, or the
--root is odd, or the argument is positive. Here's the new behaviour:
sqrt(-1::EXPR FLOAT)
sqrt(2::EXPR FLOAT)
nthRoot(-2::EXPR FLOAT, 3)
nthRoot(-2::EXPR FLOAT, 4)
--As a side-effect, this fixes the problem with numeric
-- (which was instantiating sqrt(-1)$EXPR(FLOAT)).

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/10/6/93 (Manuel Bronstein)
-------------------------------------------------------------------
)clear all
--Here is an updated efstruc.spad where ker2trigs now knows about abs.
real abs(4 + %i * 5) -- does the right thing.

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/10/5/93 (Manuel Bronstein)
-------------------------------------------------------------------
)clear all
--Here is again an enhanced elemntry.spad. What happens now is that
exp(q * %i * %pi)
--gets automatically simplifed whenever q is a rational number whose
--denominator is one of  {1,2,3,4,6}.
exp(5/3*%i*%pi)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/10/4/93 (Manuel Bronstein)
-- luczak@nag.com (Richard Luczak)
-------------------------------------------------------------------
)clear all
exp(log(-1))
sum((-1)**k * (k+m),k=0..n)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/10/4/93 (Manuel Bronstein)
-------------------------------------------------------------------
)clear all
--simplifies (abs now checks quotients instead of just retraction to R).
abs((1/2)::EXPR(INT))

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/10/4/93 (Manuel Bronstein)
-------------------------------------------------------------------
)clear all
integrate(1/(x**2 + %i*a),x)

-------------------------------------------------------------------
-- bmt@spadserv.watson.ibm.com/9/28/93 (Barry Trager)
-------------------------------------------------------------------
)clear all
limit(1/2**n,n=%plusInfinity)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/9/22/93 (Manuel Bronstein)
-------------------------------------------------------------------
)clear all
--Here is an update to efstruc.spad that handles complex constants much better.
--Negative square roots and negative logs are now recognized as complex and
--treated properly by real?, real, imag, and complexForm:
x := sqrt(-3) + sqrt 2 + sqrt(- exp a) + log(-a**2-1)
real? x
real x
imag x
--As a result, integrals involving sqrt(-2) etc... are now treated correctly
--(this was the case only for sqrt(-1) with the older version).

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/9/22/93 (Manuel Bronstein)
-------------------------------------------------------------------
)clear all
-->> haha := rule x*x == z
-->> haha(a*a + b*b + c**2)        -->  3z
-->> haha(a*a + b*b + c**2 + d*d)  -->  z
-->>
-->> the bug is that the last line returns z instead of 4z.
--
--Sorry guys, this is not a bug: haha is so general a rule that it matches
--the integer 4 (as 2 squared), so the rewrite chain for the last example is:
--
--a*a + b*b + c**2 + d*d  --->  z + z + z + z  = 4 * z ---> z * z ---> z
--
--Here is a console showing what exactly happens:
haha := rule x*x == z
haha 4
haha 3
haha(4*z)
haha(3*z)
--To see the whole rewrite chain:
a*a + b*b + c**2 + d*d
applyRules([haha], %, 1)$APPRULE(INT,INT,EXPR INT)
applyRules([haha], %, 1)$APPRULE(INT,INT,EXPR INT)
applyRules([haha], %, 1)$APPRULE(INT,INT,EXPR INT)
--I think it's actually the proper behavior for the pattern matcher. This
--example shows that rules can bite, something Mma users are quite aware about!

-------------------------------------------------------------------
-- dewar/9/16/93 (mike dewar)
-------------------------------------------------------------------
--There's an odd discontinuity about the behaviour of the following
--function:

harm(1) == 1
harm(n) == harm(n-1) + 1/n
harm : Integer -> Fraction Integer
harm(1023) -- takes a little while, as expected
harm(1024) -- takes forever?

-------------------------------------------------------------------
-- williams@inf.ethz.ch/9/1/93 (Clifton Williamson)
-------------------------------------------------------------------
--> we should probably support integrate(%,x) on series types.
--> x could be either the variable of expansion or a parameter variable.
--
--We've discussed this before with regard to the function 'differentiate'.
--The technical problem is that "generalized" power series may have
--"coefficients" involving the series variable:
series(x**x,x=0)
--In this case, as currently implemented, 'differentiate' returns incorrect
--answers:
)set mess test off
differentiate %
)set mess test on
--The solution we discussed (perhaps "resolution" is a better word) was to
--create a separate series type for generalized series.  This would be a
--carbon copy of UPXS, except that differentiate(series,variable(series)) and
--integrate(series,variable(series)) would return an error message.  The
--error message could also suggest that the user first apply 'approximate',
--then compute a derivative or integral.  We would also have a coercion from
--UPXS to this type.  The function 'series' would return a UPXS, when the
--coefficients do not involve the series variable, and a "generalized series",
--when the coefficients involve the series variable.
--
--If this sounds cool to you, I'll go ahead with it when I have time.  (I'm
--trying to meet a 16 September deadline for a MEGA '94 submission.  I'm still
--getting my results together, so time is tight.  Let me know if there is a
--time deadline for the AXIOM code!)
--
--> Also we need better tools for taking finite truncations of series objects.
--> In particular since we now mostly use series(expr,var) to create a series,
--> UPXS if over a coefficient domain which supports the symbol of expansion,
--> should be able to truncate a series to an EXPR.
--
--I thought I had you on this one!  The signature is there (in pscat.spad):
--
--    if Coef has coerce: Symbol -> Coef then
--      if Coef has "**":(Coef,Expon) -> Coef then
--        approximate: ($,Expon) -> Coef
--          ++ \spad{approximate(f)} returns a truncated power series with the
--          ++ series variable viewed as an element of the coefficient domain.
--
--It certainly works for Laurent series:
laurent(cos(a+x)/x,x=0)
approximate(%,3)
--But, unfortunately, I never implemented it for Puiseux series:
puiseux(cos(a+x)/x,x=0)
approximate(%,3)
series(cos(x**(2/3) + a),x=0)
approximate(%,2)
approximate(%% 1,7/5)


-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/10/6/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
)set message type off

-- Here is algfunc.spad with an old bug fixed: eventh-roots of negative floats
-- coerced to EXPR FLOAT (e.g. sqrt(-1::EXPR FLOAT)) used to produced an error
-- message since the sqrt from FLOAT was called. 
-- I have now fixed iroot from AF(R,)
-- to call the sqrt from R only if either R is algebraically closed, or the
-- root is odd, or the argument is positive. Here's the new behaviour:
sqrt(-1::EXPR FLOAT)
sqrt(2::EXPR FLOAT)
nthRoot(-2::EXPR FLOAT, 3)
nthRoot(-2::EXPR FLOAT, 4)
--As a side-effect, this fixes the problem with numeric
-- (which was instantiating sqrt(-1)$EXPR(FLOAT)).

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/9/22/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
--Here is a (rather major) bug fix to ffactor in FSUPFACT.NRLIB. It causes
--a large family of integrals to return 0, because ffactor(?**2+expr) returned
--?**2 when expr involved a parameter. This is fixed now.

-------------------------------------------------------------------
-- jhd@maths.bath.ac.uk/8/15/93 James Davenport
-------------------------------------------------------------------
)clear all
integrate(1/(x*(log(x)**2+a**2-1)),x)
--   (1)
--   VCONCAT
--                             log(x)
--                       atan(---------)
--                             +------+
--                             | 2
--           2                \|a  - 1
--      ((- a  + 1<0) -> ---------------)
--                           +------+
--                           | 2
--                          \|a  - 1
--  ,
--      PAREN
--               2
--           (- a  + 1>0)
--         ->
--                                      +--------+
--                          2    2      |   2           2
--               log((log(x)  - a  + 1)\|- a  + 1  + (2a  - 2)log(x))
--             +
--                           2    2
--               - log(log(x)  + a  - 1)
--          /
--               +--------+
--               |   2
--             2\|- a  + 1
--  ,
--                            1
--      ( Or (a - 1=0) -> - ------)
--           (a + 1=0)      log(x)
--                                          Type: ConditionalExpression Integer
--
--is still fairly ugly.

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/8/9/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
--Here is efstruc.spad with a change to normalize so that
normalize(2**(1/2) + 2**(1/4)) now returns  2**(1/4) ** 2 + 2**(1/4).


-------------------------------------------------------------------
-- satoshi@yktvmv.vnet.ibm.com/8/7/93 satoshi hamaguchi
-------------------------------------------------------------------
)clear all
-- In axiom, both
  integrate(%e**x,x=0..1)
--and
  integrate(log(x),x=1..2)
--yield a warning  potentialPole    .
--If I use "noPole", they give right answers. Isn't this a bug?
--... Satoshi

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/8/5/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
simplify(2**(1/3)*2**(1/2)) -- will return 2**(5/6).

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/8/4/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
integrate(1/sqrt(1+cos(x)), x)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/8/4/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
normalize atan(cos(x)/sin(x))  -- used to be division by 0 error

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/8/4/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
a := 2**(1/6)
[a**n for n in 2..13]

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/8/4/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
int:=sqrt(a*(1-u**2)/(1+u**2))/u
integrate(eval(int,a=1),u) -- works
integrate(eval(int,a=sqrt(-1)),u) -- seems to run forever
integrate(eval(int,a=1)*(-1)**(1/4),u) -- dies after a long time
 -- with an elt index error.

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/8/4/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
sqrt((1-x**2)*(1-k**2*x**2))
integrate(x/%,x)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/7/26/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
rad:=last zerosOf((2+y)**8-3,y)
k:=first kernels %
eval(rad,k,rad)

-------------------------------------------------------------------
-- nagttt@vax.ox.ac.uk/7/15/93 themos tsikas
-------------------------------------------------------------------
)clear all
f := (x - y) / (x + y)
eval(f,x=1/x)

-------------------------------------------------------------------
-- nagttt@vax.ox.ac.uk/7/15/93 themos tsikas
-------------------------------------------------------------------
)clear all
digits 200
--I claim that these two should give the same result (they do in Reduce):
--(precision was 200)
a:=4*sin(2*%pi/9)*sin(5*%pi/9)/sqrt(3)
b:=1/(2*sin(%pi/9))
a::EXPR FLOAT
b::EXPR FLOAT                                 

-------------------------------------------------------------------
-- bmt@spadserv.watson.ibm.com/9/15/93 (Barry Trager)
-------------------------------------------------------------------
)clear all
limit(tanh(x),x=%plusInfinity)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/6/14/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
a := x::EXPR INT
b := x::EXPR COMPLEX INT
zeroOf(a**4+1,x)
zeroOf(b**4+1,x)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/6/11/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
normalize(0**a) -- now returns 0  (was crashing before)

-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/6/11/93 (manuel bronstein)
-------------------------------------------------------------------
)clear all
-  new export  complexForm(f)  which returns the rectangular complex form of f:
log(a+%i * b)
complexForm %
complexIntegrate(1/(x-%i),x)
complexForm %

-------------------------------------------------------------------
-- jcf@posso.ibp.fr
-------------------------------------------------------------------
--  normalForm est faux pour les polynomes a coef dans un anneau:
D:=DMP([a,b],INT)
l:List D:=[2*a-1,b-a]
g:=groebner l
p:D:=a
normalForm(p,g)
-- was  1
--   au lieu de 1/2 !!

-------------------------------------------------------------------
-- nagttt@vax.ox.ac.uk/6/9/93 themos tsikas
-------------------------------------------------------------------
)clear all
squareFreePart(50)

-------------------------------------------------------------------
-- mcd@maths.bath.ac.uk/6/4/92 mike dewar
-------------------------------------------------------------------
)clear all
f:=complexIntegrate(1/((x-%i)*(x-2*%i)),x)
g:=subst(f,x=1)
imag g

-------------------------------------------------------------------
-- daly
-------------------------------------------------------------------
)clear all
-- not wrong but the order is obscure...
primes(1,10)

-------------------------------------------------------------------
-- nagttt@vax.ox.ac.uk/1/11/93 themos tsikas
-------------------------------------------------------------------
)clear all
m:=matrix[[1,a,0],[b,0,a],[0,b,c]]
ll:=radicalEigenvectors m
-- ends in error:
-- Internal system problem in function appInfix : fewer than 2 arguments to
-- an infix function
--
-- but this works
ll 1
ll 2
ll 3


-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/11/24/93 Manuel Bronstein
-------------------------------------------------------------------
--)clear all
--Here's a bug that appears with the )math option only.
--
--)tr MATRIX )ops determinant )math
--
--   Parameterized constructors traced:
--      MATRIX
-- this integrate is commented out since it runs forever (at least many hours)
--integrate(sqrt(1/x+1)*sqrt(x)-sqrt(x+1),x)
--1<enter Matrix.determinant,26 :
--
-- Internal system problem in function appInfix :
-- fewer than 2 arguments to an infix function

-------------------------------------------------------------------
-- copper@yktvmv/12/1/93 Don Coppersmith
-------------------------------------------------------------------
)clear all
--Attached is the list of integers n such that 2**512-n is prime
--and n is between 0 and 5000:
--
--   4893,4653,4475,4005,3893,3669,3459,3143,2967,
--   2807,2529,1827,1695,975,875,629,569

--It was gotten from OpenAxiom by issuing the commands
qrimes : Stream Integer := generate(nextPrime,2**512-5000)
rrimes := [ 2**512-p for p in qrimes while p < 2**512 ]
srimes := complete rrimes
[srimes.i for i in [1..18]]
[srimes.i for i in [10..18]]

-------------------------------------------------------------------
-- miker@num-alg-grp.co.uk/12/13/93 (Mike Richardson)
-------------------------------------------------------------------
)clear all
--
--This is fixed now. You get an Expression Float.
--
--> 
--> I suspect that this is an example of the same thing - but it looks
--> worse:
--> 
--> (145) ->X := log(0.7*%i*x)
--> 
-->    >> Error detected within library code:
-->    negative root
--> 
--> Mike
--> 


-------------------------------------------------------------------
--  sutor@yktpub.watson.ibm.com/12/15/93 (Robert S. Sutor)
-------------------------------------------------------------------
)clear all
--enter Integer in input box
--choose domains
--choose Description
--click on Convertible(String)
--
--page cannot be formatted

-------------------------------------------------------------------
-- bmt@posso.dm.unipi.it12/15/93 (barry trager)
-------------------------------------------------------------------
)clear all
log2() --> Float
exp1() --> DoubleFLoat ???


-------------------------------------------------------------------
--  quitte@knuth.univ-poitiers.fr/12/14/93 Claude Quitte
-------------------------------------------------------------------
)clear all
--There is a bug in the function delete_! for the type List ; this function,
--has two arguments, a list L and an integer n : delete_!(L, n) returns
--a list L with the n nth-item deleted, MODIFYING its argument L.
--
--But for n = 1, the argument L is not MODIFIED !!! (exactly for n =
--minIndex(L)). An extract of the package ListAggregate (abbreviation LSAGG,
--file aggcat.spad).
--
--
----------------------------------------------------------------------
--ListAggregate(S:Object): Category == Join(StreamAggregate S,
--   FiniteLinearAggregate S, ExtensibleLinearAggregate S) with
--      list: S -> $
--	++ list(x) returns the list of one element x.
-- add
--   ..........
--
--   delete_!(x:$, i:I) ==
--     i < (m := minIndex x) => error "index out of range"
--     i = m => rest x   			****** <- here is the error
--     y := rest(x, (i - 1 - m)::N)
--     setrest_!(y, rest(y, 2))
--     x
--
--A little session :
--
L : List(String) := ["There is", "it seems", "a real bug", "here"]
--   (1)  ["There is","it seems","a real bug","here"]
--                                                            Type: List String
--
L1 := delete!(L, 1)
--   (2)  ["it seems","a real bug","here"]
--                                                            Type: List String
L
--   (3)  ["There is","it seems","a real bug","here"]
--                                                            Type: List String
--
L2 := delete!(L, 2)
--   (4)  ["There is","a real bug","here"]
--                                                            Type: List String
L
--   (5)  ["There is","a real bug","here"]
--
--Is there a way to fix this bug quickly ?? (I am ready to modify the source
--OpenAxiom aggcat.spad).
--The problem with ! functions lies in the documentation not in the code. We
--need to explain more clearly that ! functions are ALLOWED to update their
--arguments, not REQUIRED to do so. 
--Whenever using ! functions you should always
--use the returned result, not try to use the updated argument. This is how
--the code is designed. 
--For example if you try to do insert! into an empty list,
--the argument is not updatable and thus a new list is created. The reason for
--having ! operations is for efficiency to do less copying, not to guarantee
--that given arguments will always be updated. So the correct use is to only
--depend on the results returned by updating functions.


-------------------------------------------------------------------
--  quitte@knuth.univ-poitiers.fr/1/13/94 Claude Quitte
-------------------------------------------------------------------
)clear all
K := Fraction(Integer)   
PolK := UP('X, K) 
X : PolK := monomial(1, 1) 
n : PositiveInteger := 15 
E := SimpleAlgebraicExtension(K, PolK, X**n + X**(n-3) -1) 
y : E := X::E
minimalPolynomial(y)$E
--
--Internal Error
--The function minimalPolynomial with signature SimpleAlgebraicExtension(
--Fraction Integer,UnivariatePolynomial(X,Fraction Integer),X**15+X**12-1)
---> UnivariatePolynomial(X,Fraction Integer) is missing from domain
--SimpleAlgebraicExtension(Fraction (Integer))
--(UnivariatePolynomial X (Fraction (Integer)))((15 1 . 1) (12 1 . 1) (0 -1 . 1))

-------------------------------------------------------------------
-- luczak@nag.com/02/20/94 (Richard Luczak)
-------------------------------------------------------------------
)clear all
--
--I'm trying to define the following rule
--
tr := rule cos(x)**(n | integer? n and even? n)==(1-sin(x)**2)**(n/2)
--
--OpenAxiom returns
--
--   Cannot find a definition or library operation named even? with argument
--   types
--                             Expression Integer
--
--   We will attempt to step through and interpret the code.
--
--                                  n
--                                  -
--              n             2     2
--   (2)  cos(x)  == (- sin(x)  + 1)
--                        Type: RewriteRule(Integer,Integer,Expression Integer)
--
--
--and henceforth the rule does not work as desired. Variations on the same bug
--occur involving positive?. Of course, I get what I want using integer?(n/2)
--but ... Is this a bug?
--
-- from barry:
--
--Should use integer? n and even? integer n  instead.
--
--I toyed with the idea of implementing:
--    even? x == integer? x and even? integer x
--    odd? x == integer? x and odd? integer x
--in EXPR, but this would mean that   (x +-> even? x or odd? x)  is not always
--true anymore and I'm not sure it's a good idea. If you guys don't think it's
--a problem, it takes 2 minutes to add it to EXPR R when R has RETRACT INT
--(don't use integer? but retractIfCan@Union(Integer, "failed")).

-------------------------------------------------------------------
-- copper@yktvmv.vnet.ibm.com/02/22/94 don coppersmith
-------------------------------------------------------------------
)clear all
--Questions that came up using OpenAxiom on Aixproj:
--
--How does one find a Pade approximation?  I tried  pade(2,2,y)
--pade(y,2,2)  Pade(2,2,y)  PadeApproximation(2,2,y)
--(where y is a series in x) to no avail.
--
--Complaint: The error message
--
--    Could not find library function "pade" with arguments
--         (integer, integer, series).
--
--is misleading.  If it can't find any "pade", then say
--
--    Could not find library function "pade".
--
--If there is such a function, but I've got the arguments wrong,
--why not say
--
--    The library function "pade" takes arguments
--          (series, integer, integer)
--
--or at least make that information available somehow.
--
-- from bob sutor:
--
--Don,
--  Try using )what op pade  to see what is available. This will give
--you all operations with "pade" in their names. You can also search for
--*pade* via the HyperDoc Browse utility. If you don't use Browse, following
--up with, say, )display op pade, will give you a more information, though
--it is badly formatted (I will look at that).
--
--Aixproj is back-level, so I'm not sure exactly what is up their
--relative to the current system.
--
--I agree we could use better messages regarding unknown functions
--and I will look at that.
--  Bob
--
-- from bob sutor:
--
--The pade operation requires a taylor series as its third object while
--"series" returns a UnivariatePuiseuxSeries. If you set
--  y := taylor(1+x)**(1/5)
--then things work. You can also do
--  pade(1,1,y :: UTS(EXPR INT, x, 0) )
--but this is more work if you really have a Taylor series.
--
--This is unintuitive and there should be a "pade" for Puiseux series that
--fails if the series is not appropriate. Barry? I will now look at adding
--a bit more coercion support to the interpreter to allow Puiseux series
--to possibly retract to Laurent series and then possibly to Taylor series.
--
--Don,
--  I added the retractions I mentioned before to the development system.
--When the system is unsuccessful in matching a UnivariatePuiseuxSeries
--to an operation, it will first try to retract it to a
--UnivariateLaurentSeries. If that fails still, it will try to retract to
--a UnivariateTaylorSeries.
--
--  Be aware that this sort of thing happens when you are evaluating
--an expression. If you are defining a function, you should add
--explicit coercions or create objects of the correct type for
--efficiency. You function may still work, but it may not run in
--a compiled form.

-------------------------------------------------------------------
-- barry@num-alg-grp.co.uk/02/24/94 Barry Trager
-------------------------------------------------------------------
)clear all
sqrt(2)*2.0 -- fails
sqrt(2)::EXPR INT * 2.0 -- works

-- tpd: this now fails with sample : () -> Reference Integer is missing
-- from domain: Integer
-- Internal Error


-------------------------------------------------------------------
-- bronstei@inf.ethz.ch/02/28/94 Manuel Bronstein
-------------------------------------------------------------------
)clear all
f := exp(exp(x)*exp(1/exp(x))) / exp exp x
--)tr LIMITPS )math
--)tr DFINTTLS )math
--)tr DFINTTLS_;findLimit
-- the following exits back to the shell after limit returns "failed":
computeInt(first tower f, f, 0, %plusInfinity, false)$DFINTTLS(INT, EXPR INT)

-------------------------------------------------------------------
-- jhd@maths.bath.ac.uk/03/02/94 James Davenport
-------------------------------------------------------------------
)clear all
--I've fixed coercions of List and Vector to DirectProduct.
--You can now do things like
[1,2,3] :: DirectProduct(3, Fraction Integer)
--
--I also fixed coercions from Vector to List, which seems to have
--been removed via a library change.

-------------------------------------------------------------------
-- jhd@maths.bath.ac.uk/03/09/94 James Davenport
-------------------------------------------------------------------
)clear all
-- used to give a lisp error
x**10+1::Polynomial PrimeField 2

-------------------------------------------------------------------
-- themos@num-alg-grp.co.uk/03/11/94 Themos Tsikas
-------------------------------------------------------------------
)clear all
f(x)==x**2
draw(f,-1..1)

-------------------------------------------------------------------
-- themos@num-alg-grp.co.uk/03/16/94 Themos Tsikas
-------------------------------------------------------------------
)clear all
--The problem was that nobody was around to catch the coerceOrCroaker.
--Just comment out the first line of intCodeGenCoerce1 in i-code.boot
--(the one with the explicit THROW). Bothe the riemannSphereDraw and JHD's
--original example will then work. -- bob sutor
)read conformal
riemannSphereDraw(1..2,0..%pi,4,4,"polar")