c_______________________________________________________________________________
c_______________________________________________________________________________
c
c     Este programa calcula el factorial de un numero n almacenado en una
c     variable de tipo integer*4.
c
c     Se trabaja en base nbase = 10**logb (logb = 1, 2, 3, 4, 5 o 6).
c     Cada una de las cifras (en base nbase) del resultado (n!) se almacena
c     en la componente correspondiente de un vector integer*4 lc(0:MXC),
c     de forma que
c                                   
c                  n!        =    SUMA  [ lc(ic) * nbase**ic ].
c                                ic=0,nc
c
c     Con el fin de evitar un posible overflow al realizar las operaciones
c     aritmeticas en base nbase = 10**logb, como valor de logb se elige el
c     mayor posible en el rango {1,2,3,4,5,6} que satisface la condicion
c
c                  ( n * (10**logb - 1) .LE. INT_MAX ),
c
c     donde INT_MAX = 2**31 - 1 es el valor del entero mas grande que se
c     puede almacenar en una variable de tipo integer*4.
c
c
c     El numero de cifras necesarias en base 10 viene dado por la formula
c
c                  1 + nc    =   1 + int( log  (n!) ).
c                        10                  10
c
c     Para predecir este valor se utiliza la Aproximacion de Stirling
c
c                                      1/2      n
c                  n!        ~ (2 pi n)    (n/e) ,
c
c     de forma que
c
c                  1 + nc    ~  1 + (1/2) log  (2 pi n) + n log  (n/e).
c                        10                  10                10
c
c     Este valor se redondea al siguiente multiplo de 60 para facilitar
c     la impresion del resultado final (60 cifras por linea).
c
c     Finalmente, el numero de cifras necesarias en base 10**logb sera
c
c                  1 + nc    =  (1 + nc  ) / logb.
c                                      10
c
c     Teniendo en cuenta todo lo anterior, se obtiene la siguiente tabla de
c     valores posibles
c
c           min(n)        max(n)  =>  logb    nbase      min(nc) -       max(nc)
c                0 -       2.147  =>     6   1000000           9 -         1.039
c            2.148 -      21.475  =>     5    100000       1.247 -        16.751
c           21.476 -     214.769  =>     4     10000      20.939 -       262.979
c          214.770 -   2.149.633  =>     3      1000     350.639 -     4.226.239
c        2.149.634 -  21.691.754  =>     2       100   6.339.359 -    74.858.249
c       21.691.755 - 238.609.294  =>     1        10 149.716.559 - 1.895.367.419
c
c     En todo caso (1 + nc) sera menor o igual que (1 + MXC), que a su vez debera
c     cumplir la limitacion de direccionamiento
c
c                  ( 4 * (1 + MXC) .LE. (UINT_MAX + 1) ),
c
c     donde UINT_MAX + 1 = 2**32 es el numero maximo de bytes direccionables
c     mediante una palabra de un ordenador de 32 bits. Por tanto, teniendo en
c     cuenta estas limitaciones el valor maximo de nc sera 1.073.741.823,
c     que corresponde a un valor maximo de n = 139.273.802.
c 
c     El sistema operativo y/o el compilador pueden imponer restricciones
c     adicionales. El compilador gfortran utilizado en windows para compilar este
c     programa, por ejemplo, no permite dimensionar un vector entero de mas de
c     490.000.000 componentes (aproximadamente), por lo que el numero mas grande
c     cuyo factorial se puede calcular con esta version sera n = 67.351.429.
c     Finalmente, para que el programa pueda ser ejecutado en un ordenador,
c     el sistema debera disponer de suficiente memoria (central o virtual).
c
c     El tiempo de calculo se estima que crece como T(n^2 log(n)) por lo que
c     puede ser muy elevado para numeros grandes. Un pequeño
c     cambio en el valor de n puede producir un incremento apreciable del
c     tiempo de computacion necesario para realizar los calculos, debido
c     al importante incremento en el tamaño del vector lc() que se produce
c     al reducir la base.
c
c_______________________________________________________________________________
c_______________________________________________________________________________


      implicit real*8 (a-h,o-z)

c.....Establecemos el valor de algunas constantes
      parameter(INT_MAX=1+2*(2**30-1))           ! Maximo entero en un INTEGER*4
      parameter(PI=3.14159265358979324d+00,      ! Pi
     .          E =2.71828182845904524d+00)      ! e

c.....Vector donde se almacenara el resultado
      parameter(MXC=490000000-1)                 ! Limite de este sistema
      dimension lc(0:MXC)

c.....Valor de n
      write(0,'(/a,$)') ' n = '
      read(5,*) n

c.....Filtro de valores inadecuados
      if (n.lt.0) call FError(0,'n < 0', 1)
      write(6,'(/a,i10     )') ' Dato         :     n      = ',n 

c.....Seleccion de base 
      logb=0          ! logb = maximo entre 1 y 6 tal que (n*nbase.LE.INT_MAX),
      nbase=1         ! nbase=10**logb
      do while((n.le.(INT_MAX/(nbase*10-1))).and.logb.le.5)
        logb=logb+1
        nbase=nbase*10 
      enddo

      if (logb.eq.0) call FError(0, 'n demasiado grande para base 10',3)
      write(6,'( a,i10     )') ' Base         :     nbase  = ',nbase

c.....Estimacion del numero de cifras
      if (n.le.1) then
        nc10=0
      else
        dn=dble(n)
        nc10=int(0.5d+00*dlog10(2.d+00*PI*dn)+dn*dlog10(dn/E))
      endif
      nc10=(nc10/60)*60 +59     ! (1+nc10) se redondea al siguiente multiplo de 60
      nc=(1+nc10)/logb-1        ! (1+nc) = numero de cifras en base nbase

      write(6,'( a,i10,a,i8 )') ' No. de cifras:     1+nc   = ',1+nc,
     .                          ' en base ',nbase
      if (nbase.ne.10) 
     .write(6,'( a,i10,a,i8/)') '                    1+nc10 = ',1+nc10,
     .                          ' en base ',10

c.....Comprobacion de Memoria

      if (nc.gt.MXC) call FError(0,'memoria insuficiente',2)

c.....Inicializamos el factorial como 0!
      call Inicializar(nc,nc0,nca,lc)

c.....Calculamos el factorial
      do i=2,n
        if (((i/10000)*10000).eq.i) then
          write(0,100) i, 1+nca
  100     format(' ... calculando ',i10,'!     -> no. de cifras = ',i10)
        endif
        call Multiplicar(i,nbase,nc,nc0,nca,lc)
      enddo

c.....Escribimos el resultado
      call Escribir(n,logb,nc,nc0,nca,lc)

      write(0,'(/a,$)') ' Fin de Ejecucion. Oprima RETURN para salir.'
      read(5,'()')
      call exit(0)

      end


c____________________________________________________________________Inicializar
      Subroutine Inicializar(nc,nc0,nca,lc)
      dimension lc(0:nc)

      nc0=-1             ! ninguna cifra de la derecha es nula al principio
      nca=0              ! el resultado inicial (0!) cabe en la cifra de orden 0

      lc(0)=1            ! 0! = 1
      do ic=1,nc
        lc(ic)=0
      enddo

      return
      end


c____________________________________________________________________Multiplicar
      Subroutine Multiplicar(i,nbase,nc,nc0,nca,lc)
      dimension lc(0:nc)

      kp=0
      do ic=nc0+1,nca
        lc(ic)=kp+i*lc(ic)
        kp=lc(ic)/nbase
        lc(ic)=lc(ic)-kp*nbase;
      enddo
      do while(kp.gt.0.and.nca.lt.nc)
        nca=nca+1
        lc(nca)=kp
        kp=lc(nca)/nbase
        lc(nca)=lc(nca)-kp*nbase
      enddo
      if (kp.gt.0) call FError(0,'memoria insuficiente',2)

      do while(lc(nc0+1).eq.0)
        nc0=nc0+1
      enddo

      return
      end


c_______________________________________________________________________Escribir
      Subroutine Escribir(n,logb,nc,nc0,nca,lc)
      dimension lc(0:nc)
      character formato*8, buffer*60, linea*115

      write(6,100) n
  100 format(//i11,'! = '/)

      write(formato,'(a,i1,a,i1,a)') '(60i',logb,'.',logb,')'

      ncl=60/logb
      nclm=ncl-1
      do ic0=nc,nclm,-ncl
        ic1=ic0-nclm
        write(buffer,formato) (lc(ic),ic=ic0,ic1,-1)
        write(linea,200) (buffer(i:i),i=1,60), ic1*logb
  200   format(16x,6a1,9('.',6a1),10x,' * 10 ^ (',i10,')')
        write(6,'(a)') linea
      enddo
      write(6,'()')

      return
      end


c_________________________________________________________________________FError
      Subroutine FError(logerr,error,ierr)
      character*(*) error

      write(logerr,100) error
  100 format(/,' hiperfactorial:: Error: ',a)

      write(logerr,'(/a,$)')
     .      ' Fin de Ejecucion. Oprima RETURN para salir.'
      read(5,'()')
      call exit(ierr)

      return
      end
