C Compute A064097 and derived information.
C Compare with lengths of shortest addition chains
C Author: Hugo Pfoertner http:/www.pfoertner.org/
C
C Version history:
C Jan 31 2006 Initial version
C
C Definition of A064097
C a(1)=0, a(p)=1+a(p-1) if p = prime
C a(n*m) = a(n)+a(m) if m,n>1
C
      implicit integer (a-z)
      parameter (mmax =40000000 ,p2=40000000)
      dimension a(0:mmax), f(32), aclen(1:p2), a76091(mmax)
C
C Read table with lengths of shortest addition chains
      open ( unit=10, file='nclift.txt', form='formatted',
     &       status='old' )
      do 3 i = 1, p2
3     aclen(i) = 10000
C 25 numbers per line
      ilast = 0
      do 5 i = 1, p2, 25
      read (10,*,end=100) (aclen(k),k=i,i+24)
      ilast = i
5     continue
C
100   continue
      write (*,*)
     &  ' Shortest addition chain lengths computed by Neill Clift.'
      write (*,*) ' Chain length table read until:', ilast
C Initialize factorization
      call facini
C Terms of a064097
      a(1) = 0
      smax = 0
C
      write (*,*) ' Records in A064097:'
      do 10 n = 2, ilast
      call factor ( n, m, f )
      if ( m .eq. 1 ) then
C n is prime
        a(n) = 1 + a(n-1)
      else
C n is composite
        a(n) = a(f(1)) + a(n/f(1))
C Diagnostic output
C        write (*,*) n, f(1), n/f(1), a(f(1)), a(n/f(1))
      endif
C Records in A064097
      if ( a(n) .gt. smax ) then
        smax = a(n)
        write (*,1000) smax, n
1000    format ( i4, i9 )
      endif
     
10    continue
C Write first 100 sequence terms
      write (*,*) ' A064097:'
      write (*,1001) (a(k),k=1,100)
1001  format ( 31 i2 )
C
C Compare with length of shortest addition chain
      write (*,*)
     &  ' Increasing length difference to shortest addition chain:'
      diffm = -1
      do 20 i = 1, ilast
      diff = a(i) - aclen(i)
      if ( diff .gt. diffm ) then
        diffm = diff
	write (*,*) i, diff, a(i), aclen(i)
      endif
20    continue
C More terms for A076091
      onecnt = 0
      do 30 i = 1, ilast
      if ( a(i) - aclen(i) .eq. 1 ) then
        onecnt = onecnt + 1
        a76091(onecnt) = i
      endif
30    continue
40    continue
      write (*,*) ' A076091:'
      write (*,1002) (a76091(k),k=1,100)
1002  format ( 15 I5 )
C Check B. Cloitre's conjecture at n=2^k
      write (*,*) ' Check if limit approaches c=16.'
      write (*,*) ' A076091(n)*ln(n)/n:'
      do 50 j = 1, 24
      k = 2**j
      if ( k .gt. onecnt ) goto 55
      write (*,*) k, a76091(k), dble(a76091(k))*log(dble(k))/dble(k)
50    continue
55    continue
C largest available terms
      do 60 k = onecnt-9, onecnt
      write (*,*) k, a76091(k), dble(a76091(k))*log(dble(k))/dble(k)
60    continue
      end