%\documentclass[a4paper,10pt]{article}
  \documentclass[a4paper,english]{article}
  
  \usepackage[utf8]{inputenc}
  \usepackage{a4wide}
  \usepackage{eurosym}
  \usepackage{url}

  \usepackage[colorlinks=true,hyperfigures=true]{hyperref}

  %\usepackage{aeguill}
  
  \usepackage{graphicx}
  \usepackage{babel}
  \makeatother
  
  
  %opening
  \title{FVN Documentation}
  \author{William Daniau}
  
  
  \begin{document}
  
  \maketitle
  
  %\begin{abstract}
  
  %\end{abstract}
  \tableofcontents
  
  \section{Whatis fvn,licence,disclaimer etc}
  \subsection{Whatis fvn}

  fvn is a Fortran95 mathematical module. It provides various usefull subroutine covering linear algebra, numerical integration, least square polynomial, spline interpolation, zero finding, special functions etc.

  Most of the work is done by interfacing Lapack \url{http://www.netlib.org/lapack} which means that Lapack and Blas \url{http://www.netlib.org/blas} must be available on your system for linking fvn. If you use an AMD microprocessor, the good idea is to use ACML ( AMD Core Math Library \url{http://developer.amd.com/acml.jsp} which contains an optimized Blas/Lapack. Fvn also contains a slightly modified version of Quadpack \url{http://www.netlib.org/quadpack} for performing the numerical integration tasks. Finally the fnlib library \url{http://www.netlib.org/fn} has been added for special functions.

  This module has been initially written for the use of the ``Acoustic and microsonic'' group leaded by Sylvain Ballandras in the Time and Frequency Department of institute Femto-ST \url{http://www.femto-st.fr/}.

  
  \subsection{Licence}
  The licence of fvn is free. You can do whatever you want with this code as far as you credit the authors.
  
  \subsubsection*{Authors}
  As of the day this manuel is written there's only one author of fvn :
  ewline
  William Daniau
  ewline
  william.daniau@femto-st.fr
  ewline
  
  \subsection{Disclaimer}
  The usual disclaimer applied : This software is provided AS IS in the hope it will be usefull. Use it at your own risks. The authors should not be taken responsible of anything that may result by the use of this software.
  
  \section{Naming scheme and convention}
  The naming scheme of the routines is as follow :
  \begin{verbatim}

      fvn_*_name()

  \end{verbatim} 

  where * can be s,d,c or z. 

  \begin{itemize}
   \item s is for single precision real (real,real*4,real(4),real(kind=4))
   \item d for double precision real (double precision,real*8,real(8),real(kind=8))
   \item c for single precision complex (complex,complex*8,complex(4),complex(kind=4))
   \item z for double precision complex (double complex,complex*16,complex(8),complex(kind=8))
  \end{itemize}
  In the following description of subroutines parameters, input parameters are followed by (in), output parameters by (out) and parameters which are used as input and modified by the subroutine are followed by (inout).

  For each routine, there is a generic interface (simply remove \verb'_*' in the name), so using the specific routine is not mandatory.

  \section{Linear algebra}
  The linear algebra routines of fvn are an interface to lapack, which make it easier to use.
  \subsection{Matrix inversion}
  \begin{verbatim}

  call fvn_matinv(d,a,inva,status)

  \end{verbatim}
  \begin{itemize}
   \item d (in) is an integer equal to the matrix rank

   \item a (in) is a real or complex matrix. It will remain untouched.
   \item inva (out) is a real or complex matrix which contain the inverse of a at the end of the routine
   \item status (out) is an optional integer equal to zero if something went wrong

  \end{itemize}
  
  \subsubsection*{Example}
  \begin{verbatim}
  program inv
   use fvn 
   implicit none
  
   real(8),dimension(3,3) :: a,inva

  
   call random_number(a)
   a=a*100

   call fvn_matinv(3,a,inva)

   write (*,*) a
   write (*,*)
   write (*,*) inva
   write (*,*)
   write (*,*) matmul(a,inva)
  end program
  \end{verbatim}
  
  
  
  \subsection{Matrix determinants}
  \begin{verbatim}

  det=fvn_det(d,a,status)

  \end{verbatim}
  \begin{itemize}
   \item d (in) is an integer equal to the matrix rank

   \item a (in) is a real or complex matrix. It will remain untouched.
   \item status (out) is an optional integer equal to zero if something went wrong

  \end{itemize}
  
  \subsubsection*{Example}
  \begin{verbatim}
  program det
   use fvn 
   implicit none
  
   real(8),dimension(3,3) :: a
   real(8) :: deta
   integer :: status
  
   call random_number(a)
   a=a*100

   deta=fvn_det(3,a,status)

   write (*,*) a
   write (*,*)
   write (*,*) "Det = ",deta
  end program
  
  \end{verbatim}
  
  
  
  \subsection{Matrix condition}
  \begin{verbatim}

  call fvn_matcon(d,a,rcond,status)

  \end{verbatim}
  \begin{itemize}
   \item d (in) is an integer equal to the matrix rank

   \item a (in) is a real or complex matrix. It will remain untouched.

   \item rcond (out) is a real of same kind as matrix a, it will contain the reciprocal condition number of the matrix

   \item status (out) is an optional integer equal to zero if something went wrong

  \end{itemize}
  
  The reciprocal condition number is evaluated using the 1-norm and is define as in equation \ref{rconddef}
  \begin{equation}
   R = \frac{1}{norm(A)*norm(invA)}
   \label{rconddef}
  \end{equation}
  
  The 1-norm itself is defined as the maximum value of the columns absolute values (modulus for complex) sum as in equation \ref{l1norm}
  \begin{equation}
   L1 = max_j ( \sum_i{\mid A(i,j)\mid}  )
   \label{l1norm}
  \end{equation}
  
  \subsubsection*{Example}
  \begin{verbatim}
  program cond
   use fvn 
   implicit none
  
   real(8),dimension(3,3) :: a
   real(8) :: rcond
   integer :: status
  
   call random_number(a)
   a=a*100
  
   call fvn_d_matcon(3,a,rcond,status)
   write (*,*) a
   write (*,*)
   write (*,*) "Cond = ",rcond
  end program
  
  \end{verbatim}
  
  
  \subsection{Eigenvalues/Eigenvectors}
  \begin{verbatim}

  call fvn_matev(d,a,evala,eveca,status)

  \end{verbatim}
  \begin{itemize}
   \item d (in) is an integer equal to the matrix rank

   \item a (in) is a real or complex matrix. It will remain untouched.

   \item evala (out) is a complex array of same kind as a. It contains the eigenvalues of matrix a
   \item eveca (out) is a complex matrix of same kind as a. Its columns are the eigenvectors of matrix a : eveca(:,j)=jth eigenvector associated with eigenvalue evala(j).

   \item status (out) is an optional integer equal to zero if something went wrong

  \end{itemize}
  
  \subsubsection*{Example}
  \begin{verbatim}
  program eigen
   use fvn 
   implicit none
  
   real(8),dimension(3,3) :: a
   complex(8),dimension(3) :: evala
   complex(8),dimension(3,3) :: eveca
   integer :: status,i,j
  
   call random_number(a)
   a=a*100

   call fvn_matev(3,a,evala,eveca,status)

   write (*,*) a
   write (*,*)
   do i=1,3
      write(*,*) "Eigenvalue ",i,evala(i)
      write(*,*) "Associated Eigenvector :"
      do j=1,3
          write(*,*) eveca(j,i)
      end do
      write(*,*)
   end do
  
  end program
  
  \end{verbatim}
  
  
  \subsection{Sparse solving}
  By interfacing Tim Davis's SuiteSparse from university of Florida \url{http://www.cise.ufl.edu/research/sparse/SuiteSparse/} which is a reference for this kind of problems, fvn provides simple subroutines for solving linear sparse systems.
  
  The provided routines solves the equation $Ax=B$ where A is sparse and given in its triplet form.
  
  \begin{verbatim}

  call fvn_sparse_solve(n,nz,T,Ti,Tj,B,x,status)

  \end{verbatim}
  \begin{itemize}

   \item For this family of subroutine the two letters (zl,zi,dl,di) of the specific interface name decribe the arguments's type. z is for complex(8), d for real(8), l for integer(8) and i for integer(4)

   \item n (in) is an integer equal to the matrix rank
   \item nz (in) is an integer equal to the number of non-zero elements
   \item T(nz) (in) is a complex/real array containing the non-zero elements
   \item Ti(nz),Tj(nz) (in) are the indexes of the corresponding element of T in the original matrix.
   \item B(n) (in) is a complex/real array containing the second member of the equation.
   \item x(n) (out) is a complex/real array containing the solution
   \item status (out) is an integer which contain non-zero is something went wrong
  \end{itemize}
  
  \subsubsection*{Example}
  \begin{verbatim}
  program test_sparse
  
   use fvn
   implicit none
  
   integer(8), parameter :: nz=12
   integer(8), parameter :: n=5
   complex(8),dimension(nz) :: A
   integer(8),dimension(nz) :: Ti,Tj
   complex(8),dimension(n) :: B,x
   integer(8) :: status
  
   A = (/ (2.,0.),(3.,0.),(3.,0.),(-1.,0.),(4.,0.),(4.,0.),(-3.,0.),&
            (1.,0.),(2.,0.),(2.,0.),(6.,0.),(1.,0.) /)
   B = (/ (8.,0.), (45.,0.), (-3.,0.), (3.,0.), (19.,0.)/)
   Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
   Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)

   !specific routine that will be used here
   !call fvn_zl_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
   call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)

   write(*,*) x
  
  end program
  
  
  program test_sparse
  
  use fvn
  implicit none
  
  integer(4), parameter :: nz=12
  integer(4), parameter :: n=5
  real(8),dimension(nz) :: A
  integer(4),dimension(nz) :: Ti,Tj
  real(8),dimension(n) :: B,x
  integer(4) :: status
  
  A = (/ 2.,3.,3.,-1.,4.,4.,-3.,1.,2.,2.,6.,1. /)
  B = (/ 8., 45., -3., 3., 19./)
  Ti = (/ 1,2,1,3,5,2,3,4,5,3,2,5 /)
  Tj = (/ 1,1,2,2,2,3,3,3,3,4,5,5 /)

  !specific routine that will be used here
  !call fvn_di_sparse_solve(n,nz,A,Ti,Tj,B,x,status)
  call fvn_sparse_solve(n,nz,A,Ti,Tj,B,x,status)

  write(*,*) x
  
  end program
  
  
  
  \end{verbatim}

  \subsection{Identity matrix}
  \begin{verbatim}

   I=fvn_*_ident(n)     (*=s,d,c,z)

  \end{verbatim}
  \begin{itemize}
   \item n (in) is an integer equal to the matrix rank
  \end{itemize}

  This function return the identity matrix of rank n, in the type of the left hand side. No generic interface for this one.

  
  
  \section{Interpolation}
  
  \subsection{Quadratic Interpolation}
  fvn provide function for interpolating values of a tabulated function of 1, 2 or 3 variables, for both single and double precision.
  \subsubsection{One variable function}
  \begin{verbatim}

   value=fvn_quad_interpol(x,n,xdata,ydata)

  \end{verbatim}
  \begin{itemize}
   \item x is the real where we want to evaluate the function
   \item n is the number of tabulated values
   \item xdata(n) contains the tabulated coordinates
   \item ydata(n) contains the tabulated function values ydata(i)=y(xdata(i))
  \end{itemize}
  xdata must be strictly increasingly ordered.
  x must be within the range of xdata to actually perform an interpolation, otherwise the resulting value is an extrapolation
  \paragraph*{Example}
  \begin{verbatim}
  program inter1d
  
  use fvn
  implicit none
  
  integer(kind=4),parameter :: ndata=33
  integer(kind=4) :: i,nout
  real(kind=8) :: f,fdata(ndata),h,pi,q,sin,x,xdata(ndata)
  real(kind=8) ::tv
  
  intrinsic sin
  
  f(x)=sin(x)
  
  xdata(1)=0.
  fdata(1)=f(xdata(1))
  h=1./32.
  do i=2,ndata
        xdata(i)=xdata(i-1)+h
        fdata(i)=f(xdata(i))
  end do
  call random_seed()
  call random_number(x)
  
  q=fvn_d_quad_interpol(x,ndata,xdata,fdata)
  
  tv=f(x)
  write(*,*) "x ",x
  write(*,*) "Calculated (real) value :",tv
  write(*,*) "fvn interpolation :",q
  write(*,*) "Relative fvn error :",abs((q-tv)/tv)
  
  end program
  
  \end{verbatim}
  
  
  \subsubsection{Two variables function}
  \begin{verbatim}

  value=fvn_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)

  \end{verbatim}
  \begin{itemize}
   \item x,y are the real coordinates where we want to evaluate the function
   \item nx is the number of tabulated values along x axis
   \item xdata(nx) contains the tabulated x
   \item ny is the number of tabulated values along y axis
   \item ydata(ny) contains the tabulated y
   \item zdata(nx,ny) contains the tabulated function values zdata(i,j)=z(xdata(i),ydata(j))
  \end{itemize}
  xdata and ydata must be strictly increasingly ordered.
  (x,y) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation
  
  \paragraph*{Example}
  
  \begin{verbatim}
  program inter2d
  use fvn
  implicit none
  
  integer(kind=4),parameter  :: nx=21,ny=42
  integer(kind=4) :: i,j
  real(kind=8) :: f,fdata(nx,ny),dble,pi,q,sin,x,xdata(nx),y,ydata(ny)
  real(kind=8) :: tv
  
  intrinsic dble,sin
  
  f(x,y)=sin(x+2.*y)
  do i=1,nx
        xdata(i)=dble(i-1)/dble(nx-1)
  end do
  do i=1,ny
        ydata(i)=dble(i-1)/dble(ny-1)
  end do
  do i=1,nx
        do j=1,ny
              fdata(i,j)=f(xdata(i),ydata(j))
        end do
  end do
  call random_seed()
  call random_number(x)
  call random_number(y)
  
  q=fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,fdata)
  tv=f(x,y)
  
  write(*,*) "x y",x,y
  write(*,*) "Calculated (real) value :",tv
  write(*,*) "fvn interpolation :",q
  write(*,*) "Relative fvn error :",abs((q-tv)/tv)
  
  end program
  
  \end{verbatim}
  
  
  
  \subsubsection{Three variables function}
  \begin{verbatim}

  value=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)

  \end{verbatim}
  \begin{itemize}
   \item x,y,z are the real coordinates where we want to evaluate the function
   \item nx is the number of tabulated values along x axis
   \item xdata(nx) contains the tabulated x
   \item ny is the number of tabulated values along y axis
   \item ydata(ny) contains the tabulated y
   \item nz is the number of tabulated values along z axis
   \item zdata(ny) contains the tabulated z
   \item tdata(nx,ny,nz) contains the tabulated function values tdata(i,j,k)=t(xdata(i),ydata(j),zdata(k))
  \end{itemize}
  xdata, ydata and zdata must be strictly increasingly ordered.
  (x,y,z) must be within the range of xdata and ydata to actually perform an interpolation, otherwise the resulting value is an extrapolation
  
  \paragraph*{Example}
  \begin{verbatim}
  program inter3d
  use fvn
  
  implicit none
  
  integer(kind=4),parameter  :: nx=21,ny=42,nz=18
  integer(kind=4) :: i,j,k
  real(kind=8) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
  real(kind=8) :: tv
  
  intrinsic dble,sin
  
  f(x,y,z)=sin(x+2.*y+3.*z)
  do i=1,nx
        xdata(i)=2.*(dble(i-1)/dble(nx-1))
  end do
  do i=1,ny
        ydata(i)=2.*(dble(i-1)/dble(ny-1))
  end do
  do i=1,nz
        zdata(i)=2.*(dble(i-1)/dble(nz-1))
  end do
  do i=1,nx
        do j=1,ny
              do k=1,nz
                    fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
              end do
        end do
  end do
  call random_seed()
  call random_number(x)
  call random_number(y)
  call random_number(z)
  
  q=fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
  tv=f(x,y,z)
  
  write(*,*) "x y z",x,y,z
  write(*,*) "Calculated (real) value :",tv
  write(*,*) "fvn interpolation :",q
  write(*,*) "Relative fvn error :",abs((q-tv)/tv)
  
  end program
  
  \end{verbatim}
  
  \subsubsection{Utility procedure}
  fvn provides a simple utility procedure to locate the interval in which a value is located in an increasingly ordered array.
  \begin{verbatim}

  call fvn_find_interval(x,i,xdata,n)

  \end{verbatim}
  \begin{itemize}
   \item x (in) the real value to locate
   \item i (out) the resulting indice
   \item xdata(n) (in) increasingly ordered array
   \item n (in) size of the array
  \end{itemize}
  The resulting integer i is as : $xdata(i) <= x < xdata(i+1)$. If $x < xdata(1)$ then $i=0$ is returned. If $x > xdata(n)$ then $i=n$ is returned. Finally if $x=xdata(n)$ then $i=n-1$ is returned.
  
  
  
  \subsection{Akima spline}
  fvn provides Akima spline interpolation and evaluation for both single and double precision real.
  \subsubsection{Interpolation}
  \begin{verbatim}

  call fvn_akima(n,x,y,br,co)

  \end{verbatim}
  \begin{itemize}
   \item n (in) is an integer equal to the number of points
   \item x(n) (in) ,y(n) (in) are the known couples of coordinates
   \item br (out) on output contains a copy of x
   \item co(4,n) (out) is a real matrix containing the 4 coefficients of the Akima interpolation spline for a given interval.
  \end{itemize}
  
  \subsubsection{Evaluation}
  \begin{verbatim}

  y=fvn_spline_eval(x,n,br,co)

  \end{verbatim}
  \begin{itemize}
   \item x (in) is the point where we want to evaluate
   \item n (in) is the number of known points and br(n) (in), co(4,n) (in) \\
  are the outputs of fvn\_x\_akima(n,x,y,br,co) 
  \end{itemize}
  
  \subsubsection{Example}
  In the following example we will use Akima splines to interpolate a sinus function with 30 points between -10 and 10. We then use the evaluation function to calculate the coordinates of 1000 points between -11 and 11, and write a 3 columns file containing : x, calculated sin(x), interpolation evaluation of sin(x).
  
  One can see that the interpolation is very efficient even with only 30 points. Of course as soon as we leave the -10 to 10 interval, the values are extrapolated and thus can lead to very inacurrate values.
  
  \begin{verbatim}
  program akima
   use fvn
   implicit none
  
   integer :: nbpoints,nppoints,i
   real(8),dimension(:),allocatable :: x_d,y_d,breakpoints_d
   real(8),dimension(:,:),allocatable :: coeff_fvn_d
   real(8) :: xstep_d,xp_d,ty_d,fvn_y_d
  
   open(2,file='fvn_akima_double.dat')
   open(3,file='fvn_akima_breakpoints_double.dat')
   nbpoints=30
   allocate(x_d(nbpoints))
   allocate(y_d(nbpoints))
   allocate(breakpoints_d(nbpoints))
   allocate(coeff_fvn_d(4,nbpoints))
  
   xstep_d=20./dfloat(nbpoints)
   do i=1,nbpoints
      x_d(i)=-10.+dfloat(i)*xstep_d
      y_d(i)=dsin(x_d(i))
      write(3,44) (x_d(i),y_d(i))
   end do
   close(3)
  
   call fvn_d_akima(nbpoints,x_d,y_d,breakpoints_d,coeff_fvn_d)
  
   nppoints=1000 
   xstep_d=22./dfloat(nppoints)
   do i=1,nppoints
      xp_d=-11.+dfloat(i)*xstep_d
      ty_d=dsin(xp_d)
      fvn_y_d=fvn_d_spline_eval(xp_d,nbpoints-1,breakpoints_d,coeff_fvn_d)
      write(2,44) (xp_d,ty_d,fvn_y_d)
   end do
  
   close(2)
  
  44      FORMAT(4(1X,1PE22.14))
  
  end program
  
  \end{verbatim}
  Results are plotted on figure \ref{akima}
  
  \begin{figure}
   \begin{center}
   \includegraphics[width=0.9\textwidth]{akima.pdf}
   % akima.pdf: 504x720 pixel, 72dpi, 17.78x25.40 cm, bb=0 0 504 720
   \caption{Akima Spline Interpolation}
   \label{akima}
  \end{center}
  
  \end{figure}
  
  
  
  \section{Least square polynomial}
  fvn provide a function to find a least square polynomial of a given degree, for real in single or double precision. It is performed using Lapack subroutine sgelss (dgelss), which solve this problem using singular value decomposition.
  
  \begin{verbatim}

  call fvn_lspoly(np,x,y,deg,coeff,status)

  \end{verbatim}
  \begin{itemize}
   \item np (in) is an integer equal to the number of points
   \item x(np) (in),y(np) (in) are the known coordinates
   \item deg (in) is an integer equal to the degree of the desired polynomial, it must be lower than np.
   \item coeff(deg+1) (out) on output contains the polynomial coefficients
   \item status (out) is an integer containing 0 if a problem occured.
  \end{itemize}
  
  \subsection*{Example}
  Here's a simple example : we've got 13 measurement points and we want to find the least square degree 3 polynomial for these points :
  \begin{verbatim}
   program lsp
   use fvn
   implicit none
  
   integer,parameter :: npoints=13,deg=3
   integer :: status,i
   real(kind=8) :: xm(npoints),ym(npoints),xstep,xc,yc
   real(kind=8) :: coeff(deg+1)
  
   xm = (/ -3.8,-2.7,-2.2,-1.9,-1.1,-0.7,0.5,1.7,2.,2.8,3.2,3.8,4. /)
   ym = (/ -3.1,-2.,-0.9,0.8,1.8,0.4,2.1,1.8,3.2,2.8,3.9,5.2,7.5  /)
  
   open(2,file='fvn_lsp_double_mesure.dat')
   open(3,file='fvn_lsp_double_poly.dat')
  
   do i=1,npoints
      write(2,44) xm(i),ym(i)
   end do
   close(2)
  
  
   call fvn_d_lspoly(npoints,xm,ym,deg,coeff,status)
  
   xstep=(xm(npoints)-xm(1))/1000.
   do i=1,1000
      xc=xm(1)+(i-1)*xstep
      yc=poly(xc,coeff)
      write(3,44) xc,yc
   end do
   close(3)
  
  44      FORMAT(4(1X,1PE22.14))
  
  contains
  function poly(x,coeff)
      implicit none
      real(8) :: x
      real(8) :: coeff(deg+1)
      real(8) :: poly
      integer :: i
  
      poly=0.
  
      do i=1,deg+1
          poly=poly+coeff(i)*x**(i-1)
      end do
  
  end function
  end program
  \end{verbatim}
  The results are plotted on figure \ref{lsp} .
  
  \begin{figure}
   \begin{center}
   \includegraphics[width=0.9\textwidth]{lsp.pdf}
   \caption{Least Square Polynomial}
   \label{lsp}
   \end{center}
  \end{figure}
  
  
  
  \section{Zero finding}
  fvn provide a routine for finding zeros of a complex function using Muller algorithm (only for double complex type). It is based on a version provided on the web by Hans D Mittelmann \url{http://plato.asu.edu/ftp/other\_software/muller.f}.
  
  \begin{verbatim}

   call fvn_muller(f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)

  \end{verbatim}
  \begin{itemize}
   \item f (in) is the complex function (kind=8) for which we search zeros
   \item eps (in) is a real(8) corresponding to the first stopping criterion : let fp(z)=f(z)/p where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) and z(1),...,z(k-1) are previously found roots.  if ((cdabs(f(z)).le.eps) .and. (cdabs(fp(z)).le.eps)), then z is accepted as a root.
   \item eps1 (in) is a real(8) corresponding to the second stopping criterion : a root is accepted if two successive approximations to a given root agree within eps1. Note that if either or both of the stopping criteria are fulfilled, the root is accepted.
   \item kn (in) is an integer equal to the number of known roots, which must be stored in x(1),...,x(kn), prior to entry in the subroutine.
   \item nguess (in) is the number of initial guesses provided. These guesses must be stored in x(kn+1),...,x(kn+nguess). nguess must be set equal to zero if no guesses are provided.
   \item n (in) is an integer equal to the number of new roots to be found.
   \item x (inout) is a complex(8) vector of length kn+n. x(1),...,x(kn) on input must contain any known roots.  x(kn+1),..., x(kn+n) on input may, on user option, contain initial guesses for the n new roots which are to be computed. If the user does not provide an initial guess, zero is used. On output, x(kn+1),...,x(kn+n) contain the approximate roots found by the subroutine.
   \item itmax (in) is an integer equal to the maximum allowable number of iterations per root.
   \item infer (out) is an integer vector of size kn+n. On output infer(j) contains the number of iterations used in finding the j-th root when convergence was achieved.  If convergence was not obtained in itmax iterations, infer(j) will be greater than itmax
   \item ier (out) is an integer used as an error parameter. ier = 33 indicates failure to converge within itmax iterations for at least one of the (n) new roots.
  \end{itemize}
  This subroutine always returns the last approximation for root j in x(j). if the convergence criterion is satisfied, then infer(j) is less than or equal to itmax. if the convergence criterion is not satisified, then infer(j) is set to either itmax+1 or itmax+k, with k greater than 1. infer(j) = itmax+1 indicates that muller did not obtain convergence in the allowed number of iterations. in this case, the user may wish to set itmax to a larger value. infer(j) = itmax+k means that convergence was obtained (on iteration k) for the deflated function fp(z) = f(z)/((z-z(1)...(z-z(j-1))) but failed for f(z). in this case, better initial guesses might help or, it might be necessary to relax the convergence criterion.
  
  \subsection*{Example}
  Example to find the ten roots of $x^{10}-1$
  \begin{verbatim}
   program muller
   use fvn
   implicit none
  
   integer :: i,info
   complex(8),dimension(10) :: roots
   integer,dimension(10) :: infer
   complex(8), external :: f
  
   call fvn_z_muller(f,1.d-12,1.d-10,0,0,10,roots,200,infer,info)
  
   write(*,*) "Error code :",info
   do i=1,10
      write(*,*) roots(i),infer(i)
   enddo
   end program
  
   function f(x)
      complex(8) :: x,f
      f=x**10-1
   end function
  
  \end{verbatim}

  
  \section{Numerical integration}
  Using an integrated slightly modified version of quadpack \url{http://www.netlib.org/quadpack}, fvn provide adaptative numerical integration (Gauss Kronrod) of real functions of 1 and 2 variables. fvn also provide a function to calculate Gauss-Legendre abscissas and weight, and a simple non adaptative integration subroutine. All routines exists only in fvn for double precision real.
  
  \subsection{Gauss Legendre Abscissas and Weigth}
  This subroutine was inspired by Numerical Recipes routine gauleg.
  \begin{verbatim}

  call fvn_gauss_legendre(n,qx,qw)

  \end{verbatim}
  \begin{itemize}
   \item n (in) is an integer equal to the number of Gauss Legendre points
   \item qx (out) is a real(8) vector of length n containing the abscissas.
   \item qw (out) is a real(8) vector of length n containing the weigths.
  \end{itemize}
  This subroutine computes n Gauss-Legendre abscissas and weigths
  
  \subsection{Gauss Legendre Numerical Integration}
  \begin{verbatim}

  call fvn_gl_integ(f,a,b,n,res)

  \end{verbatim}
  \begin{itemize}
   \item f (in) is a real(8) function to integrate
   \item a (in) and b (in) are real(8) respectively lower and higher bound of integration
   \item n (in) is an integer equal to the number of Gauss Legendre points to use
   \item res (out) is a real(8) containing the result
  \end{itemize}
  This function is a simple Gauss Legendre integration subroutine, which evaluate the integral of function f as in equation \ref{intsple} using n Gauss-Legendre pairs.
  
  \subsection{Gauss Kronrod Adaptative Integration}
  This kind of numerical integration is an iterative procedure which try to achieve a given precision.
  \subsubsection{Numerical integration of a one variable function}
  \begin{verbatim}

  call fvn_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

  \end{verbatim}
  This routine evaluate the integral of function f as in equation \ref{intsple}
  \begin{itemize}
   \item f (in) is an external real(8) function of one variable
   \item a (in) and b (in) are real(8) respectively lower an higher bound of integration
   \item epsabs (in) and epsrel (in) are real(8) respectively desired absolute and relative error
   \item key (in) is an integer between 1 and 6 correspondind to the Gauss-Kronrod rule to use :
      \begin{itemize}
          \item 1 : 7 - 15 points
          \item 2 : 10 - 21 points
          \item 3 : 15 - 31 points
          \item 4 : 20 - 41 points
          \item 5 : 25 - 51 points
          \item 6 : 30 - 61 points
      \end{itemize}
   \item res (out) is a real(8) containing the estimation of the integration.
   \item abserr (out) is a real(8) equal to the estimated absolute error
   \item ier (out) is an integer used as an error flag
      \begin{itemize}
          \item 0 : no error
          \item 1 : maximum number of subdivisions allowed has been achieved. one can allow more subdivisions by increasing the value of limit (and taking the according dimension adjustments into account). however, if this yield no improvement it is advised to analyze the integrand in order to determine the integration difficulaties. If the position of a local difficulty can be determined (i.e.singularity, discontinuity within the interval) one will probably gain from splitting up the interval at this point and calling the integrator on the subranges. If possible, an appropriate special-purpose integrator should be used which is designed for handling the type of difficulty involved.
          \item 2 : the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved.
          \item 3 : extremely bad integrand behaviour occurs at some points of the integration interval.
          \item 6 : the input is invalid, because (epsabs.le.0 and epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) or limit.lt.1 or lenw.lt.limit*4. result, abserr, neval, last are set to zero. Except when lenw is invalid, iwork(1), work(limit*2+1) and work(limit*3+1) are set to zero, work(1) is set to a and work(limit+1) to b.
      \end{itemize}

   \item limit (in) is an optional integer equal to maximum number of subintervals in the partition of the given integration interval (a,b). If the parameter is not present a default value of 500 will be used.

  \end{itemize}
  
  \begin{equation}
   \int_a^b f(x)~dx
   \label{intsple}
  \end{equation}
  
  
  
  
  \subsubsection{Numerical integration of a two variable function}
  \begin{verbatim}

  call fvn_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)

  \end{verbatim}
  This function evaluate the integral of a function f(x,y) as defined in equation \ref{intdble}. The parameters of same name as in the previous paragraph have exactly the same function and behaviour thus only what differs is decribed here
  \begin{itemize}
   \item a (in) and b (in) are real(8) corresponding respectively to lower and higher bound of integration for the x variable.
   \item g(x) (in) and h(x) (in) are external functions describing the lower and higher bound of integration for the y variable as a function of x.
  \end{itemize}
  
  \begin{equation}
   \int_a^b \int_{g(x)}^{h(x)} f(x,y)~dy~dx
   \label{intdble}
  \end{equation}
  
  \subsubsection*{Example}
  \begin{verbatim}
  program integ
   use fvn
   implicit none
  
   real(8), external :: f1,f2,g,h
   real(8) :: a,b,epsabs,epsrel,abserr,res
   integer :: key,ier
  
   a=0.
   b=1.
   epsabs=1d-8
   epsrel=1d-8
   key=2
   call fvn_d_integ_1_gk(f1,a,b,epsabs,epsrel,key,res,abserr,ier,500)
   write(*,*) "Integration of x*x between 0 and 1 : "
   write(*,*) res
  
   call fvn_d_integ_2_gk(f2,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,500)
   write(*,*) "Integration of x*y between 0 and 1 on both x and y : "
   write(*,*) res
   
  
  end program
  
  function f1(x)
   implicit none
      real(8) :: x,f1
      f1=x*x
  end function
  
  function f2(x,y)
   implicit none
      real(8) :: x,y,f2
      f2=x*y
  end function
  
  function g(x)
   implicit none
      real(8) :: x,g
      g=0.
  end function
  
  function h(x)
   implicit none
      real(8) :: x,h
      h=1.
  end function
  \end{verbatim}

  \section{Special functions}

  Specials functions are available in fvn by using an implementation of fnlib \url{http://www.netlib.org/fn}. This can be used separatly from the rest of fvn by using the module \verb'fvn_fnlib' and linking the library \verb'libfvn_fnlib.a' . The module provides a generic interfaces to all the routines. Specific names of the routines are given in the description. The double complex versions of the routines are not present in the web version of fnlib, so these have been added, but not intensely tested.
  
  \paragraph{Important Note}
  Due to the addition of fnlib to fvn, some functions that were in fvn and are redondant will be removed from fvn, so update your code now and replace them with the fnlib version. These are listed here after :
  \begin{itemize}
   \item \verb'fvn_z_acos' replaced by \verb'acos'
   \item \verb'fvn_z_asin' replaced by \verb'asin'
   \item \verb'fvn_d_asinh' replaced by \verb'asinh'
   \item \verb'fvn_d_acosh' replaced by \verb'acosh'
   \item \verb'fvn_s_csevl' replaced by \verb'csevl'
   \item \verb'fvn_d_csevl' replaced by \verb'csevl'
   \item \verb'fvn_d_factorial' replaced by \verb'fac'
   \item \verb'fvn_d_lngamma' replaced by \verb'alngam'
  \end{itemize}
  
  
  \subsection{Elementary functions}
  \subsubsection{carg}
  \begin{verbatim}
   carg(z)
  \end{verbatim}
  \begin{itemize}
   \item z (in) is a complex
  \end{itemize}
  This function evaluates the argument of the complex z. That is $\theta$ for $z=\rho e^{i\theta}$.
  
  Specific interfaces : \verb'carg,zarg'
  
  
  \subsubsection{cbrt}
  \begin{verbatim}
   cbrt(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the cubic root of the argument x.
  
  Specific interfaces : \verb'cbrt,dcbrt,ccbrt,zcbrt'
  
  
  \subsubsection{exprl}
  \begin{verbatim}
   exprl(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates ${e^x-1}\over x$.
  
  Specific interfaces : \verb'exprel,dexprl,cexprl,zexprl'
  
  \subsubsection{log10}
  \begin{verbatim}
   log10(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function is an extension of the intrinsic function log10 to complex arguments.
  
  Specific interfaces : \verb'clog10,zlog10'
  
  
  \subsubsection{alnrel}
  \begin{verbatim}
   alnrel(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates $ln(1+x)$.
  
  Specific interfaces : \verb'alnrel,dlnrel,clnrel,zlnrel'
  
  
  \subsection{Trigonometry}
  \subsubsection{tan}
  \begin{verbatim}
   tan(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the tangent of the argument. It is an extension of the intrinsic function tan to complex arguments.
  
  Specific interfaces : \verb'ctan,ztan'
  
  
  \subsubsection{cot}
  \begin{verbatim}
   cot(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluate the cotangent of the argument.
  
  Specific interfaces : \verb'cot,dcot,ccot,zcot'
  
  \subsubsection{sindg}
  \begin{verbatim}
   sindg(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluate the sinus of the argument expressed in degrees.
  
  Specific interfaces : \verb'sindg,dsindg'
  
  
  \subsubsection{cosdg}
  \begin{verbatim}
   cosdg(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluate the cosinus of the argument expressed in degrees.
  
  Specific interfaces : \verb'cosdg,dcosdg'
  
  
  \subsubsection{asin}
  \begin{verbatim}
   asin(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the arc sine of the argument. It is an extension of the intrinsic function asin to complex arguments.
  
  Specific interfaces : \verb'casin,zasin'
  
  \subsubsection{acos}
  \begin{verbatim}
   acos(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the arc cosine of the argument. It is an extension of the intrinsic function acos to complex arguments.
  
  Specific interfaces : \verb'cacos,zacos'
  
  
  \subsubsection{atan}
  \begin{verbatim}
   atan(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the arc tangent of the argument. It is an extension of the intrinsic function atan to complex arguments.
  
  Specific interfaces : \verb'catan,zatan'
  
  \subsubsection{atan2}
  \begin{verbatim}
   atan2(x,y)
  \end{verbatim}
  \begin{itemize}
   \item x,y are real or complex
  \end{itemize}
  This function evaluates the arc tangent of $x \over y$. It is an extension of the intrinsic function atan2 to complex arguments.
  
  Specific interfaces : \verb'catan2,zatan2'
  
  \subsubsection{sinh}
  \begin{verbatim}
   sinh(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the hyperbolic sine of the argument. It is an extension of the intrinsic function sinh to complex arguments.
  
  Specific interfaces : \verb'csinh,zsinh'
  
  
  \subsubsection{cosh}
  \begin{verbatim}
   cosh(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the hyperbolic cosine of the argument. It is an extension of the intrinsic function cosh to complex arguments.
  
  Specific interfaces : \verb'ccosh,zcosh'
  
  \subsubsection{tanh}
  \begin{verbatim}
   tanh(x)
  \end{verbatim}
  This function evaluates the hyperbolic tangent of the argument. It is an extension of the intrinsic function tanh to complex arguments.
  
  Specific interfaces : \verb'ctanh,ztanh'
  
  \subsubsection{asinh}
  \begin{verbatim}
   asinh(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the arc hyperbolic sine of the argument.
  
  Specific interfaces : \verb'asinh,dasinh,casinh,zasinh'
  
  \subsubsection{acosh}
  \begin{verbatim}
   acosh(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the arc hyperbolic cosine of the argument.
  
  Specific interfaces : \verb'acosh,dacosh,cacosh,zacosh'
  
  \subsubsection{atanh}
  \begin{verbatim}
   atanh(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the arc hyperbolic tangent of the argument.
  
  Specific interfaces : \verb'atanh,datanh,catanh,zatanh'
  
  \subsection{Exponential Integral and related}
  \subsubsection{ei}
  \begin{verbatim}
   ei(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the exponential integral for argument greater then 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{ei} for $x 
  eq 0$.
  \begin{equation}
  \label{ei}
   ei(x)= - \int _{-x} ^\infty {e^{-t}\over t}dt
  \end{equation}
  
  Specific interfaces : \verb'ei,dei'
  
  
  \subsubsection{e1}
  \begin{verbatim}
   e1(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x 
  eq 0$.
  \begin{equation}
  \label{e1}
   e1(x)= \int _{x} ^\infty {e^{-t}\over t}dt
  \end{equation}
  
  Specific interfaces : \verb'e1,de1'
  
  \subsubsection{ali}
  \begin{verbatim}
   ali(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the logarithm integral. it is define by equation \ref{ali} for $x > 0$ and $x 
  eq 1$.
  \begin{equation}
   \label{ali}
   ali(x)= - \int _0 ^x {dt \over ln(x)}
  \end{equation}
  
  Specific interfaces : \verb'ali,dli'
  
  \subsubsection{si}
  \begin{verbatim}
   si(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the sine integral defined by equation \ref{si}.
  \begin{equation}
   \label{si}
   si(x)= \int _0 ^x {sin(t) \over t }dt
  \end{equation}
  
  Specific interfaces : \verb'si,dsi'
  
  
  \subsubsection{ci}
  \begin{verbatim}
   ci(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the cosine integral defined by equation \ref{ci} where $\gamma \approx 0.57721566$ represent Euler's constant.
  \begin{equation}
   \label{ci}
   ci(x)= \gamma + ln(x) + \int _0 ^x {{1-cos(t)} \over t} dt
  \end{equation}
  
  Specific interfaces : \verb'ci,dci'
  
  \subsubsection{cin}

  \begin{verbatim}

   cin(x)

  \end{verbatim}
  \begin{itemize}

   \item x is a real
  \end{itemize}
  This function evaluates the cosine integral alternate definition given by equation \ref{cin}.
  \begin{equation}
   \label{cin}
   cin(x)= \int _0 ^x {{1-cos(t)} \over t} dt
  \end{equation}
  
  Specific interface : \verb'cin,dcin'
  
  \subsubsection{shi}
  \begin{equation}
   shi(x)
  \end{equation}
  \begin{itemize}
   \item x is a real

  \end{itemize}

  This function evaluates the hyperbolic sine integral defined by equation \ref{shi}.
  \begin{equation}
   \label{shi}
   shi(x) = \int _0 ^x {sinh(t) \over t}dt
  \end{equation}
  
  Specific interfaces : \verb'shi,dshi'

  \subsubsection{chi}

  \begin{verbatim}

   chi(x)

  \end{verbatim}
  \begin{itemize}

   \item x is a real

  \end{itemize}

  This function evaluates the hyperbolic cosine integral defined by equation \ref{chi} where $\gamma \approx 0.57721566$ represent Euler's constant.
  \begin{equation}
   \label{chi}
   chi(x)= \gamma + ln(x) + \int _0 ^x {{cosh(t) -1} \over t}dt
  \end{equation}
  
  Specific interfaces : chi,dchi

  
  \subsubsection{cinh}

  \begin{verbatim}

   cinh(x)

  \end{verbatim}
  \begin{itemize}
   \item x is a real

  \end{itemize}
  This function evaluates the hyperbolic cosine integral alternate definition given by equation \ref{cinh}.
  \begin{equation}
   \label{cinh}
   cinh(x) = \int _0 ^x {{cosh(t) -1} \over t}dt
  \end{equation}
  
  Specific interfaces : cinh,dcinh
  
  
  \subsection{Gamma function and related}
  \subsubsection{fac}
  \begin{verbatim}
   fac(n)

   dfac(n)

  \end{verbatim}
  \begin{itemize}

   \item n is an integer
  \end{itemize}

  This function return $n!$ as a real(4) or real(8) for dfac. There's no generic interface for this one.

  
  Specific interfaces : \verb'fac,dfac'
  
  \subsubsection{binom}
  \begin{verbatim}
   binom(n,m)

   dbinom(n,m)

  \end{verbatim}
  \begin{itemize}
   \item n,m are integers
  \end{itemize}

  This function return the binomial coefficient defined by equation \ref{binom} with $n \geq m \geq 0$. binom returns a real(4), dbinom a real(8). There's no generic interface for this one.

  \begin{equation}
   \label{binom}
   binom(n,m) = C_n^m = {{n!} \over {m!(n-m)!}}
  \end{equation}
  
  Specific interfaces : \verb'binom,dbinom'
  
  
  \subsubsection{gamma}
  \begin{verbatim}
   gamma(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates $ \Gamma (x) $ defined by equation \ref{gamma}.
  \begin{equation}
   \label{gamma}
   \Gamma (x) = \int _0 ^{\infty} t^{x-1}e^{-t}dt
  \end{equation}
  Note that $n!=\Gamma (n+1)$.
  
  Specific interfaces :\verb'gamma,dgamma,cgamma,zgamm'
  
  \subsubsection{gamr}
  \begin{verbatim}
   gamr(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the reciprocal gamma function $gamr(x)= {1 \over \Gamma(x)}$
  
  
  \subsubsection{alngam}
  \begin{verbatim}
   alngam(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates $ln(|\Gamma(x)|)$
  
  Specific interfaces : \verb'alngam,dlngam,clngam,zlngam'
  
  
  \subsubsection{algams}
  \begin{verbatim}
   call algams(x,algam,sgngam)
  \end{verbatim}
  \begin{itemize}
   \item x (in) is a real
   \item algam (out) is a real
   \item sgngam (out) is a real
  \end{itemize}
  This subroutine evaluates the logarithm of the absolute value of gamma and the sign of gamma. 
  $algam=ln(|\Gamma(x)|)$ and $sgngam=1.0$ or $-1.0$ according to the sign of $\Gamma(x)$.
  
  Specific interfaces : \verb'algams,dlgams'
  
  \subsubsection{gami}
  \begin{verbatim}
   gami(a,x)
  \end{verbatim}
  \begin{itemize}
   \item x is a positive real
   \item a is a strictly positive real
  \end{itemize}
  This function evaluates the incomplete gamma function defined by equation \ref{gami}.
  \begin{equation}
   \label{gami}
   gami(a,x)=\gamma(a,x)=\int _0 ^x t^{a-1} e^{-t}dt
  \end{equation}
  
  Specific interfaces : \verb'gami,dgami'
  
  \subsubsection{gamic}
  \begin{verbatim}
   gamic(a,x)
  \end{verbatim}
  \begin{itemize}
   \item x is a positive real
   \item a is a real
  \end{itemize}
  This function evaluates the complementary incomplete gamma function defined by equation \ref{gamic}.
  \begin{equation}
   \label{gamic}
   gamic(a,x)=\Gamma(a,x)=\int _x ^\infty t^{a-1} e^{-t}dt
  \end{equation}
  
  Specific interfaces : \verb'gamic,dgamic'
  
  \subsubsection{gamit}
  \begin{verbatim}
   gamit(a,x)
  \end{verbatim}
  \begin{itemize}
   \item x is a positive real
   \item a is a real
  \end{itemize}
  This function evaluates the Tricomi's incomplete gamma function defined by equation \ref{gamit}.
  \begin{equation}
   \label{gamit}
   gamit(a,x)=\gamma^* (a,x)= {{x^{-a}\gamma(a,x)}\over \Gamma(a)}
  \end{equation}
  
  Specific interfaces : \verb'gamit,dgamit'
  
  
  \subsubsection{psi}
  \begin{verbatim}
   psi(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real or complex
  \end{itemize}
  This function evaluates the psi function which is the logarithm derivative of the gamma function as defined in equation \ref{psi}.
  \begin{equation}
   \label{psi}
    psi(x)= \psi(x) = {d\over dx} ln(\Gamma(x))
  \end{equation}
  x must not be zero or a negative integer.
  
  Specific interfaces : \verb'psi,dpsi,cpsi,zpsi'
  
  
  \subsubsection{poch}
  \begin{verbatim}
   poch(a,x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
   \item a is a real
  \end{itemize}
  This function evaluates a generalization of Pochhammer's symbol.
  
  Pochhammer's symbol for n a positive integer is given by equation \ref{poch_int}
  \begin{equation}
   \label{poch_int}
  (a)_n = a(a-1)(a-2)...(a-n+1)
  \end{equation}
  
  The generalization of Pochhammer's symbol is given by equation \ref{poch}
  \begin{equation}
   \label{poch}
   poch(a,x)= (a)_x = {\Gamma(a+x) \over \Gamma(a)}
  \end{equation}
  
  Specific interfaces : \verb'poch,dpoch'
  
  
  \subsubsection{poch1}
  \begin{verbatim}
   poch1(a,x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
   \item a is a real
  \end{itemize}
  This function is defined by equation \ref{poch1}. It is usefull for certains situations, especially when x is small.
  
  \begin{equation}
   \label{poch1}
   poch1(a,x)={{(a)_x-1} \over x}
  \end{equation}
  
  Specific interfaces : \verb'poch1,dpoch1'
  
  \subsubsection{beta}
  \begin{verbatim}
   beta(a,b)
  \end{verbatim}
  \begin{itemize}
   \item a,b are real positive or complex
  \end{itemize}
  This function evaluates $\beta$ function defined by equation \ref{beta}.
  \begin{equation}
   \label{beta}
  beta(a,b)=\beta(a,b)={  {\Gamma(a) \Gamma(b)} \over \Gamma(a+b) }
  \end{equation}
  
  Specific interfaces : \verb'beta,dbeta,cbeta,zbeta'
  
  
  \subsubsection{albeta}
  \begin{verbatim}
   albeta(a,b)
  \end{verbatim}
  \begin{itemize}
   \item a,b are real positive or complex
  \end{itemize}
  This function evaluates the natural logarithm of beta function : $ln(\beta(a,b))$
  
  Specific interfaces : \verb'albeta,dlbeta,clbeta,zlbeta'
  
  \subsubsection{betai}
  \begin{verbatim}
   betai(x,pin,qin)
  \end{verbatim}
  \begin{itemize}
   \item x is a real in [0,1]
   \item pin and qin are strictly positive real
  \end{itemize}
  This function evaluates the incomplete beta function ratio, that is the probability that a random variable from a beta distribution having parameters pin and qin will be less than or equal to x. It is defined by equation \ref{betai}.
  
  \begin{equation}
   \label{betai}
  betai(x,pin,qin)=I_x(pin,qin)={1 \over \beta(pin,qin)} \int _0 ^x t^{pin-1}(1-t)^{qin-1}dt
  \end{equation}
  
  Specific interfaces : \verb'betai,dbetai'
  
  \subsection{Error function and related}
  \subsubsection{erf}
  \begin{verbatim}
   erf(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the error function defined by equation \ref{erf}.
  \begin{equation}
   \label{erf}
   erf(x)={2\over \sqrt{ \pi}} \int _0 ^x e^{-t^2}dt
  \end{equation}
  
  Specific interfaces : \verb'erf,derf'
  
  \subsubsection{erfc}
  \begin{verbatim}
   erfc(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the complimentary error function defined by equation \ref{erfc}.
  \begin{equation}
   \label{erfc}
   erfc(x)={2\over \sqrt{ \pi}} \int _x ^\infty e^{-t^2}dt
  \end{equation}
  
  Specific interfaces : \verb'erfc,derfc'
  
  
  \subsubsection{daws}
  \begin{verbatim}
   daws(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates Dawson's function defined by equation \ref{daws}.
  \begin{equation}
   \label{daws}
  daws(x)=e^{-x^2} \int _0 ^x e^{t^2}dt
  \end{equation}
  
  Specific interfaces : \verb'daws,ddaws'
  
  \subsection{Bessel functions and related}
  \subsubsection{bsj0}
  \begin{verbatim}
   bsj0(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates Bessel function of the first kind of order 0 defined by equation \ref{bsj0}.
  \begin{equation}
   \label{bsj0}
   bsj0(x)=J_0(x)= {1 \over \pi} \int _0 ^\pi cos(x sin(\theta)) d\theta
  \end{equation}
  
  Specific interfaces : \verb'besj0,dbesj0'
  
  \subsubsection{bsj1}
  \begin{verbatim}
   bsj1(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates Bessel function of the first kind of order 1 defined by equation \ref{bsj1}.
  \begin{equation}
   \label{bsj1}
  bsj1(x)=J_1(x)={1 \over \pi} \int _0 ^\pi cos(x sin(\theta)-\theta) d\theta
  \end{equation}
  
  Specific interfaces : \verb'besj1,dbesj1'
  
  \subsubsection{bsy0}
  \begin{verbatim}
   bsy0(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the Bessel function of the second kind of order 0 defined by equation \ref{bsy0}
  \begin{equation}
   \label{bsy0}
   bsy0(x)=Y_0(x)={1 \over \pi} \int _0 ^\pi sin(x sin(\theta))d\theta -{2 \over \pi} \int _0 ^\infty e^{-x sinh(t)}dt
  \end{equation}
  
  Specific interfaces : \verb'besy0,dbesy0'
  
  \subsubsection{bsy1}
  \begin{verbatim}
   bsy1(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the Bessel function of the second kind of order 1 defined by equation \ref{bsy1}.
  \begin{equation}
   \label{bsy1}
   bsy1(x)=Y_1(x)=-{1 \over \pi} \int _0 ^\pi sin (\theta - x sin(\theta)) d\theta 
  - {1 \over \pi} \int _0 ^\infty (e^t -e^{-t})e^{-x sinh(t)}dt
  \end{equation}
  
  Specific interfaces : \verb'besy1,dbesy1'
  
  \subsubsection{bsi0}
  \begin{verbatim}
   bsi0(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the Bessel function of the third kind of order 0 defined by equation \ref{bsi0}.
  \begin{equation}
   \label{bsi0}
  bsi0(x)=I_0(x)={1 \over \pi} \int _0 ^\pi cosh(x cos(\theta))d\theta
  \end{equation}
  
  Specific interfaces : \verb'besi0,dbesi0'
  
  \subsubsection{bsi1}
  \begin{verbatim}
   bsi1(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the Bessel function of the third kind of order 1 defined by equation \ref{bsi1}.
  \begin{equation}
   \label{bsi1}
  bsi1(x)=I_1(x)={1 \over \pi} \int _0 ^\pi e^{x cos(\theta)} cos(\theta)d\theta
  \end{equation}
  
  Specific interfaces : \verb'besi1,dbesi1'
  
  \subsubsection{bsk0}
  \begin{verbatim}
   bsk0(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a strictly positive real
  \end{itemize}
  This function evaluates the modified Bessel function of the second kind of order 0 defined by equation \ref{bsk0}
  \begin{equation}
   \label{bsk0}
  bsk0(x)=K_0(x)=\int _0 ^\infty cos(x sinh(t))dt
  \end{equation}
  
  Specific interfaces : \verb'besk0,dbesk0'
  
  \subsubsection{bsk1}
  \begin{verbatim}
   bsk1(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a strictly positive real
  \end{itemize}
  This function evaluates the modified Bessel function of the second kind of order 1 defined by equation \ref{bsk1}
  \begin{equation}
   \label{bsk1}
  bsk1(x)=K_1(x)=\int _0 ^\infty sin(x sinh(t))sinh(t)dt
  \end{equation}
  
  Specific interfaces : \verb'besk1,dbesk1'
  
  \subsubsection{bsi0e}
  \begin{verbatim}
   bsi0e(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates $e^{-|x|}I_0(x)$
  
  Specific interfaces : \verb'besi0e,dbsi0e'
  
  \subsubsection{bsi1e}
  \begin{verbatim}
   bsi1e(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates $e^{-|x|}I_1(x)$
  
  Specific interfaces : \verb'besi1e,dbsi1e'
  
  \subsubsection{bsk0e}
  \begin{verbatim}
   bsk0e(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a strictly positive real
  \end{itemize}
  This function evaluates $e^x K_0(x)$
  
  Specific interfaces : \verb'besk0e,dbsk0e'
  
  \subsubsection{bsk1e}
  \begin{verbatim}
   bsk1e(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a strictly positive real
  \end{itemize}
  This function evaluates $e^x K_1(x)$
  
  Specific interfaces : \verb'besk1e,dbsk1e'
  
  \subsubsection{bsks}
  \begin{verbatim}
   call bsks(xnu,x,nin,bk)
  \end{verbatim}
  \begin{itemize}
   \item xnu (in) is a real with $|xnu|<1$. It's the fractional order
   \item x (in) is a real. The value for which the sequence of Bessel functions is to be evaluated.
   \item nin (in) is an integer.
   \item bk (out) is a real vector of length abs(nin), containing the values of the function.
  \end{itemize}
  This subroutine evaluates a sequence of modified Bessel function of the second kind of fractional order.
  
  If nin is positive, on completion $bk(1)=K_
  u(x)$,$bk(2)=K_{
  u+1}(x)$,...,$bk(nin)=K_{
  u+nin-1}(x)$. If nin is negative, on completion $bk(1)=K_
  u(x)$,$bk(2)=K_{
  u-1}(x)$,...,$bk(|nin|)=K_{
  u+nin+1}(x)$.
  
  Specific interfaces : \verb'besks,dbesks'
  
  \subsubsection{bskes}
  \begin{verbatim}
   call bskes(xnu,x,nin,bke)
  \end{verbatim}
  \begin{itemize}
   \item xnu (in) is a real with $|xnu|<1$. It's the fractional order
   \item x (in) is a real. The value for which the sequence of exponentialy scaled Bessel functions is to be evaluated.
   \item nin (in) is an integer. Number of elements in the sequence.
   \item bke (out) is a real vector of length abs(nin), containing the values of the function.
  \end{itemize}
  This subroutine evaluates a sequence of exponentially scaled modified Bessel function of the second kind of fractional order.
  
  If nin is positive, on completion $bk(1)=e^x K_
  u(x)$,$bk(2)=e^x K_{
  u+1}(x)$,...,$bk(nin)=e^x K_{
  u+nin-1}(x)$. If nin is negative, on completion $bk(1)=e^x K_
  u(x)$,$bk(2)=e^x K_{
  u-1}(x)$,...,$bk(|nin|)=e^x K_{
  u+nin+1}(x)$.
  
  Specific interfaces : \verb'beskes,dbskes'
  
  \subsection{Airy function and related}
  \subsubsection{ai}
  \begin{verbatim}
   ai(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the airy function defined by equation \ref{ai}
  \begin{equation}
   \label{ai}
  Ai(x)={1 \over \pi} \int _0 ^\infty cos(xt+ {1 \over 3}t^3)dt
  \end{equation}
  
  Specific interfaces : \verb'ai,dai'
  
  \subsubsection{bi}
  \begin{verbatim}
   bi(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the Airy function of the second kind defined by equation \ref{bi}
  \begin{equation}
   \label{bi}
  Bi(x)={1 \over \pi} \int _0 ^\infty e^{xt- {1 \over 3}t^3}dt +
  {1 \over \pi} \int _0 ^\infty sin(xt+ {1 \over 3}t^3)dt
  \end{equation}
  
  Specific interfaces : \verb'bi,dbi'
  
  
  \subsubsection{aid}
  \begin{verbatim}
   aid(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the derivative of the Airy function, $aid(x)={d \over {dx}}Ai(x)$.
  
  Specific interface : \verb'aid,daid'
  
  \subsubsection{bid}
  \begin{verbatim}
   bid(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the derivative of the Airy function of the second kind, $bid(x)={d \over {dx}}Bi(x)$.
  
  Specific interfaces : \verb'bid,dbid'
  
  
  \subsubsection{aie}
  \begin{verbatim}
   aie(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  
  This function evaluates the exponentially scaled Airy function defined in equation \ref{aie}.
  \begin{equation}
   \label{aie}
  aie(x)=Ai(x) \textrm{~if~}x\leq0 \qquad \qquad aie(x)=e^{{2\over3}x^{3\over2}}Ai(x)\rm{~if~}x>0
  \end{equation}

  %\begin{equation}
  %  \left{\begin{array}{ll}
  %  Ai(x) & \textrm{~if~}x\leq0 \\
  %  e^{{2\over3}x^{3\over2}}Ai(x) & \textrm{~if~}x>0
  %\end{array}\right. 
  %\label{aie}
  %\end{equation}

  
  Specific interfaces : \verb'aie,daie'
  
  \subsubsection{bie}
  
  \begin{verbatim}
   bie(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the exponentially scaled Airy function of the second kind defined in equation \ref{bie}.
  \begin{equation}
   \label{bie}
  bie(x)=Bi(x)\textrm{~if~}x\leq0 \qquad \qquad bie(x)=e^{-{2\over3}x^{3\over2}}Bi(x)\rm{~if~}x>0
  \end{equation}
  
  Specific interfaces : \verb'bie,dbie'
  
  \subsubsection{aide}
  \begin{verbatim}
   aide(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the exponentially scaled derivative of the Airy function as defined in equation \ref{aide}.
  \begin{equation}
   \label{aide}
  aie(x)=Ai^\prime(x) \textrm{~if~}x\leq0 \qquad \qquad aie(x)=e^{{2\over3}x^{3\over2}}Ai^\prime(x)\rm{~if~}x>0
  \end{equation}
  
  Specific interfaces : \verb'aide,daide'
  
  
  \subsubsection{bide}
  \begin{verbatim}
   bide(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates the exponentially scaled derivative of the Airy function of the second kind as defined in equation \ref{bide}.
  \begin{equation}
   \label{bide}
  bie(x)=Bi^\prime(x)\textrm{~if~}x\leq0 \qquad \qquad bie(x)=e^{-{2\over3}x^{3\over2}}Bi^\prime(x)\rm{~if~}x>0
  \end{equation}
  
  Specific interfaces : \verb'bide,dbide'
  
  
  \subsection{Miscellanous functions}
  \subsubsection{spenc}
  \begin{verbatim}
   spenc(x)
  \end{verbatim}
  \begin{itemize}
   \item x is a real
  \end{itemize}
  This function evaluates Spence function defined in equation \ref{spenc}.
  \begin{equation}
   \label{spenc}
  spenc(x)=s(x)=- \int_0^x {{ln(|1-t|)}\over t}dt
  \end{equation}
  
  Specific interfaces : \verb'spenc,dspenc'
  
  
  \subsubsection{inits}
  \begin{verbatim}
   inits(os,nos,eta)
  \end{verbatim}
  \begin{itemize}
   \item os is a real vector of length nos, containing the coefficients in an orthogonal series.
   \item nos is an integer
   \item eta is a real (Warning eta is a real(4) even with the double precision version) representing the requested accuracy.
  \end{itemize}
  This function initialize the orthogonal series so that inits is the number of terms needed to insure the error is no larger than eta.
  
  Specific interfaces : \verb'inits,initds'
  
  
  \subsubsection{csevl}
  \begin{verbatim}
   csevl(x,cs,n)
  \end{verbatim}
  \begin{itemize}
   \item x is a real in [-1,1]
   \item cs is a real vector of length n containing the coefficients of the Chebyshev serie.
   \item n is an integer
  \end{itemize}
  This function evaluates the Chebyshev series whose coefficients are stored in cs.
  
  Specific interfaces : \verb'csevl,dcsevl'

  
  
  \end{document}