Re: Turning an "Expression Float" into a Float?

The attached may be useful


(1) -> )r erf2

Nerf3(a:FLOAT,ndigs:PI,nterms:PI):FLOAT ==
-- evaluate the error function
  olds:= digits ndigs
  table:=empty()$TABLE(NNI,FLOAT)
-- expansion 7.1.5 of Abramowitz
  table.0:=a
  for i in 1..nterms repeat table.i:=(-1)*a^2*(2*i-1)*table.(i-1)/(i*(2*i+1))
  ans:=(2/sqrt(%pi)) * reduce(+,parts table)
  digits olds
  ans

   Function declaration Nerf3 : (Float,PositiveInteger,PositiveInteger)
       -> Float has been added to workspace.
                                                                   Type: Void


Nerff(a:FLOAT,ndigs:PI,nterms:PI):FLOAT ==
-- in terms of onefone 7.1.21 of Abramowitz
  olds:= digits ndigs
  ans:= 2*a*exp(-a*a)*onefone(1.0,1.5,a*a,ndigs,nterms)/sqrt(%pi)
  digits olds
  ans

   Function declaration Nerff : (Float,PositiveInteger,PositiveInteger)
       -> Float has been added to workspace.
                                                                   Type: Void



onefone(a:FLOAT,b:FLOAT,z:FLOAT,ndigs:PI,nterms:PI):FLOAT ==
-- the 1F1 confluent hypergeometric function
  olds:= digits ndigs
  table:=empty()$TABLE(NNI,FLOAT)
-- expansion 13.1.2 of Abramowitz
  table.0:=1
  table.1:=a*z/b
  for i in 2..nterms repeat table.i:=(a+i-1)*z*table.(i-1)/((b+i-1)*i)
  ans:=reduce(+,parts table)
  digits olds
  ans

   Function declaration onefone : (Float,Float,Float,PositiveInteger,
      PositiveInteger) -> Float has been added to workspace.
                                                                   Type: Void

expr:=erf(1.0)


   (4)  erf(1.0)
                                                       Type: Expression Float
tower %


   (5)  [erf(1.0)]
                                           Type: List Kernel Expression Float
% 1


   (6)  erf(1.0)
                                                Type: Kernel Expression Float
erfop:=operator %


   (7)  erf
                                                          Type: BasicOperator
evl:LIST EXPR FLOAT -> EXPR FLOAT

                                                                   Type: Void
-- insert an evaluation rule that always uses 50 digits and 500 terms
evl(x) == Nerff(x(1),50,500)

                                                                   Type: Void
--attach it to the erf operator
evaluate(erfop,evl)$BOP1(EXPR FLOAT)

   Compiling function onefone with type (Float,Float,Float,
      PositiveInteger,PositiveInteger) -> Float
   Compiling function Nerff with type (Float,PositiveInteger,
      PositiveInteger) -> Float
   Compiling function evl with type List Expression Float -> Expression
      Float

   (10)  erf
                                                          Type: BasicOperator
erf (1.0)


   (11)  0.8427007929 4971486934
                                                       Type: Expression Float
expr


   (12)  erf(1.0)
                                                       Type: Expression Float
numeric %


   (13)  0.8427007929 4971486934
                                                                  Type: Float

integrate(exp(-x^2),x=0..1)


                +---+
         erf(1)\|%pi
   (14)  ------------
               2
                    Type: Union(f1: OrderedCompletion Expression Integer,...)
numeric %


   (15)  0.7468241328 124270254
                                                                  Type: Float


________________________________________________________________________
The Numerical Algorithms Group Ltd is a company registered in England
and Wales with company number 1249803. The registered office is:
Wilkinson House, Jordan Hill Road, Oxford OX2 8DR, United Kingdom.

This e-mail has been scanned for all viruses by Star. The service is
powered by MessageLabs. 
________________________________________________________________________

_______________________________________________
Axiom-mail mailing list
[email protected]
http://lists.nongnu.org/mailman/listinfo/axiom-mail

Nerf3(a:FLOAT,ndigs:PI,nterms:PI):FLOAT == 
-- evaluate the error function 
  olds:= digits ndigs
  table:=empty()$TABLE(NNI,FLOAT)
-- expansion 7.1.5 of Abramowitz
  table.0:=a
  for i in 1..nterms repeat table.i:=(-1)*a^2*(2*i-1)*table.(i-1)/(i*(2*i+1))
  ans:=(2/sqrt(%pi)) * reduce(+,parts table)
  digits olds
  ans


Nerff(a:FLOAT,ndigs:PI,nterms:PI):FLOAT ==
-- in terms of onefone 7.1.21 of Abramowitz
  olds:= digits ndigs
  ans:= 2*a*exp(-a*a)*onefone(1.0,1.5,a*a,ndigs,nterms)/sqrt(%pi)
  digits olds
  ans



onefone(a:FLOAT,b:FLOAT,z:FLOAT,ndigs:PI,nterms:PI):FLOAT == 
-- the 1F1 confluent hypergeometric function
  olds:= digits ndigs
  table:=empty()$TABLE(NNI,FLOAT)
-- expansion 13.1.2 of Abramowitz
  table.0:=1
  table.1:=a*z/b
  for i in 2..nterms repeat table.i:=(a+i-1)*z*table.(i-1)/((b+i-1)*i)
  ans:=reduce(+,parts table)
  digits olds
  ans

expr:=erf(1.0)
tower %
% 1
erfop:=operator %
evl:LIST EXPR FLOAT -> EXPR FLOAT
-- insert an evaluation rule that always uses 50 digits and 500 terms
evl(x) == Nerff(x(1),50,500)
--attach it to the erf operator
evaluate(erfop,evl)$BOP1(EXPR FLOAT)
erf (1.0)
expr
numeric %

integrate(exp(-x^2),x=0..1)
numeric %