c      implicit real*8 (a-h,o-z)
      parameter (ndim=2)
      real p(ndim+1,ndim),y(ndim+1),pb(ndim),pom(ndim),yb,ftol,temptr
      real f
      external f
      temptr=0.1d0
      iter=100
      ftol=1.0d-4
      yb=1.0d6
      do i=1,ndim
        do j=1,ndim+1
          p(j,i)=1.0d1
        enddo
      enddo

      p(1,1)=-100.0d0
      p(1,2)=-50.0d0
      p(2,1)=-80.0d0
      p(2,2)=6.0d1
      p(3,1)=4.0d0
      p(3,2)=20.0d0
      do i=1,ndim+1
        do j=1,ndim
          pom(j)=p(i,j)
        enddo
        y(i)=f(pom)
      enddo

      call amebsa(p,y,ndim+1,ndim,ndim,pb,yb,ftol,f,iter,temptr)

      print *,"pb",pb
      print *,"yb",yb
      print *,"iter",iter

      stop
      end

      real function f(x)
      real x(2)
      f=0.5*((x(1)-1)**2+(x(2)+1)**2)+10*(1-cos(x(1)-1))*(1-cos(x(2)+1))
      return
      end

CcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC
      REAL FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *     NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
C     "Minimal" random number generator of Park and Miller with Bays-Durham shuffle and
C     added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of
C     the endpoint values). Call with idum a negative integer to initialize; thereafter, do not
C     alter idum between successive deviates in a sequence. RNMX should approximate the largest
C     floating value that is less than 1. Recommended for seqences smaller than 100e6 (i.e. 5%
C     of ran1's period.
C     (from Press etal, Numerical Recipes in Fortran 77)
CcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then  ! initialize
         idum=max(-idum,1)            ! load the shuffle table (after 8 warmups)
         do j=NTAB+8,1,-1
            k=idum/IQ
            idum=IA*(idum-k*IQ)-IR*k
            if (idum.lt.0) idum=idum+IM
            if (j.le.NTAB) iv(j)=idum
         enddo
         iy=iv(1)
      endif
      k=idum/IQ                       ! start here when not initializing
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END