This routine provides Simpson's method for numerical integration. Although FriCAS already provides a Simpson's method, this version has a syntax that will be intuitive to anyone who has used the integrate() function.

)abbrev package SIMPINT SimpsonIntegration
SimpsonIntegration(): Exports == Implementation where
  F        ==> Float
  SF       ==> Segment F
  EF       ==> Expression F
  SBF      ==> SegmentBinding F
  Ans      ==> Record(value:F, error:F)

  Exports ==> with
   simpson : (EF,SBF,F) -> Ans
   simpson : (EF,SBF) -> Ans

  Implementation ==> add
   simpson(func:EF, sbf:SBF, tol:F) ==
     a : F := low(segment(sbf))
     b : F := high(segment(sbf))
     x : EF := variable(sbf) :: EF

     h : F
     k : Integer
     n : Integer

     simps : F
     newsimps : F

     oe : F
     ne : F
     err : F

     sumend : F := retract(eval(func, x, a::EF) + eval(func, x, b::EF))
     sumodd : F := 0.0 :: F
     sumeven : F := 0.0 :: F

     -- First base case -- 2 intervals ----------------
     n := 2
     h := (b-a)/n
     sumeven := sumeven + sumodd
     sumodd := 0.0 :: F

     for k in 1..(n-1) by 2 repeat
       sumodd := sumodd + retract(eval( func, x, (k*h+a)::EF ))

     simps := ( sumend + 2.0*sumeven + 4.0*sumodd )*(h/3.0)

     -- Second base case -- 4 intervals ---------------
     n := n*2
     h := (b-a)/n
     sumeven := sumeven + sumodd
     sumodd := 0.0 :: F

     for k in 1..(n-1) by 2 repeat
       sumodd := sumodd + retract(eval( func, x, (k*h+a)::EF ))

     newsimps := ( sumend + 2.0*sumeven + 4.0*sumodd )*(h/3.0)

     oe := abs(newsimps-simps)                  -- old error
     simps := newsimps

     -- general case -----------------------------------
     while true repeat
       n := n*2
       h := (b-a)/n

       sumeven := sumeven + sumodd
       sumodd := 0.0 :: F

       for k in 1..(n-1) by 2 repeat
         sumodd := sumodd + retract(eval( func, x, (k*h+a)::EF ))

       newsimps := ( sumend + 2.0*sumeven + 4.0*sumodd )*(h/3.0)

       -- This is a check of Richardson's error estimate.
       -- Usually p is approximately 4 for Simpson's rule, but
       -- occasionally convergence is slower

       ne := abs( newsimps - simps )            -- new error

       if ( (ne<oe*2.0) and (oe<ne*16.5) ) then -- Richardson should be ok
         -- p := log(oe/ne)/log(2.0)
         err := ne/(oe/ne-1.0::F)              -- ne/(2^p-1)
       else
         err := ne                              -- otherwise estimate crudely

       oe := ne
       simps := newsimps

       if( err < tol ) then
         break

     [ newsimps, err ]

   simpson(func:EF, sbf:SBF) ==
     simpson( func, sbf, 1.e-6::F )

   Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/6914547373735257695-25px001.spad
      using old system compiler.
   SIMPINT abbreviates package SimpsonIntegration 
------------------------------------------------------------------------
   initializing NRLIB SIMPINT for SimpsonIntegration 
   compiling into NRLIB SIMPINT 
   compiling exported simpson : (Expression Float,SegmentBinding Float,Float) -> Record(value: Float,error: Float)
Time: 0.66 SEC.

   compiling exported simpson : (Expression Float,SegmentBinding Float) -> Record(value: Float,error: Float)
Time: 0.01 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |SimpsonIntegration| REDEFINED

;;;     ***       |SimpsonIntegration| REDEFINED
Time: 0 SEC.

   Cumulative Statistics for Constructor SimpsonIntegration
      Time: 0.67 seconds

   finalizing NRLIB SIMPINT 
   Processing SimpsonIntegration for Browser database:
--->-->SimpsonIntegration(constructor): Not documented!!!!
--->-->SimpsonIntegration((simpson ((Record (: value (Float)) (: error (Float))) (Expression (Float)) (SegmentBinding (Float)) (Float)))): Not documented!!!!
--->-->SimpsonIntegration((simpson ((Record (: value (Float)) (: error (Float))) (Expression (Float)) (SegmentBinding (Float))))): Not documented!!!!
--->-->SimpsonIntegration(): Missing Description
; compiling file "/var/aw/var/LatexWiki/SIMPINT.NRLIB/SIMPINT.lsp" (written 25 JUN 2019 02:33:51 PM):

; /var/aw/var/LatexWiki/SIMPINT.NRLIB/SIMPINT.fasl written
; compilation finished in 0:00:00.128
------------------------------------------------------------------------
   SimpsonIntegration is now explicitly exposed in frame initial 
   SimpsonIntegration will be automatically loaded when needed from 
      /var/aw/var/LatexWiki/SIMPINT.NRLIB/SIMPINT

This simpson() function overloads the already existing function and either may be used. To see available simpson() functions, do:

)display op simpson

There are 3 exposed functions called simpson :
   [1] (Expression(Float),SegmentBinding(Float),Float) -> Record(value
            : Float,error: Float)
            from SimpsonIntegration
   [2] (Expression(Float),SegmentBinding(Float)) -> Record(value: Float
            ,error: Float)
            from SimpsonIntegration
   [3] ((D3 -> D3),D3,D3,D3,D3,Integer,Integer) -> Record(value: D3,
            error: D3,totalpts: Integer,success: Boolean)
            from NumericalQuadrature(D3) if D3 has FPS

To compute an integral using Simpson's rule, pass an expression and a BindingSegment? with the limits. Optionally, you may include a third argument to specify the acceptable error.

The exact integral:

integrate( sin(x), x=0..1 ) :: Expression Float

(1)

Our approximations:

simpson( sin(x), x=0..1 )

(2)

simpson( sin(x), x=0..1, 1.e-10 )

(3)