anchor.diff:        return 
         END
   
  -      DOUBLE PRECISION FUNCTION factrl(n) 
  +      REAL FUNCTION factrl(n) 
   C Returns the value n! as a floating-point number. 
   
         INTEGER n 
   
         INTEGER j,ntop 
   C Table to be filled in only as required. 
  -      DOUBLE PRECISION a(33),gammln 
  +      REAL a(33),gammln 
         SAVE ntop,a 
   C Table initialized with 0! only. 
         DATA ntop,a(1)/0,1./
  
positive.code:       
  c matrix multiplication
  
  cf2py double precision dimension(m,q),intent(out) :: prod
  cf2py double precision dimension(m,n),intent(in) :: mat1
  cf2py double precision dimension(p,q),intent(in) :: mat2
  cf2py integer intent(hide),depend(mat1) :: m=len(mat1),n=shape(mat1,1)
  cf2py integer intent(hide),depend(mat2) :: p=len(mat2),q=shape(mat2,1)
  
  
        INTEGER i,j,k,m,n,p,q
        DOUBLE PRECISION mat1(m,n), mat2(p,q), prod(m,q)
        DOUBLE PRECISION sum
        
        if (n.eq.p) then
          do 30 i = 1,m
  
negative.code: C USES factln Returns the binomial coefficient as a 
  C floating point number.
        INTEGER k,n
        REAL bico
  C The nearest-integer function cleans up roundoff error
  C for smaller values of n and k. 
        bico=nint(exp(factln(n)-factln(k)-factln(n-k))) 
  