intpakX v1.0
written by Grimmer, Markus, Department of Mathematics, University of Wuppertal, Germany, http://www.math.uni-wuppertal.de/wrswt
<© 1999-2002 Scientific Computing/Software Engineering Research Group, University of Wuppertal, Germany>
intpakX
by Krämer, Walter, Geulig, Ilse and Grimmer, Markus, Department of Mathematics, University of Wuppertal, Germany, http://www.math.uni-wuppertal.de/wrswt
<© 1999-2002 Scientific Computing/Software Engineering Research Group, University of Wuppertal, Germany>
intpak
by Corless, R. and Connell, A., University of Western Ontario, Canada
<© 1992-1993 R. Corless and A. Connell>
| > | restart; |
###########################################################################
intpakX v1.0 - Interval Arithmetic in Maple
This is a new, integrated version of the Maple package intpakX for interval computations.
It integrates and updates the package
- intpak by R. Corless and A. Connell
and its extension
- intpakX by I. Geulig and W. Kraemer
intpakX v1.0 was written/updated by M. Grimmer.
intpakX v1.0 provides basic data types and operations for interval arithmetic as well as additional
features for further interval computation.
Through this package, the features of the above mentioned packages are available in a
combined version as a Maple 7 module.
A couple of updates have been made to provide a version working with Maple 7 properly.
For information how to create a Maple package from this document, see the end of this worksheet.
intpakX v1.0 (c) Scientific Computing/Software Engineering Research Group,
Department of Mathematics, University of Wuppertal, Germany
Contact:
kraemer@math.uni-wuppertal.de
markus.grimmer@math.uni-wuppertal.de
http://www.math.uni-wuppertal.de/wrswt
Disclaimer: see end of module
###########################################################################
###########################################################################
Module header
###########################################################################
| > | intpakX:=module() export init,Evalf,min,intpakX_min,max,intpakX_max,ilog10,intpakX_ilog10, Interval_ulp,ulp, Interval_Round_Up,Interval_Round_Down,ru,rd, intpakX_greater,is_in,construct,expinfinity, Interval_exp,infinityln,`&exp`,`&ln`, Interval_ln,addinfinity,Interval_add,`&+`,subtractinfinity, Interval_subtract,`&-`,timesinfinity,Interval_times,`&*`, Interval_reciprocal,inv,Interval_divide,`&/`,sqrtinfinity,Interval_sqrt,`&sqrt`, sqrinfinity,Interval_sqr,`&sqr`, Interval_option_zero,powerinfinity,Interval_Integerpower,`&intpower`, Interval_power,`&**`,Interval_midpoint,midpoint, Interval_width,width,Interval_intersect,`&intersect`, Interval_union,`&union`,Interval_trig_ru,Interval_trig_rd, Interval_scale,Interval_range_values,Interval_sin,Interval_cos,Interval_tan, Interval_arcsin,Interval_arccos,Interval_arctan,Interval_hyp_rd,Interval_hyp_ru, Interval_cosh,Interval_sinh,Interval_tanh,`&sin`,`&cos`,`&tan`,`&arcsin`, `&arccos`,`&arctan`,`&sinh`,`&cosh`,`&tanh`, ext_int_div,rel_diam,mid,compute_naive_interval_range,compute_mean_value_range, compute_monotonic_range,compute_combined_range,compute_taylor_form_range,`&Convex_Hull`, subdivide_equidistant,subdivide_adaptive,interval_list_plot,compute_range, max_abs_error,cni_range3d, subdivide_equi3d,interval_list_plot3d,compute_range3d,newton,compute_all_zeros,newton_plot, newton_with_plot,compute_all_zeros_with_plot, `&cabs`,complex_disc_plot,complex_polynom_plot,`&cadd`,`&csub`,`&cmult`,`&cdiv`,x0_start, `&cmult_opt`,`&cdiv_opt`,horner_eval_cent,horner_eval_opt,centred_form_eval,cexp; option package; |
###########################################################################
Name definitions (see below for body of corresponding procedures)
###########################################################################
| > | min:=intpakX_min; max:=intpakX_max; ilog10:=intpakX_ilog10; width:=Interval_width; ulp:=Interval_ulp; ru:=Interval_Round_Up; rd:=Interval_Round_Down; `&exp`:=Interval_exp; `&ln`:=Interval_ln; `&+`:=Interval_add; `&-`:=Interval_subtract; `&*`:=Interval_times; inv:=Interval_reciprocal; `&/`:=Interval_divide; `&sqrt`:=Interval_sqrt; `&sqr`:=Interval_sqr; `&intpower`:=Interval_Integerpower; `&**`:=Interval_power; midpoint:=Interval_midpoint; `&intersect`:=Interval_intersect; `&union`:=Interval_union; `&sin`:=Interval_sin; `&cos`:=Interval_cos; `&tan`:=Interval_tan; `&arcsin`:=Interval_arcsin; `&arccos`:=Interval_arccos; `&arctan`:=Interval_arctan; `&cosh`:=Interval_cosh; `&sinh`:=Interval_sinh; `&tanh`:=Interval_tanh; |
###########################################################################
Initialization
The init procedure contains type definitions for the various (procedural) interval types.
In Maple, type defintions must be global; they seemingly cannot be module exports.
###########################################################################
| > | init:=proc() |
| > | global Interval_fnlist,`type/num_or_FAIL`,`type/interval_comp`,`type/interval`, `type/complex_disc`,`type/complex_interval`,`convert/interval`,inapply; |
| > |
| > | Interval_fnlist:=[sqrt=`&sqrt`,ln=`&ln`,exp=`&exp`,sin=`&sin`,cos=`&cos`, |
| > | tan=`&tan`,arcsin=`&arcsin`,arccos=`&arccos`,arctan=`&arctan`,sinh=`&sinh`, |
| > | cosh=`&cosh`]; |
- - in proc init - - - - - - - - - - - - - - -
This type is included for ease of writing and clarity of code.
| > | `type/num_or_FAIL`:=proc(a) |
| > | local bool,Constants: |
| > | Constants:={constants}: |
| > | bool:=type(a,numeric) or a=-infinity or a=infinity or a=FAIL |
| > | or member(a,Constants): |
| > | bool: |
| > | end: |
- - in proc init - - - - - - - - - - - - - - -
Simplified version of type interval_comp (cf. intpak/intpakX).
Pi, gamma, Catalan, false, true are not allowed as interval borders.
| > | `type/interval_comp`:=proc(x) |
| > | local bool; |
| > | bool:=type(x,float) or x=FAIL or x=-infinity or x=infinity or x=0; |
| > | bool; |
| > | end: |
- - in proc init - - - - - - - - - - - - - - -
PROCEDURE - TYPE/INTERVAL
This procedure defines the procedural type interval.
It tests an argument to see if it is an interval. An interval
is defined here to be a list with either zero elements,+/- infinity,FAIL
or a list with two floating point members.
`type/interval`:=proc(x)
| > | local bool; |
| > | bool:=false; |
| > | if type(x,list) then |
| > | if nops(x) = 2 then |
| > | if type(x[1],interval_comp) and type(x[2],interval_comp) then |
| > | if intpakX_greater(x[2],x[1]) then bool:=true |
| > | elif x[1] = FAIL or x[2] = FAIL then bool:=true |
| > | else ERROR(`enter the lowest endpoint first`) |
| > | fi |
| > | fi |
| > | elif nops(x)=0 then bool:=true |
| > | fi; |
| > | fi; |
| > | bool; |
| > | end: |
- - in proc init - - - - - - - - - - - - - - -
PROCEDURE - TYPE/COMPLEX_DISC
This procedure implements (i.e. tests an argument to be of) the procedural type complex disc
for use in complex interval arithmetic.
| > | `type/complex_disc`:=proc(z) |
| > | local bool; |
| > | bool:=false; |
| > | if type(z,list) then |
| > | if nops(z)=3 then |
| > | if type(z[3],numeric) then |
| > | if z[3] >= 0 then |
| > | bool:= type(z[1],numeric) and |
| > | type(z[2],numeric); |
| > | fi; |
| > | fi; |
| > | fi; |
| > | fi; |
| > | bool; |
| > | end: |
- - in proc init - - - - - - - - - - - - - - -
PROCEDURE - TYPE/COMPLEX_INTERVAL
This procedure implements (i.e. tests an argument to be of) the procedural type complex interval
(~rectangular interval)
for use in complex interval arithmetic.
| > | `type/complex_interval`:=proc(z) |
| > | local bool; |
| > | bool:=false; |
| > | if type(z,list(interval)) then |
| > | if nops(z)=2 then |
| > | if z[1]=[] then |
| > | ERROR(`First Element of List must be a real interval.`); |
| > | elif z[2]=[] then |
| > | ERROR(`Second Element of List must be a real interval.`); |
| > | else |
| > | bool:=type(z[1][1],numeric) and type(z[1][2],numeric) and |
| > | type(z[2][1],numeric) and type(z[2][2],numeric); |
| > | fi; |
| > | fi; |
| > | fi; |
| > | bool; |
| > | end: |
- - in proc init - - - - - - - - - - - - - - -
PROCEDURE `CONVERT/INTERVAL`
A utility program to convert Maple expressions to interval arithmetic.
convert(1+x + x^2,'interval') returns (1 &+ x) &+ (x &^ 2),
whereas inapply(1+x+x^2,x) yields the operator x -> (1 &+ x) etc.
| > | `convert/interval`:=proc(e) |
| > | local ope, mope, fn; |
| > | # option system: |
| > |
| > | if type(e,interval) or type(e,interval_comp) then e; |
| > | elif type(e,`+`) then |
| > | if type(op(2,e),`*`) and op(1,op(2,e)) = (-1) then |
| > | `convert/interval`(op(1,e)) &- `convert/interval`(op(1,e)-e); |
| > | else |
| > | `convert/interval`(op(1,e)) &+ `convert/interval`(e-op(1,e)); |
| > | fi; |
| > | elif type(e,`*`) then |
| > | if op(1,e)=1.0 then |
| > | `convert/interval`(subs(op(1,e)=1,e)); |
| > | else |
| > | `convert/interval`(op(1,e)) &* |
| > | (`convert/interval`(e/op(1,e))); |
| > | fi; |
| > | ### WARNING: note that `I` is no longer of type `^` |
| > | elif type(e,`^`) then |
| > | if type(op(2,e),posint) then |
| > | `convert/interval`(op(1,e)) &intpower op(2,e); |
| > | elif type(op(2,e),integer) then |
| > | inv(`convert/interval`(op(1,e)) &intpower (-op(2,e))); |
| > | elif op(2,e)=(1/2) then |
| > | &sqrt(`convert/interval`(op(1,e))); |
| > | else (`convert/interval`(op(1,e)) &** `convert/interval`(op(2,e))) |
| > | fi; |
| > | elif type(e,function) then |
| > | ope:=[op(e)]; |
| > | mope:=op(map(`convert/interval`,ope)); |
| > | fn:=subs(Interval_fnlist,op(0,e)); |
| > | fn(mope): |
| > | else e |
| > |
| > | fi: |
| > | end: |
- - in proc init - - - - - - - - - - - - - - -
PROCEDURE INAPPLY
Converts the argument to an interval function.
| > | inapply := proc(); unapply(convert(args[1],'interval'),args[2..nargs]); end: |
- - end proc init global type definitions - - - - - - - - - - - - - - -
| > | with(plots): |
| > | with(geometry): |
| > | macro(evalf=Evalf); |
| > | RETURN(); |
| > | end: |
| > | init(); |
#### end initialization
###########################################################################
| > | Evalf := proc(x) |
| > | local e; |
| > | e := evalf(x); |
| > | if e = -1.*infinity then -infinity elif e = 1.*infinity then infinity else e fi |
| > | end: |
| > |
###########################################################################
BASIC INTERVAL DEFINITIONS
-----------------------------------------------------------------------------
Definition of `type/num_or_FAIL`:
## see init procedure!
-----------------------------------------------------------------------------
Definition of `type/interval_comp`:
## see init procedure!
-----------------------------------------------------------------------------
| > | Interval_ulp := proc(x): |
| > | if x=0 then (0) else |
| > | Float(1,length(op(1,x))+op(2,x)-Digits) |
| > | fi: |
| > | end: |
-----------------------------------------------------------------------------
| > | Interval_Round_Up := proc(x): |
| > | if x=-infinity then x |
| > | elif x=infinity then infinity |
| > | elif x=FAIL then FAIL |
| > | else x + ulp(x) |
| > | fi: |
| > | end: |
-----------------------------------------------------------------------------
| > | Interval_Round_Down := proc(x): |
| > | if x=-infinity then x |
| > | elif x=infinity then infinity |
| > | elif x=FAIL then FAIL |
| > | else x - ulp(x) |
| > | fi: |
| > | end: |
-----------------------------------------------------------------------------
Definition of `type/interval`:
## moved to init procedure!
-----------------------------------------------------------------------------
OPERATOR &GREATER/INTERVAL_GREATER
Checks if (Arg1>=Arg2) applies.
| > | intpakX_greater:=proc(a,b) |
| > | local a_rational, b_rational; |
| > |
| > | if (a=FAIL) or (b=FAIL) then FAIL; |
| > | elif not(type(a,numeric) or member(a,{infinity,-infinity})) then |
| > | ERROR(`first argument must be a numeric, FAIL, infinity or |
| > | -infinity`); |
| > | elif not(type(b,numeric) or member(b,{infinity,-infinity})) then |
| > | ERROR(`second argument must be a numeric, FAIL, infinity or |
| > | -infinity`); |
| > | elif (a=infinity) then true; |
| > | elif (b=-infinity) then true; |
| > | elif (a = -infinity) then false; |
| > | elif (b=infinity) then false; |
| > | elif a > b then true |
| > | elif type(a,rational) and type(b,rational) then |
| > | if a >= b then true else false fi; |
| > | elif type(a,rational) and type(b,float) then |
| > | b_rational:=op(1,b)*(10^(op(2,b))); |
| > | if a >= b_rational then true else false fi; |
| > | elif type(a,float) and type(b,rational) then |
| > | a_rational:=op(1,a)*(10^(op(2,a))); |
| > | if a_rational >= b then true else false fi; |
| > | elif type(a,float) and type(b,float) then |
| > | a_rational:=op(1,a)*(10^(op(2,a))); |
| > | b_rational:=op(1,b)*(10^(op(2,b))); |
| > | if a_rational >= b_rational then true else false fi; |
| > | fi; |
| > | end: |
-----------------------------------------------------------------------------
PROCEDURE - IS_IN
This procedure takes two parameters. It tests to see whether the interval,
or numeric x is contained in the interval a. If a is a float then
the procedure constructs an interval to test.
The operator &greater/Interval_greater is used here.
| > | is_in:=proc(x,a) |
| > |
| > | if type(a,interval) then |
| > | if x=FAIL then FAIL |
| > | elif x=infinity then member(infinity,{a[1],a[2]}) |
| > | elif x=-infinity then member(-infinity,{a[1],a[2]}) |
| > | elif not(type(x,{numeric,interval})) then |
| > | ERROR(`first argument must be a numeric or an interval`) |
| > | elif a=[] and x<>[] then false |
| > | elif a[1]=FAIL or a[2]=FAIL then FAIL |
| > | elif type(x,numeric) then |
| > | intpakX_greater(x,a[1]) and intpakX_greater(a[2],x) |
| > | elif type(x,interval) then |
| > | if x=[] then false #???????????? |
| > | elif x[1]=FAIL or x[2]=FAIL then FAIL |
| > | else intpakX_greater(x[1],a[1]) and intpakX_greater(a[2],x[2]) |
| > | fi; |
| > | fi; |
| > | elif a=infinity then evalb(x=infinity) |
| > | elif a=-infinity then evalb(x=-infinity) |
| > | elif a=FAIL then FAIL |
| > | elif type(a,numeric) then |
| > | if x=FAIL then FAIL |
| > | elif x=infinity or x=-infinity then false |
| > | elif type(x,numeric) then |
| > | intpakX_greater(a,x) and intpakX_greater(x,a); |
| > | elif type(x,interval) then |
| > | if x=[] then false |
| > | else intpakX_greater(x[1],a) and intpakX_greater(a,x[2]) |
| > | fi; |
| > | else ERROR(`first argument must be a numeric or an interval`) |
| > | fi; |
| > | else ERROR(`second argument must be an interval or a numeric`) |
| > | fi; |
| > | end: |
-----------------------------------------------------------------------------
PROCEDURE - CONSTRUCT
This procedure can accept a single argument, to construct a degenerate
interval, two arguments, to construct an interval form the low/high
endpoints. The option `rounded` can be entered as the last argument in
each case to construct a rounded interval.
Special type checking is performed to find (Evalf(-infinity)). Where it
occurs as an argument, Evalf(-infinity) is always rounded to -infinity, to
eliminate Evalf(-infinity) from being an interval endpoint.
Suitable arguments are numeric, FAIL, +-infinity and constants.
| > | construct:=proc() |
| > | local p,q,Constants; |
| > | if member(false,{args}) then |
| > | ERROR(`false and true ar not accepted as arguments`); |
| > | fi; |
| > |
| > | if member(FAIL,{args}) then RETURN([FAIL,FAIL]) fi; |
| > | Constants:={Pi,gamma,Catalan}; |
| > |
| > | if nargs = 3 then |
| > | if args[3] = 'rounded' and |
| > | type(args[1],num_or_FAIL) and |
| > | type(args[2],num_or_FAIL) then |
| > | p:=Evalf(args[1]); |
| > | q:=Evalf(args[2]); |
| > | [rd(intpakX_min(p,q)),ru(intpakX_max(p,q))] |
| > | else ERROR(`first,second args must be numeric,third arg must be "rounded"`) |
| > | fi |
| > |
| > | elif nargs=2 then |
| > | if args[2]=rounded and |
| > | type(args[1],num_or_FAIL) then |
| > | [rd(Evalf(args[1])),ru(Evalf(args[1]))] |
| > | elif type(args[1],num_or_FAIL) and |
| > | type(args[2],num_or_FAIL) then |
| > |
| > | p:=(args[1]); |
| > | q:=(args[2]); |
| > | if member(p,Constants) and member(q,Constants) |
| > | then construct(p,q,rounded); |
| > | elif member(p,Constants) then |
| > | if Evalf(p) < q then |
| > | if intpakX_greater(Evalf(q),q) then |
| > | [rd(Evalf(p)),Evalf(q)] |
| > | else [rd(Evalf(p)),ru(Evalf(q))] |
| > | fi; |
| > | elif Evalf(p) > q then |
| > | if intpakX_greater(q,Evalf(q)) then |
| > | [Evalf(q),ru(Evalf(p))] |
| > | else [rd(Evalf(q)),ru(Evalf(p))] |
| > | fi; |
| > | else construct(p,q,rounded); |
| > | fi; |
| > | elif member(q,Constants) then |
| > | if Evalf(q) < p then |
| > | if intpakX_greater(Evalf(p),p) then |
| > | [rd(Evalf(q)),Evalf(p)] |
| > | else [rd(Evalf(q)),ru(Evalf(p))] |
| > | fi; |
| > | elif Evalf(q) > p then |
| > | if intpakX_greater(p,Evalf(p)) then |
| > | [Evalf(p),ru(Evalf(q))] |
| > | else [rd(Evalf(p)),ru(Evalf(q))] |
| > | fi; |
| > | else construct(p,q,rounded); |
| > | fi; |
| > | elif intpakX_greater(q,p) then |
| > | if intpakX_greater(Evalf(q),q) and |
| > | intpakX_greater(p,Evalf(p)) then |
| > | [Evalf(p),Evalf(q)]; |
| > | elif intpakX_greater(Evalf(q),q) then |
| > | [rd(Evalf(p)),Evalf(q)]; |
| > | elif intpakX_greater(p,Evalf(p)) then |
| > | [Evalf(p),ru(Evalf(q))]; |
| > | else [rd(Evalf(p)),ru(Evalf(q))]; |
| > | fi; |
| > | else |
| > | if intpakX_greater(Evalf(p),p) and |
| > | intpakX_greater(q,Evalf(q)) then |
| > | [Evalf(q),Evalf(p)]; |
| > | elif intpakX_greater(Evalf(p),p) then |
| > | [rd(Evalf(q)),Evalf(p)]; |
| > | elif intpakX_greater(q,Evalf(q)) then |
| > | [Evalf(q),ru(Evalf(p))]; |
| > | else [rd(Evalf(q)),ru(Evalf(p))]; |
| > | fi; |
| > |
| > | fi |
| > | else ERROR(`incorrect arguments entered`) |
| > | fi |
| > |
| > | elif nargs=1 then |
| > | if type(args[1],num_or_FAIL) then |
| > | p:=args[1]; |
| > | if member(p,Constants) then |
| > | [rd(Evalf(p)),ru(Evalf(p))]; |
| > | elif is_in(p,Evalf(p)) then |
| > | [Evalf(p),Evalf(p)]; |
| > | elif intpakX_greater(p,Evalf(p)) then |
| > | [Evalf(p),ru(Evalf(p))]; |
| > | else [rd(Evalf(p)),Evalf(p)]; |
| > | fi; |
| > | else ERROR(`incorrect arguments entered`) |
| > | fi |
| > | else ERROR(`maximum of three arguments accepted`) |
| > | fi |
| > | end: |
###########################################################################
Logarithm and Exponential Functions:
PROCEDURE &EXP/INTERVAL_EXP
expinfinity is only called from &exp. It deals with FAIL and +/- infinity.
Like most of the other subroutines &exp takes floating point intervals
or numerics (which are converted into intervals).
| > | expinfinity:=proc(x): |
| > | if x=FAIL then FAIL |
| > | elif x=infinity then infinity |
| > | elif x=-infinity then 0 |
| > | else Evalf(exp(x)) |
| > | fi: |
| > | end: |
| > | Interval_exp:=proc(x): |
| > | if type(x,'interval') then |
| > | if x=[] then [] |
| > | elif x[1]=FAIL or x[2]=FAIL then [FAIL,FAIL] |
| > | else [rd(expinfinity(x[1])),ru(expinfinity(x[2]))] |
| > | fi: |
| > | elif type(x,'num_or_FAIL') then |
| > | Interval_exp(construct(x)) |
| > | else |
Return unevaluated.
| > | 'Interval_exp(x)' |
ERROR(`floating point interval or scalar arguments are required`)
| > | fi: end: |
---------------------------------------------------------------------------
PROCEDURE &LN
| > | infinityln:=proc(x) |
| > | if x=FAIL then FAIL |
| > | elif x=infinity then infinity |
| > | elif x=0 then -infinity |
| > | elif (min(x,0)=x ) then FAIL |
The above line returns FAIL, as opposed to Maple`s ln function which
returns an ERROR message stating that a singularity has been encountered.
| > | else ln(x) |
| > | fi: |
| > | end: |
| > | Interval_ln:=proc(x) |
| > | if type(x,'interval') then |
| > | if x=[] then [] |
| > | else [rd(infinityln(x[1])),ru(infinityln(x[2]))] |
| > | fi: |
| > | elif type(x,'num_or_FAIL') then Interval_ln(construct(x)) |
| > | else |
Return unevaluated
| > | 'Interval_ln(x)' |
ERROR (`floating point intervals or scalars required`)
| > | fi: |
| > | end: |
---------------------------------------------------------------------------
###########################################################################
BASIC INTERVAL ARITHMETIC OPERATIONS
All the interval arithmetic subroutines perform type checking.
They accept scalars (of type numeric) or intervals (floating point).
Floating point and integer scalars are made into intervals.
-------------------------------------------------------------------------------
The ***infinity subroutines correct any problems that may occur with
infinite and FAIL results
-------------------------------------------------------------------------------
PROCEDURE &+
| > | addinfinity:=proc(x,y) |
| > | if (x=infinity and y=(-infinity)) or (x=(-infinity) and y=infinity) |
| > | then FAIL |
| > | elif x=FAIL or y=FAIL then FAIL |
| > | elif x=infinity or y=infinity then infinity |
| > | elif x=-infinity or y=-infinity then -infinity |
| > | else x+y |
| > | fi: |
| > | end: |
| > | Interval_add:=proc(a,b) |
| > | if type(a,'interval') and type(b,'interval') then |
| > | if a=[] or b=[] then [] |
| > | else [rd(addinfinity(a[1],b[1])),ru(addinfinity(a[2],b[2]))] |
| > | fi: |
| > | elif type(a,'interval') and type(b,'num_or_FAIL') then Interval_add(a,construct(b)) |
| > | elif type(a,'num_or_FAIL') and type(b,'interval') then Interval_add(construct(a),b) |
| > | elif type(a,'num_or_FAIL') and type(b,'num_or_FAIL') |
| > | then Interval_add(construct(a),construct(b)) |
| > | else |
Return unevaluated
| > | 'Interval_add'(a,b) |
ERROR(`floating point interval and scalar arguments required`)
| > | fi: |
| > | end: |
-------------------------------------------------------------------------------
PROCEDURE &-
| > | subtractinfinity:=proc(x,y) |
| > | if (x=infinity and y=(infinity)) or (x=(-infinity) and y=(-infinity)) |
| > | then FAIL |
| > | elif x=FAIL or y=FAIL then FAIL |
| > | elif x=infinity then infinity |
| > | elif y=infinity then -infinity |
| > | elif x=-infinity then -infinity |
| > | elif y=-infinity then infinity |
| > | else x-y |
| > | fi: |
| > |