380 lines
9.4 KiB
Fortran
380 lines
9.4 KiB
Fortran
! Max Tegmark 171119, 190128-31, 190506, May 2020
|
|
|
|
! Binary:
|
|
! +: add
|
|
! *: multiply
|
|
! -: subtract
|
|
! /: divide (Put "D" instead of "/" in file, since f77 can't load backslash
|
|
! Unary:
|
|
! >: increment (x -> x+1)
|
|
! <: decrement (x -> x-1)
|
|
! ~: negate (x-> -x)
|
|
! \: invert (x->1/x) (Put "I" instead of "\" in file, since f77 can't load backslash
|
|
! L: logaritm (x-> ln(x)
|
|
! E: exponentiate (x->exp(x))
|
|
! S: sin: (x->sin(x))
|
|
! C: cos: (x->cos(x))
|
|
! A: abs: (x->abs(x))
|
|
! N: arcsin (x->arcsin(x))
|
|
! T: arctan (x->arctan(x))
|
|
! R: sqrt (x->sqrt(x))
|
|
! O: double (x->2*x); note that this is the letter "O", not zero
|
|
! J: double+1 (x->2*x+1)
|
|
! nonary:
|
|
! 0
|
|
! 1
|
|
! P: pi
|
|
real*8 function f(n,arities,ops,x) ! n=number of ops, x=arg vector
|
|
implicit none
|
|
integer nmax, n, i, j, arities(n), arity, lnblnk
|
|
character*60 ops
|
|
parameter(nmax=100)
|
|
real*8 x(nmax), y, stack(nmax)
|
|
character op
|
|
!write(*,*) 'Evaluating function with ops = ',ops(1:n)
|
|
!write(*,'(3f10.5,99i3)') (x(i),i=1,3), (arities(i),i=1,n)
|
|
j = 0 ! Number of numbers on the stack
|
|
do i=1,n
|
|
arity = arities(i)
|
|
op = ops(i:i)
|
|
if (arity.eq.0) then ! This is a nonary function
|
|
if (op.eq."0") then
|
|
y = 0.
|
|
else if (op.eq."1") then
|
|
y = 1.
|
|
else if (op.eq."P") then
|
|
y = 4.*atan(1.) ! pi
|
|
else
|
|
y = x(ichar(op)-96)
|
|
end if
|
|
else if (arity.eq.1) then ! This is a unary function
|
|
if (op.eq.">") then
|
|
y = stack(j) + 1
|
|
else if (op.eq."<") then
|
|
y = stack(j) - 1
|
|
else if (op.eq."~") then
|
|
y = -stack(j)
|
|
else if (op.eq."\") then
|
|
y = 1./stack(j)
|
|
else if (op.eq."L") then
|
|
y = log(stack(j))
|
|
else if (op.eq."E") then
|
|
y = exp(stack(j))
|
|
else if (op.eq."S") then
|
|
y = sin(stack(j))
|
|
else if (op.eq."C") then
|
|
y =cos(stack(j))
|
|
else if (op.eq."A") then
|
|
y = abs(stack(j))
|
|
else if (op.eq."N") then
|
|
y = asin(stack(j))
|
|
else if (op.eq."T") then
|
|
y = atan(stack(j))
|
|
else if (op.eq."O") then
|
|
y = 2.*stack(j)
|
|
else if (op.eq."J") then
|
|
y = 1+2.*stack(j)
|
|
else
|
|
y = sqrt(stack(j))
|
|
end if
|
|
else ! This is a binary function
|
|
if (op.eq."+") then
|
|
y = stack(j-1)+stack(j)
|
|
else if (op.eq."-") then
|
|
y = stack(j-1)-stack(j)
|
|
else if (op.eq."*") then
|
|
y = stack(j-1)*stack(j)
|
|
else
|
|
y = stack(j-1)/stack(j)
|
|
end if
|
|
end if
|
|
j = j + 1 - arity
|
|
stack(j) = y
|
|
! write(*,'(9f10.5)') (stack(k),k=1,j)
|
|
end do
|
|
if (j.ne.1) stop 'DEATH ERROR: STACK UNBALANCED'
|
|
f = stack(1)
|
|
!write(*,'(9f10.5)') 666.,x(1),x(2),x(3),f
|
|
return
|
|
end
|
|
|
|
subroutine multiloop(n,bases,i,done)
|
|
! Handles <n> nested loops with loop variables i(1),...i(n).
|
|
! Example: With n=3, bases=2, repeated calls starting with i=(000) will return
|
|
! 001, 010, 011, 100, 101, 110, 111, 000 (and done=.true. the last time).
|
|
! All it's doing is counting in mixed radix specified by the array <bases>.
|
|
implicit none
|
|
integer n, bases(n), i(n), k
|
|
logical done
|
|
done = .false.
|
|
k = 1
|
|
555 i(k) = i(k) + 1
|
|
if (i(k).lt.bases(k)) return
|
|
i(k) = 0
|
|
k = k + 1
|
|
if (k.le.n) goto 555
|
|
done = .true.
|
|
return
|
|
end
|
|
|
|
real*8 function limit(x)
|
|
implicit none
|
|
real*8 x, xmax
|
|
parameter(xmax=666.)
|
|
if (abs(x).lt.xmax) then
|
|
limit = x
|
|
else
|
|
limit = sign(xmax,x)
|
|
end if
|
|
return
|
|
end
|
|
|
|
subroutine LoadMatrixTranspose(nd,n,mmax,m,A,f)
|
|
! Reads the n x m matrix A from the file named f, stored as its transpose
|
|
implicit none
|
|
integer nd,mmax,n,m,j
|
|
real*8 A(nd,mmax)
|
|
character*60 f
|
|
open(2,file=f,status='old')
|
|
m = 0
|
|
555 m = m + 1
|
|
if (m.gt.mmax) stop 'DEATH ERROR: m>mmax in LoadVectorTranspose'
|
|
read(2,*,end=666) (A(j,m),j=1,n)
|
|
goto 555
|
|
666 close(2)
|
|
m = m - 1
|
|
print *,m,' rows read from file ',f
|
|
return
|
|
end
|
|
|
|
real*8 function mymedian(n,a)
|
|
implicit none
|
|
integer n,nmax, i
|
|
parameter(nmax=10000000)
|
|
real*8 a(n), b(nmax)
|
|
if (n.gt.nmax) stop 'DEATH ERROR: n>nmax in mymedian'
|
|
do i=1,n
|
|
b(i) = a(i)
|
|
end do
|
|
call sort(n,b)
|
|
i = nint((n+.5)/2)
|
|
if (i.eq.0) i=1
|
|
mymedian = b(i)
|
|
return
|
|
end
|
|
|
|
subroutine permutation(n,iarr) ! Return a random permutation of the first n integer:
|
|
integer iarr(n), idum, nmax, i
|
|
parameter(nmax=10000000)
|
|
real*8 arr(nmax), brr(nmax), ran1
|
|
if (n.gt.nmax) stop "PERMUTATION DEATH ERROR: nmax TOO SMALL"
|
|
idum = -666
|
|
do i=1,n
|
|
arr(i) = ran1(idum)
|
|
brr(i) = i
|
|
end do
|
|
call sort2(n,arr,brr)
|
|
do i=1,n
|
|
iarr(i) = nint(brr(i))
|
|
end do
|
|
return
|
|
end
|
|
|
|
SUBROUTINE sort(n,arr) ! Numerical Recipes Quicksort:
|
|
INTEGER n,M,NSTACK
|
|
REAL*8 arr(n)
|
|
PARAMETER (M=7,NSTACK=50)
|
|
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
|
|
REAL*8 a,temp
|
|
jstack=0
|
|
l=1
|
|
ir=n
|
|
1 if(ir-l.lt.M)then
|
|
do 12 j=l+1,ir
|
|
a=arr(j)
|
|
do 11 i=j-1,1,-1
|
|
if(arr(i).le.a)goto 2
|
|
arr(i+1)=arr(i)
|
|
11 continue
|
|
i=0
|
|
2 arr(i+1)=a
|
|
12 continue
|
|
if(jstack.eq.0)return
|
|
ir=istack(jstack)
|
|
l=istack(jstack-1)
|
|
jstack=jstack-2
|
|
else
|
|
k=(l+ir)/2
|
|
temp=arr(k)
|
|
arr(k)=arr(l+1)
|
|
arr(l+1)=temp
|
|
if(arr(l+1).gt.arr(ir))then
|
|
temp=arr(l+1)
|
|
arr(l+1)=arr(ir)
|
|
arr(ir)=temp
|
|
endif
|
|
if(arr(l).gt.arr(ir))then
|
|
temp=arr(l)
|
|
arr(l)=arr(ir)
|
|
arr(ir)=temp
|
|
endif
|
|
if(arr(l+1).gt.arr(l))then
|
|
temp=arr(l+1)
|
|
arr(l+1)=arr(l)
|
|
arr(l)=temp
|
|
endif
|
|
i=l+1
|
|
j=ir
|
|
a=arr(l)
|
|
3 continue
|
|
i=i+1
|
|
if(arr(i).lt.a)goto 3
|
|
4 continue
|
|
j=j-1
|
|
if(arr(j).gt.a)goto 4
|
|
if(j.lt.i)goto 5
|
|
temp=arr(i)
|
|
arr(i)=arr(j)
|
|
arr(j)=temp
|
|
goto 3
|
|
5 arr(l)=arr(j)
|
|
arr(j)=a
|
|
jstack=jstack+2
|
|
if(jstack.gt.NSTACK) stop 'NSTACK too small in sort'
|
|
if(ir-i+1.ge.j-l)then
|
|
istack(jstack)=ir
|
|
istack(jstack-1)=i
|
|
ir=j-1
|
|
else
|
|
istack(jstack)=j-1
|
|
istack(jstack-1)=l
|
|
l=i
|
|
endif
|
|
endif
|
|
goto 1
|
|
END
|
|
|
|
SUBROUTINE sort2(n,arr,brr) ! Numerical Recipes Quicksort:
|
|
INTEGER n,M,NSTACK
|
|
REAL*8 arr(n),brr(n)
|
|
PARAMETER (M=7,NSTACK=50)
|
|
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
|
|
REAL*8 a,b,temp
|
|
jstack=0
|
|
l=1
|
|
ir=n
|
|
1 if(ir-l.lt.M)then
|
|
do 12 j=l+1,ir
|
|
a=arr(j)
|
|
b=brr(j)
|
|
do 11 i=j-1,1,-1
|
|
if(arr(i).le.a)goto 2
|
|
arr(i+1)=arr(i)
|
|
brr(i+1)=brr(i)
|
|
11 continue
|
|
i=0
|
|
2 arr(i+1)=a
|
|
brr(i+1)=b
|
|
12 continue
|
|
if(jstack.eq.0)return
|
|
ir=istack(jstack)
|
|
l=istack(jstack-1)
|
|
jstack=jstack-2
|
|
else
|
|
k=(l+ir)/2
|
|
temp=arr(k)
|
|
arr(k)=arr(l+1)
|
|
arr(l+1)=temp
|
|
temp=brr(k)
|
|
brr(k)=brr(l+1)
|
|
brr(l+1)=temp
|
|
if(arr(l+1).gt.arr(ir))then
|
|
temp=arr(l+1)
|
|
arr(l+1)=arr(ir)
|
|
arr(ir)=temp
|
|
temp=brr(l+1)
|
|
brr(l+1)=brr(ir)
|
|
brr(ir)=temp
|
|
endif
|
|
if(arr(l).gt.arr(ir))then
|
|
temp=arr(l)
|
|
arr(l)=arr(ir)
|
|
arr(ir)=temp
|
|
temp=brr(l)
|
|
brr(l)=brr(ir)
|
|
brr(ir)=temp
|
|
endif
|
|
if(arr(l+1).gt.arr(l))then
|
|
temp=arr(l+1)
|
|
arr(l+1)=arr(l)
|
|
arr(l)=temp
|
|
temp=brr(l+1)
|
|
brr(l+1)=brr(l)
|
|
brr(l)=temp
|
|
endif
|
|
i=l+1
|
|
j=ir
|
|
a=arr(l)
|
|
b=brr(l)
|
|
3 continue
|
|
i=i+1
|
|
if(arr(i).lt.a)goto 3
|
|
4 continue
|
|
j=j-1
|
|
if(arr(j).gt.a)goto 4
|
|
if(j.lt.i)goto 5
|
|
temp=arr(i)
|
|
arr(i)=arr(j)
|
|
arr(j)=temp
|
|
temp=brr(i)
|
|
brr(i)=brr(j)
|
|
brr(j)=temp
|
|
goto 3
|
|
5 arr(l)=arr(j)
|
|
arr(j)=a
|
|
brr(l)=brr(j)
|
|
brr(j)=b
|
|
jstack=jstack+2
|
|
if(jstack.gt.NSTACK)stop 'NSTACK too small in sort2'
|
|
if(ir-i+1.ge.j-l)then
|
|
istack(jstack)=ir
|
|
istack(jstack-1)=i
|
|
ir=j-1
|
|
else
|
|
istack(jstack)=j-1
|
|
istack(jstack-1)=l
|
|
l=i
|
|
endif
|
|
endif
|
|
goto 1
|
|
END
|
|
|
|
! Numerical Recipes random number generator:
|
|
FUNCTION ran1(idum)
|
|
INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
|
|
REAL*8 ran1,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)
|
|
INTEGER j,k,iv(NTAB),iy
|
|
SAVE iv,iy
|
|
DATA iv /NTAB*0/, iy /0/
|
|
if (idum.le.0.or.iy.eq.0) then
|
|
idum=max(-idum,1)
|
|
do 11 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
|
|
11 continue
|
|
iy=iv(1)
|
|
endif
|
|
k=idum/IQ
|
|
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
|
|
|