\documentclass{article}
\usepackage{open-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]
      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]
      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
      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)
      pts    : I  := 1
      four   : 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 negative? epsrel then
         output("romberg: eps_r < 0.0            eps_r  = ",epsrel::E)
         return([0.0,0.0,0,false])
      if negative? epsabs 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
      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 negative? epsrel then
         output("simpson: eps_r < 0.0 : eps_r = ",epsrel::E)
         return([0.0,0.0,0,false])
      if negative? epsabs 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
      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 negative? epsrel then
         output("trapezoidal: eps_r < 0.0 : eps_r = ",epsrel::E)
         return([0.0,0.0,0,false])
      if negative? epsabs 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
      i      : I
      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
      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
      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
      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
      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}