Initial population.

1 files changed, 600 insertions, 0 deletions

diff --git a/src/algebra/numquad.spad.pamphlet b/src/algebra/numquad.spad.pamphlet
new file mode 100644
index 00000000..79987429
--- /dev/null
+++ b/src/algebra/numquad.spad.pamphlet
@@ -0,0 +1,600 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra numquad.spad}
+\author{Yurij Baransky}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package NUMQUAD NumericalQuadrature}
+<<package NUMQUAD NumericalQuadrature>>=
+)abbrev package NUMQUAD NumericalQuadrature
+++ Author: Yurij A. Baransky
+++ Date Created: October 90
+++ Date Last Updated: October 90
+++ Basic Operations:
+++ Related Constructors:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ References:
+++ Description:
+++ This suite of routines performs numerical quadrature using
+++ algorithms derived from the basic trapezoidal rule. Because
+++ the error term of this rule contains only even powers of the
+++ step size (for open and closed versions), fast convergence
+++ can be obtained if the integrand is sufficiently smooth.
+++
+++ Each routine returns a Record of type TrapAns, which contains\indent{3}
+++ \newline value (\spadtype{Float}):\tab{20}  estimate of the integral
+++ \newline error (\spadtype{Float}):\tab{20}  estimate of the error in the computation
+++ \newline totalpts (\spadtype{Integer}):\tab{20} total number of function evaluations
+++ \newline success (\spadtype{Boolean}):\tab{20} if the integral was computed within the user specified error criterion
+++ \indent{0}\indent{0}
+++ To produce this estimate, each routine generates an internal
+++ sequence of sub-estimates, denoted by {\em S(i)}, depending on the
+++ routine, to which the various convergence criteria are applied.
+++ The user must supply a relative accuracy, \spad{eps_r}, and an absolute
+++ accuracy, \spad{eps_a}. Convergence is obtained when either
+++ \center{\spad{ABS(S(i) - S(i-1)) < eps_r * ABS(S(i-1))}}
+++ \center{or    \spad{ABS(S(i) - S(i-1)) < eps_a}}
+++ are true statements.
+++
+++ The routines come in three families and three flavors:
+++ \newline\tab{3} closed:\tab{20}romberg,\tab{30}simpson,\tab{42}trapezoidal
+++ \newline\tab{3} open:  \tab{20}rombergo,\tab{30}simpsono,\tab{42}trapezoidalo
+++ \newline\tab{3} adaptive closed:\tab{20}aromberg,\tab{30}asimpson,\tab{42}atrapezoidal
+++ \par
+++ The {\em S(i)} for the trapezoidal family is the value of the
+++ integral using an equally spaced absicca trapezoidal rule for
+++ that level of refinement.
+++ \par
+++ The {\em S(i)} for the simpson family is the value of the integral
+++ using an equally spaced absicca simpson rule for that level of
+++ refinement.
+++ \par
+++ The {\em S(i)} for the romberg family is the estimate of the integral
+++ using an equally spaced absicca romberg method. For
+++ the \spad{i}-th level, this is an appropriate combination of all the
+++ previous trapezodial estimates so that the error term starts
+++ with the \spad{2*(i+1)} power only.
+++ \par
+++ The three families come in a closed version, where the formulas
+++ include the endpoints, an open version where the formulas do not
+++ include the endpoints and an adaptive version, where the user
+++ is required to input the number of subintervals over which the
+++ appropriate closed family integrator will apply with the usual
+++ convergence parmeters for each subinterval. This is useful
+++ where a large number of points are needed only in a small fraction
+++ of the entire domain.
+++ \par
+++ Each routine takes as arguments:
+++ \newline f\tab{10} integrand
+++ \newline a\tab{10} starting point
+++ \newline b\tab{10} ending point
+++ \newline \spad{eps_r}\tab{10} relative error
+++ \newline \spad{eps_a}\tab{10} absolute error
+++ \newline \spad{nmin} \tab{10} refinement level when to start checking for convergence (> 1)
+++ \newline \spad{nmax} \tab{10} maximum level of refinement
+++ \par
+++ The adaptive routines take as an additional parameter
+++ \newline \spad{nint}\tab{10} the number of independent intervals to apply a closed
+++  family integrator of the same name.
+++ \par Notes:
+++ \newline Closed family level i uses \spad{1 + 2**i} points.
+++ \newline Open family level i uses \spad{1 + 3**i} points.
+NumericalQuadrature(): Exports == Implementation where
+  L        ==> List
+  V        ==> Vector
+  I        ==> Integer
+  B        ==> Boolean
+  E        ==> OutputForm
+  F       ==> Float
+  PI       ==> PositiveInteger
+  OFORM    ==> OutputForm
+  TrapAns  ==> Record(value:F, error:F, totalpts:I, success:B )
+
+  Exports ==> with
+   aromberg     : (F -> F,F,F,F,F,I,I,I) -> TrapAns
+     ++ aromberg(fn,a,b,epsrel,epsabs,nmin,nmax,nint)
+     ++ uses the adaptive romberg method to numerically integrate function
+     ++ \spad{fn} over the closed interval from \spad{a} to \spad{b},
+     ++ with relative accuracy \spad{epsrel} and absolute accuracy
+     ++ \spad{epsabs}, with the refinement levels for convergence checking
+     ++ vary from \spad{nmin} to \spad{nmax}, and where \spad{nint}
+     ++ is the number of independent intervals to apply the integrator.
+     ++ The value returned is a record containing the value of the integral,
+     ++ the estimate of the error in the computation, the total number of
+     ++ function evaluations, and either a boolean value which is true if
+     ++ the integral was computed within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+   asimpson     : (F -> F,F,F,F,F,I,I,I) -> TrapAns
+     ++ asimpson(fn,a,b,epsrel,epsabs,nmin,nmax,nint) uses the
+     ++ adaptive simpson method to numerically integrate function \spad{fn}
+     ++ over the closed interval from \spad{a} to \spad{b}, with relative
+     ++ accuracy \spad{epsrel} and absolute accuracy \spad{epsabs}, with the
+     ++ refinement levels for convergence checking vary from \spad{nmin}
+     ++ to \spad{nmax}, and where \spad{nint} is the number of independent
+     ++ intervals to apply the integrator. The value returned is a record
+     ++ containing the value of the integral, the estimate of the error in
+     ++ the computation, the total number of function evaluations, and
+     ++ either a boolean value which is true if the integral was computed
+     ++ within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+   atrapezoidal : (F -> F,F,F,F,F,I,I,I) -> TrapAns
+     ++ atrapezoidal(fn,a,b,epsrel,epsabs,nmin,nmax,nint) uses the
+     ++ adaptive trapezoidal method to numerically integrate function
+     ++ \spad{fn} over the closed interval from \spad{a} to \spad{b}, with
+     ++ relative accuracy \spad{epsrel} and absolute accuracy \spad{epsabs},
+     ++ with the refinement levels for convergence checking vary from
+     ++ \spad{nmin} to \spad{nmax}, and where \spad{nint} is the number
+     ++ of independent intervals to apply the integrator. The value returned
+     ++ is a record containing the value of the integral, the estimate of
+     ++ the error in the computation, the total number of function
+     ++ evaluations, and either a boolean value which is true if
+     ++ the integral was computed within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+   romberg      : (F -> F,F,F,F,F,I,I) -> TrapAns
+     ++ romberg(fn,a,b,epsrel,epsabs,nmin,nmax) uses the romberg
+     ++ method to numerically integrate function \spadvar{fn} over the closed
+     ++ interval \spad{a} to \spad{b}, with relative accuracy \spad{epsrel}
+     ++ and absolute accuracy \spad{epsabs}, with the refinement levels
+     ++ for convergence checking vary from \spad{nmin} to \spad{nmax}.
+     ++ The value returned is a record containing the value
+     ++ of the integral, the estimate of the error in the computation, the
+     ++ total number of function evaluations, and either a boolean value
+     ++ which is true if the integral was computed within the user specified
+     ++ error criterion. See \spadtype{NumericalQuadrature} for details.
+   simpson      : (F -> F,F,F,F,F,I,I) -> TrapAns
+     ++ simpson(fn,a,b,epsrel,epsabs,nmin,nmax) uses the simpson
+     ++ method to numerically integrate function \spad{fn} over the closed
+     ++ interval \spad{a} to \spad{b}, with
+     ++ relative accuracy \spad{epsrel} and absolute accuracy \spad{epsabs},
+     ++ with the refinement levels for convergence checking vary from
+     ++ \spad{nmin} to \spad{nmax}. The value returned
+     ++ is a record containing the value of the integral, the estimate of
+     ++ the error in the computation, the total number of function
+     ++ evaluations, and either a boolean value which is true if
+     ++ the integral was computed within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+   trapezoidal  : (F -> F,F,F,F,F,I,I) -> TrapAns
+     ++ trapezoidal(fn,a,b,epsrel,epsabs,nmin,nmax) uses the
+     ++ trapezoidal method to numerically integrate function \spadvar{fn} over
+     ++ the closed interval \spad{a} to \spad{b}, with relative accuracy
+     ++ \spad{epsrel} and absolute accuracy \spad{epsabs}, with the
+     ++ refinement levels for convergence checking vary
+     ++ from \spad{nmin} to \spad{nmax}. The value
+     ++ returned is a record containing the value of the integral, the
+     ++ estimate of the error in the computation, the total number of
+     ++ function evaluations, and either a boolean value which is true
+     ++ if the integral was computed within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+   rombergo     : (F -> F,F,F,F,F,I,I) -> TrapAns
+     ++ rombergo(fn,a,b,epsrel,epsabs,nmin,nmax) uses the romberg
+     ++ method to numerically integrate function \spad{fn} over
+     ++ the open interval from \spad{a} to \spad{b}, with
+     ++ relative accuracy \spad{epsrel} and absolute accuracy \spad{epsabs},
+     ++ with the refinement levels for convergence checking vary from
+     ++ \spad{nmin} to \spad{nmax}. The value returned
+     ++ is a record containing the value of the integral, the estimate of
+     ++ the error in the computation, the total number of function
+     ++ evaluations, and either a boolean value which is true if
+     ++ the integral was computed within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+   simpsono     : (F -> F,F,F,F,F,I,I) -> TrapAns
+     ++ simpsono(fn,a,b,epsrel,epsabs,nmin,nmax) uses the
+     ++ simpson method to numerically integrate function \spad{fn} over
+     ++ the open interval from \spad{a} to \spad{b}, with
+     ++ relative accuracy \spad{epsrel} and absolute accuracy \spad{epsabs},
+     ++ with the refinement levels for convergence checking vary from
+     ++ \spad{nmin} to \spad{nmax}. The value returned
+     ++ is a record containing the value of the integral, the estimate of
+     ++ the error in the computation, the total number of function
+     ++ evaluations, and either a boolean value which is true if
+     ++ the integral was computed within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+   trapezoidalo : (F -> F,F,F,F,F,I,I) -> TrapAns
+     ++ trapezoidalo(fn,a,b,epsrel,epsabs,nmin,nmax) uses the
+     ++ trapezoidal method to numerically integrate function \spad{fn}
+     ++ over the open interval from \spad{a} to \spad{b}, with
+     ++ relative accuracy \spad{epsrel} and absolute accuracy \spad{epsabs},
+     ++ with the refinement levels for convergence checking vary from
+     ++ \spad{nmin} to \spad{nmax}. The value returned
+     ++ is a record containing the value of the integral, the estimate of
+     ++ the error in the computation, the total number of function
+     ++ evaluations, and either a boolean value which is true if
+     ++ the integral was computed within the user specified error criterion.
+     ++ See \spadtype{NumericalQuadrature} for details.
+
+  Implementation ==> add
+   trapclosed : (F -> F,F,F,F,I) -> F
+   trapopen   : (F -> F,F,F,F,I) -> F
+   import OutputPackage
+
+---------------------------------------------------
+
+   aromberg(func,a,b,epsrel,epsabs,nmin,nmax,nint) ==
+      ans  : TrapAns
+      sum  : F := 0.0
+      err  : F := 0.0
+      pts  : I  := 1
+      done : B  := true
+      hh   : F := (b-a) / nint
+      x1   : F := a
+      x2   : F := a + hh
+      io   : L OFORM := [x1::E,x2::E]
+      i    : I
+      for i in 1..nint repeat
+         ans := romberg(func,x1,x2,epsrel,epsabs,nmin,nmax)
+         if (not ans.success) then
+           io.1 := x1::E
+           io.2 := x2::E
+           print blankSeparate cons("accuracy not reached in interval"::E,io)
+         sum  := sum + ans.value
+         err  := err + abs(ans.error)
+         pts  := pts + ans.totalpts-1
+         done := (done and ans.success)
+         x1   := x2
+         x2   := x2 + hh
+      return( [sum , err , pts , done] )
+
+---------------------------------------------------
+
+   asimpson(func,a,b,epsrel,epsabs,nmin,nmax,nint) ==
+      ans  : TrapAns
+      sum  : F := 0.0
+      err  : F := 0.0
+      pts  : I  := 1
+      done : B  := true
+      hh   : F := (b-a) / nint
+      x1   : F := a
+      x2   : F := a + hh
+      io   : L OFORM := [x1::E,x2::E]
+      i    : I
+      for i in 1..nint repeat
+         ans := simpson(func,x1,x2,epsrel,epsabs,nmin,nmax)
+         if (not ans.success) then
+           io.1 := x1::E
+           io.2 := x2::E
+           print blankSeparate cons("accuracy not reached in interval"::E,io)
+         sum  := sum + ans.value
+         err  := err + abs(ans.error)
+         pts  := pts + ans.totalpts-1
+         done := (done and ans.success)
+         x1   := x2
+         x2   := x2 + hh
+      return( [sum , err , pts , done] )
+
+---------------------------------------------------
+
+   atrapezoidal(func,a,b,epsrel,epsabs,nmin,nmax,nint) ==
+      ans  : TrapAns
+      sum  : F := 0.0
+      err  : F := 0.0
+      pts  : I  := 1
+      i    : I
+      done : B := true
+      hh   : F := (b-a) / nint
+      x1   : F := a
+      x2   : F := a + hh
+      io   : L OFORM := [x1::E,x2::E]
+      for i in 1..nint repeat
+         ans := trapezoidal(func,x1,x2,epsrel,epsabs,nmin,nmax)
+         if (not ans.success) then
+           io.1 := x1::E
+           io.2 := x2::E
+           print blankSeparate cons("accuracy not reached in interval"::E,io)
+         sum  := sum + ans.value
+         err  := err + abs(ans.error)
+         pts  := pts + ans.totalpts-1
+         done := (done and ans.success)
+         x1   := x2
+         x2   := x2 + hh
+      return( [sum , err , pts , done] )
+
+---------------------------------------------------
+
+   romberg(func,a,b,epsrel,epsabs,nmin,nmax) ==
+      length : F := (b-a)
+      delta  : F := length
+      newsum : F := 0.5 * length * (func(a)+func(b))
+      newest : F := 0.0
+      oldsum : F := 0.0
+      oldest : F := 0.0
+      change : F := 0.0
+      qx1    : F := newsum
+      table  : V F := new((nmax+1)::PI,0.0)
+      n      : I  := 1
+      pts    : I  := 1
+      four   : I
+      j      : I
+      i      : I
+      if (nmin < 2) then
+         output("romberg: nmin to small (nmin > 1) nmin = ",nmin::E)
+         return([0.0,0.0,0,false])
+      if (nmax < nmin) then
+         output("romberg: nmax < nmin : nmax = ",nmax::E)
+         output("                       nmin = ",nmin::E)
+         return([0.0,0.0,0,false])
+      if (a = b) then
+         output("romberg: integration limits are equal  = ",a::E)
+         return([0.0,0.0,1,true])
+      if (epsrel < 0.0) then
+         output("romberg: eps_r < 0.0            eps_r  = ",epsrel::E)
+         return([0.0,0.0,0,false])
+      if (epsabs < 0.0) then
+         output("romberg: eps_a < 0.0            eps_a  = ",epsabs::E)
+         return([0.0,0.0,0,false])
+      for n in 1..nmax repeat
+         oldsum := newsum
+         newsum := trapclosed(func,a,delta,oldsum,pts)
+         newest := (4.0 * newsum - oldsum) / 3.0
+         four   := 4
+         table(n) := newest
+         for j in 2..n repeat
+            i        := n+1-j
+            four     := four * 4
+            table(i) := table(i+1) + (table(i+1)-table(i)) / (four-1)
+         if n > nmin then
+            change := abs(table(1) - qx1)
+            if change < abs(epsrel*qx1) then
+               return( [table(1) , change , 2*pts+1 , true] )
+            if change < epsabs then
+               return( [table(1) , change , 2*pts+1 , true] )
+         oldsum := newsum
+         oldest := newest
+         delta  := 0.5*delta
+         pts    := 2*pts
+         qx1    := table(1)
+      return( [table(1) , 1.25*change , pts+1 ,false] )
+
+---------------------------------------------------
+
+   simpson(func,a,b,epsrel,epsabs,nmin,nmax) ==
+      length : F := (b-a)
+      delta  : F := length
+      newsum : F := 0.5*(b-a)*(func(a)+func(b))
+      newest : F := 0.0
+      oldsum : F := 0.0
+      oldest : F := 0.0
+      change : F := 0.0
+      n      : I  := 1
+      pts    : I  := 1
+      if (nmin < 2) then
+         output("simpson: nmin to small (nmin > 1) nmin = ",nmin::E)
+         return([0.0,0.0,0,false])
+      if (nmax < nmin) then
+         output("simpson: nmax < nmin : nmax = ",nmax::E)
+         output("                       nmin = ",nmin::E)
+         return([0.0,0.0,0,false])
+      if (a = b) then
+         output("simpson: integration limits are equal  = ",a::E)
+         return([0.0,0.0,1,true])
+      if (epsrel < 0.0) then
+         output("simpson: eps_r < 0.0 : eps_r = ",epsrel::E)
+         return([0.0,0.0,0,false])
+      if (epsabs < 0.0) then
+         output("simpson: eps_a < 0.0 : eps_a = ",epsabs::E)
+         return([0.0,0.0,0,false])
+      for n in 1..nmax repeat
+         oldsum := newsum
+         newsum := trapclosed(func,a,delta,oldsum,pts)
+         newest := (4.0 * newsum - oldsum) / 3.0
+         if n > nmin then
+            change := abs(newest-oldest)
+            if change < abs(epsrel*oldest) then
+               return( [newest , 1.25*change , 2*pts+1 , true] )
+            if change < epsabs then
+               return( [newest , 1.25*change , 2*pts+1 , true] )
+         oldsum := newsum
+         oldest := newest
+         delta  := 0.5*delta
+         pts    := 2*pts
+      return( [newest , 1.25*change , pts+1 ,false] )
+
+---------------------------------------------------
+
+   trapezoidal(func,a,b,epsrel,epsabs,nmin,nmax) ==
+      length : F := (b-a)
+      delta  : F := length
+      newsum : F := 0.5*(b-a)*(func(a)+func(b))
+      change : F := 0.0
+      oldsum : F
+      n      : I  := 1
+      pts    : I  := 1
+      if (nmin < 2) then
+         output("trapezoidal: nmin to small (nmin > 1) nmin = ",nmin::E)
+         return([0.0,0.0,0,false])
+      if (nmax < nmin) then
+         output("trapezoidal: nmax < nmin : nmax = ",nmax::E)
+         output("                           nmin = ",nmin::E)
+         return([0.0,0.0,0,false])
+      if (a = b) then
+         output("trapezoidal: integration limits are equal  = ",a::E)
+         return([0.0,0.0,1,true])
+      if (epsrel < 0.0) then
+         output("trapezoidal: eps_r < 0.0 : eps_r = ",epsrel::E)
+         return([0.0,0.0,0,false])
+      if (epsabs < 0.0) then
+         output("trapezoidal: eps_a < 0.0 : eps_a = ",epsabs::E)
+         return([0.0,0.0,0,false])
+      for n in 1..nmax repeat
+         oldsum := newsum
+         newsum := trapclosed(func,a,delta,oldsum,pts)
+         if n > nmin then
+            change := abs(newsum-oldsum)
+            if change < abs(epsrel*oldsum) then
+               return( [newsum , 1.25*change , 2*pts+1 , true] )
+            if change < epsabs then
+               return( [newsum , 1.25*change , 2*pts+1 , true] )
+         delta := 0.5*delta
+         pts   := 2*pts
+      return( [newsum , 1.25*change , pts+1 ,false] )
+
+---------------------------------------------------
+
+   rombergo(func,a,b,epsrel,epsabs,nmin,nmax) ==
+      length : F := (b-a)
+      delta  : F := length / 3.0
+      newsum : F := length * func( 0.5*(a+b) )
+      newest : F := 0.0
+      oldsum : F := 0.0
+      oldest : F := 0.0
+      change : F := 0.0
+      qx1    : F := newsum
+      table  : V F := new((nmax+1)::PI,0.0)
+      four   : I
+      j      : I
+      i      : I
+      n      : I  := 1
+      pts    : I  := 1
+      for n in 1..nmax repeat
+         oldsum   := newsum
+         newsum   := trapopen(func,a,delta,oldsum,pts)
+         newest   := (9.0 * newsum - oldsum) / 8.0
+         table(n) := newest
+         nine     := 9
+         output(newest::E)
+         for j in 2..n repeat
+            i        := n+1-j
+            nine     := nine * 9
+            table(i) := table(i+1) + (table(i+1)-table(i)) / (nine-1)
+         if n > nmin then
+            change := abs(table(1) - qx1)
+            if change < abs(epsrel*qx1) then
+               return( [table(1) , 1.5*change , 3*pts , true] )
+            if change < epsabs then
+               return( [table(1) , 1.5*change , 3*pts , true] )
+         output(table::E)
+         oldsum := newsum
+         oldest := newest
+         delta  := delta / 3.0
+         pts    := 3*pts
+         qx1    := table(1)
+      return( [table(1) , 1.5*change , pts ,false] )
+
+---------------------------------------------------
+
+   simpsono(func,a,b,epsrel,epsabs,nmin,nmax) ==
+      length : F := (b-a)
+      delta  : F := length / 3.0
+      newsum : F := length * func( 0.5*(a+b) )
+      newest : F := 0.0
+      oldsum : F := 0.0
+      oldest : F := 0.0
+      change : F := 0.0
+      n      : I  := 1
+      pts    : I  := 1
+      for n in 1..nmax repeat
+         oldsum := newsum
+         newsum := trapopen(func,a,delta,oldsum,pts)
+         newest := (9.0 * newsum - oldsum) / 8.0
+         output(newest::E)
+         if n > nmin then
+            change := abs(newest - oldest)
+            if change < abs(epsrel*oldest) then
+               return( [newest , 1.5*change , 3*pts , true] )
+            if change < epsabs then
+               return( [newest , 1.5*change , 3*pts , true] )
+         oldsum := newsum
+         oldest := newest
+         delta  := delta / 3.0
+         pts    := 3*pts
+      return( [newest , 1.5*change , pts ,false] )
+
+---------------------------------------------------
+
+   trapezoidalo(func,a,b,epsrel,epsabs,nmin,nmax) ==
+      length : F := (b-a)
+      delta  : F := length/3.0
+      newsum : F := length*func( 0.5*(a+b) )
+      change : F := 0.0
+      pts    : I  := 1
+      oldsum : F
+      n      : I
+      for n in 1..nmax repeat
+         oldsum := newsum
+         newsum := trapopen(func,a,delta,oldsum,pts)
+         output(newsum::E)
+         if n > nmin then
+            change := abs(newsum-oldsum)
+            if change < abs(epsrel*oldsum) then
+               return([newsum , 1.5*change , 3*pts , true] )
+            if change < epsabs then
+               return([newsum , 1.5*change , 3*pts , true] )
+         delta := delta / 3.0
+         pts   := 3*pts
+      return([newsum , 1.5*change , pts ,false] )
+
+---------------------------------------------------
+
+   trapclosed(func,start,h,oldsum,numpoints) ==
+      x   : F := start + 0.5*h
+      sum : F := 0.0
+      i   : I
+      for i in 1..numpoints repeat
+          sum := sum + func(x)
+          x   := x + h
+      return( 0.5*(oldsum + sum*h) )
+
+---------------------------------------------------
+
+   trapopen(func,start,del,oldsum,numpoints) ==
+      ddel : F := 2.0*del
+      x   : F := start + 0.5*del
+      sum : F := 0.0
+      i   : I
+      for i in 1..numpoints repeat
+          sum := sum + func(x)
+          x   := x + ddel
+          sum := sum + func(x)
+          x   := x + del
+      return( (oldsum/3.0 + sum*del) )
+
+@
+\section{License}
+<<license>>=
+--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+--All rights reserved.
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+--    - Redistributions of source code must retain the above copyright
+--      notice, this list of conditions and the following disclaimer.
+--
+--    - Redistributions in binary form must reproduce the above copyright
+--      notice, this list of conditions and the following disclaimer in
+--      the documentation and/or other materials provided with the
+--      distribution.
+--
+--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--      names of its contributors may be used to endorse or promote products
+--      derived from this software without specific prior written permission.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+<<*>>=
+<<license>>
+
+<<package NUMQUAD NumericalQuadrature>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}