c_______________________________________________________________________________
c_______________________________________________________________________________
c
c     Este programa calcula el factorial de un numero n almacenado en una
c     variable de tipo integer*4.
c
c     A continuacion espera a que se pulse la tecla RETURN (Entrar).
c
c     El factorial se codifica como un numero en base nbase = 10^logb (con
c     logb = 1, 2, 3, 4, 5 o 6) cuyas cifras se almacenan en un vector
c     integer*4 lc(0:MXC-1) DIMENSIONADO EN TIEMPO DE COMPILACION, de forma que
c
c       n! =   SUMA  [ lc(ic) * nbase^ic ], con lc(ic) en {0,1,2,...,nbase-1}. 
c           ic=0,nc-1
c
c     donde nc debe ser menor o igual que MXC.
c
c     Para calcular n! se procede iterativamente en la forma
c
c       1! = 1,    i! = i * (i-1)!   para i = 2, 3, ..., n.
c
c     Al calcular n * (n-1)! se multiplica n por las cifras de (n-1)!,
c     cada una de las cuales es un numero entero entre 0 y (nbase-1).
c
c     Por tanto, el valor maximo calculable esta limitado por las desigualdades
c
c       n*nbase < INT_MAX+1 = 2^31  <-- para evitar "overflows" 
c
c       n!      < 1000000^MXC       <-- para poder almacenar el resultado 
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.le.INT_MAX/nbase),   con nbase=10**logb
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 int.
c
c
c     El numero de cifras necesarias en base 10 viene dado por la formula
c
c       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       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        nc    =  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          10 -         1.040
c            2.148 -      21.474  =>     5    100000       1.248 -        16.752
c           21.475 -     214.748  =>     4     10000      20.940 -       262.950
c          214.749 -   2.147.483  =>     3      1000     350.600 -     4.221.700
c        2.147.484 -  21.474.836  =>     2       100   6.332.550 -    74.062.800
c       21.474.837 - 214.748.364  =>     1        10 148.125.600 - 1.696.004.340
c
c     En todo caso nc sera menor o igual que MXC, que a su vez debera
c     cumplir la limitacion de direccionamiento
c
c       (4*MXC.le.UINT_MAX+1),
c
c     donde UINT_MAX + 1 = 2**NBITS es el numero maximo de bytes direccionables
c     directamente en un procesador y con un compilador con longitudes de
c     palabra de NBITS bits. En un entorno de 32 bits esta restriccion limitaria
c     el valor maximo de los numeros cuyo factorial se podria calcular a
c     n = 139.273.802. En entornos de 64 bits esta restriccion no tendria
c     ningun efecto practico.
c
c     El sistema operativo y/o el compilador pueden imponer restricciones
c     adicionales. El compilador gfortran utilizado en windows para compilar
c     este programa no permite dimensionar en tiempo de compilacion un vector
c     entero de mas de 490.000.000 componentes (aproximadamente).
c     En consecuencia, el valor maximo de nc (multiplo de 60) seria
c     489.999.960, lo que corresponde a un valor maximo de n = 66.329.064.
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 (en 2022) 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(PI=3.14159265358979324d+00,      ! Pi
     .          E =2.71828182845904524d+00)      ! e
      parameter(INT_MAX=1+2*(2**30-1))           ! Maximo entero en un INTEGER*4


c.....Vector donde se almacenara el resultado
      parameter(MXC=490000000)        ! Limite aproximado para esta instalacion.

      dimension lc(0:MXC-1)

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
      nbase=1         ! (n*(nbase-1).LE.INT_MAX), con nbase=10**logb
      do while((n.le.(INT_MAX/(nbase*10))).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=1
      else
        dn=dble(n)
        nc10=1+int(0.5d+00*dlog10(2.d+00*PI*dn)+dn*dlog10(dn/E))
      endif
      nc10=((nc10-1)/60+1)*60   ! nc10 se redondea al siguiente multiplo de 60
      nc=nc10/logb              ! nc = numero de cifras en base nbase

      write(6,'( a,i10,a,i8 )') ' No. de cifras:     nc     = ',nc,
     .                          ' en base ',nbase
      if (nbase.ne.10) 
     .write(6,'( a,i10,a,i8/)') '                    nc10   = ',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,ic0,ica,lc)

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

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

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

      end


c____________________________________________________________________Inicializar
      Subroutine Inicializar(nc,ic0,ica,lc)
      dimension lc(0:nc-1)

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

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

      return
      end


c____________________________________________________________________Multiplicar
      Subroutine Multiplicar(i,nbase,nc,ic0,ica,lc)
      dimension lc(0:nc-1)

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

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

      return
      end


c_______________________________________________________________________Escribir
      Subroutine Escribir(n,logb,nc,lc)
      dimension lc(0:nc-1)
      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 ic_ini=nc-1,nclm,-ncl
        ic_fin=ic_ini-nclm
        write(buffer,formato) (lc(ic),ic=ic_ini,ic_fin,-1)
        write(linea,200) (buffer(i:i),i=1,60), ic_fin*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
