c_______________________________________________________________________________
c
c     Este programa calcula el factorial de un numero n, relativamente grande,
c     almacenado en una 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^6 = 1.000.000,
c     cuyas cifras se almacenan en un vector integer*4 lc(NC), de forma que
c
c       n! = SUMA  [ lc(ic) * NBASE^(ic-1) ], con lc(ic) en {0,1,2,...,NBASE-1}. 
c           ic=1,NC
c
c     Al escribir el resultado se imprimen 60 cifras en base 10 (o lo que es
c     lo mismo: 10 cifras en base 10^6) por linea, por lo que se exige que NC
c     sea multiplo de 10.
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^NC        <-- para poder almacenar el resultado 
c
c     Debido a la primera limitacion el valor maximo de n sera
c
c       (INT_MAX+1) / NBASE = 2^31 / 1000000 = 2147  <-- division entera
c
c     cuyo factorial es aproximadamente
c
c       2147 ! ~  11,972078 * 10^6222 = 11,972078 * 1000000^1037
c
c     que consta de 1038 cifras en base 1.000.000.
c
c     En consecuencia, aunque se aumente el valor de NC por encima de 1038
c     no se podran calcular factoriales de numeros mayores que 2147.
c
c_______________________________________________________________________________

      parameter(NC=1040,                ! Numero de cifras (multiplo de 10)
     .          NBASE=1000000,          ! en base 1.000.000.
     .          INT_MAX=2*(2**30-1)+1)  ! Se calcula (2**31-1) sin sobrepasarlo

      dimension lc(NC)

      write(6,'(/a,$)') ' n = '
      read(5,*) n

      if (n.lt.0) then
        write(6,'(/a/)') ' Error -> n no puede ser negativo.'
        read(5,'()')
        call exit(1)
      else if (n.gt.INT_MAX/NBASE) then  ! Se evita (n*NBASE.gt.INT_MAX)
        write(6,'(/a,i12,a/)')
     .    ' El valor de n es demasiado grande para trabajar en base ',
     .    NBASE,' .'
        read(5,'()')
        call exit(3)
      endif

      call Inicializar(NC,lc)

      do i=2,n
        call Multiplicar(i,NBASE,NC,lc)
      enddo

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

      call Escribir(NC,lc)

      read(5,'()')
      call exit(0)

      end


c_______________________________________________________________________________
      Subroutine Inicializar(nc,lc)
      dimension lc(nc)

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

      return
      end

c_______________________________________________________________________________
      Subroutine Multiplicar(i,nbase,nc,lc)
      dimension lc(nc)

      kp=0
      do ic=1,nc
        lc(ic)=kp+i*lc(ic)
        kp=lc(ic)/nbase
        lc(ic)=lc(ic)-kp*nbase
      enddo
      if (kp.gt.0) then
        write(6,'(/a/)') ' Error -> El factorial es demasiado grande.'
        read(5,'()')
        call exit(2)
      endif

      return
      end

c_______________________________________________________________________________
      Subroutine Escribir(nc,lc)
      dimension lc(nc)

      do ic=nc,1,-1
        if (mod(nc-ic,10).eq.0) then
          write(6,100)
  100     format(/5x,$)
        endif
        write(6,200) lc(ic)
  200   format(i6.6,'.',$)
      enddo
      write(6,'()')

      return
      end
