bigatom.e

-- -------------------------------------------------------------------------------------------
--    Copyright (C) 2014 Carlos Gómez Andreu (cargoan)
--
--    This program is free software: you can redistribute it and/or modify
--    it under the terms of the GNU General Public License as published by
--    the Free Software Foundation, either version 3 of the License, or
--    (at your option) any later version.
--
--    This program is distributed in the hope that it will be useful,
--    but WITHOUT ANY WARRANTY; without even the implied warranty of
--    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
--    GNU General Public License for more details.
--
--    You should have received a copy of the GNU General Public License
--    along with this program.  If not, see <http://www.gnu.org/licenses/>.
-- -------------------------------------------------------------------------------------------
--
-- bigatom.e
--    por Carlos J. Gómez Andreu (cargoan)
--
--    Mis números en coma flotante en base 10.
--
--    La función log() y relacionadas están adaptadas de bigfixedmath.e (Lucius L. Hilley III)
--    del Archivo de Euphoria.
--    Más precisa que la interna de euphoria y la que mejor ha ido con bigatoms.
--
--    La función exp() adaptada de las manpages de bc.
--
--    La función euler(), sobre el algoritmo no tengo ni idea, sacada de un ejemplo en Modula-2
--    de un viejo compilador de Modula-Oberon que, muuucho tiempo ha, utilizaba en OS/2 llamado
--    XDS. El algoritmo original fué escrito en Algol por Serge Batalov, reescrito en Modula-2
--    y modificado por Eugene Nalimov y Pavel Zemtsov.
--    Adaptado a Euphoria por mí (sin límite de decimales).
--
-- -------------------------------------------------------------------------------------------
--
-- Un bigatom es una sequence de tres elementos, el signo, el exponente y la mantisa.
-- El signo es 1 si el número es positivo, -1 si es negativo ó 0 si es cero.
--
-- El valor de un bigatom se evalúa en la forma:
--       valor = signo * mantisa * 10^exponente
-- El punto decimal va implícito tras el primer dígito de la mantisa.
--
-- Un mismo valor puede tener diferentes representaciones.
-- ej. el número -23.456 quedaría representado así:
--        { -1, 1, {2, 3, 4, 5, 6} } = -1 * 2.3456 * 10^1
--     o así:
--        { -1, 3, {0, 0, 2, 3, 4, 5, 6} } = -1 * 0.023456 * 10^3
--     también así:
--        { -1, 1, {2, 3, 4, 5, 6, 0, 0} } = -1 * 2.345600 * 10^1
--     etc.
--
-- Reduciéndose siempre a la forma más corta y almacenando únicamente los dígitos
-- significativos del número; eliminando los ceros a derecha e izquierda hasta el
-- primer dígito distinto de cero y ajustando el exponente en consecuencia.
--
-- La función normalize() es la encargada de realizar esta tarea.
--
-- -------------------------------------------------------------------------------------------

namespace bigatom

without warning &= { override }
without type_check

-- -------------------------------------------------------------------------------------------

-- máximo número de dígitos de un atom (para conversión binario => decimal)
ifdef BITS64 then
   constant ATOM_DIGS = 19   -- entre 18 y 20 más o menos (a gusto de cada uno)
elsedef
   constant ATOM_DIGS = 15   -- ??? (no sé si en 32 bits sean más)
end ifdef


-- estructura de un bigatom
public enum SIGN, EXPONENT, DIGITS
--
public type bigatom(object x)
   if length(x) = DIGITS
      and t_digits(x[DIGITS])
      and t_sign(x[SIGN])
      and integer(x[EXPONENT])
   then
      return 1
   end if

   return 0
end type
--


--
constant tdigits = {0,1,2,3,4,5,6,7,8,9}
--
type t_digits(object x)
   if sequence(x) then
      for i = 1 to length(x) do
         if not find(x[i], tdigits) then
            return 0
         end if
      end for
      return 1
   end if
   return 0
end type
--


-- signos válidos
type enum t_sign
   SG_NOVALUE = -2,
   SG_MINUS,
   SG_ZERO,
   SG_PLUS
end type
--


-- -------------------------------------------------------------------------------------------


-- escala
integer SCALE = 25    -- nº de decimales (posiciones o dígitos)
integer SC_MODE = 1   -- 0 son posiciones, no 0 son dígitos

-- NO_VALUE: no es un número real.
-- devuelto cuando un resultado es indefinido, indeterminado, un complejo, ...
-- ej. división por cero, logaritmo de cero o negativo, raíz de un negativo, ...
constant NO_VALUE = {SG_NOVALUE, 0, {}}

-- definición del cero
constant ZERO = {SG_ZERO, -1, {}}  -- ZERO = normalize({0,0,{0}})
-- de utilidad en
constant ONE  = {SG_PLUS, 0, {1}}  -- ONE  = ba_new(1)
--constant TWO  = {SG_PLUS, 0, {2}}  -- TWO  = ba_new(2)

-- caracteres permitidos en la entrada
constant SPACE     = ' '    -- carácter espacio
constant UNDERLINE = '_'    -- carácter subrayado

-- caracteres signos
constant SPLUS     = '+'    -- carácter signo positivo
constant SMINUS    = '-'    -- carácter signo negativo
constant SZERO     = SPACE  -- carácter signo de cero

integer DFCHAR = SPACE  -- carácter de relleno decimales (formato 'c')


-- -------------------------------------------------------------------------------------------


--
-- Establece el número de decimales y el modo (posiciones o dígitos)
--    decs: (atom) el número de decimales
--          (sequence) si es nula, devuelve los valores actuales
--                     si no, el primer elemento es el número de decimales
--                     y el segundo, si existe, es el modo
--
--    mode: si es cero, establece que son posiciones a partir del punto decimal
--          mayor que cero, dígitos (si 0 < n < 1 no cuenta ceros delanteros) esto permite
--                          manejar números muy pequeños sin tener que usar una escala muy grande
--          menor que cero, mantiene el modo anterior
--
--    Devuelve una sequence con los valores anteriores: { escala, modo }
--
export function scale(object decs = -1, integer mode = -1)

   if sequence(decs) then
      if length(decs) then
         if length(decs) > 1 then
            mode = not (not decs[2])
         end if
         decs = decs[1]
      else
         decs = -1
         mode = -1
      end if
   end if

   sequence prev = {SCALE, SC_MODE}
   if decs >= 0 then
      SCALE = floor(decs)
   end if
   if mode >= 0 then
      SC_MODE = not (not mode)
   end if

   return prev
end function
--


--
-- Devuelve la escala de un bigatom
-- (el número de decimales)
--
export function scale_of(bigatom N)
   integer decs = length(N[DIGITS]) - N[EXPONENT] - 1
   return decs * (decs > 0)
end function
--


--
-- Devuelve un bigatom con el valor de un número entero, un atom o, mejor aún,
-- una cadena con la representación de un número en coma fija o flotante que
-- no tiene pérdida de precisión por la conversión binario <=> decimal.
--
-- Permite el uso de separadores (espacio y subrayado) en la cadena para mejorar
-- la legibilidad... o puede que no.
--    Ej.: s = ba_new(1.2345)              -- s = {1,1,{1,2,3,4,5}}
--         s = ba_new("1 234 567.891_123") -- s = {1,6,{1,2,3,4,5,6,7,8,9,1,1,2,3}}
--         s = ba_new("-.2_3   45e-1 _3")  -- s = {-1,-14,{2,3,4,5}}
--
-- Devuelve NO_VALUE si la cadena es nula o no empieza con una representación de un
-- número válida. Ej.:
--    s = ba_new("e-231.2345 hola")  -- s = {-2,0,{}} = NO_VALUE
--    s = ba_new("-231.2345hola")    -- s = {-1,2,{2,3,1,2,3,4,5}}
--    s = ba_new("-231.2345e-23hola")              -- s = {1,-21,{2,3,1,2,3,4,5}}
--    s = ba_new("__ -231.2345 e-12e++..--Eholae") -- s = {1,-10,{2,3,1,2,3,4,5}}
--
export function ba_new(object number)

   if bigatom(number) then
      return number
   elsif atom(number) then
      if integer(number) then
         return int_to_bigatom(number)
      end if

      atom absn = number
      if absn < 0 then
         absn = -absn
      end if

      integer exp, len, ndigits
      -- ndigits = -floor(-logb(absn), 10)
      ndigits = -floor(-eu:log(absn)/eu:log(10))  -- suficiente
      if ndigits >= 0 then
         if ndigits > ATOM_DIGS then
            ndigits = ATOM_DIGS
         end if
         number = sprintf(sprintf("%%+.%df", ATOM_DIGS - ndigits), number)
         len = find('.', number)
         if not len then
            len = length(number) - 1
         end if
         if len > ATOM_DIGS then
            number  = int_value(number[1..ATOM_DIGS + 1])
            number  = sprintf("%+.f", number)
            number &= repeat('0', len - length(number) + 1)
         end if
      else
         sequence snum = sprintf("%.e", number)
         exp = find('e', snum)
         if exp then
            exp = int_value(snum[exp+1..$])
         end if
         number = sprintf(sprintf("%%+.%df", ATOM_DIGS - exp), number)
      end if
   end if

   sequence big = NO_VALUE

   -- elimina separadores válidos (espacio y subrayado)
   -- se detiene con el primer carácter no válido y descarta el resto
   object c
   integer pos = 0, sflag = 0
   for i = 1 to length(number) do
      c = number[i]
      if integer(c) then
         if find(c, "+-") then
            if sflag then
               exit
            end if
            pos += 1
            number[pos] = c
            sflag = 1
         elsif find(c, ".0123456789eE") then
            pos += 1
            number[pos] = c
         elsif c = SPACE or c = UNDERLINE then
            continue
         else
            exit
         end if
      else
         exit
      end if
   end for
   if pos then
      while pos and find(number[pos], "+-.eE") do
         pos -= 1
      end while
   end if
   number = number[1..pos]   -- remove(number, pos+1, length(number))

   -- coma flotante
   pos = find('e', number) + find('E', number)
   if pos then
      big = ba_new(number[1..pos-1])
      big[EXPONENT] += int_value(number[pos+1..$])

      return normalize(big)
   end if

   -- coma fija
   if length(number) then
      if number[1] = SMINUS then
         big[SIGN] = SG_MINUS
         number = remove(number, 1)
      else
         big[SIGN] = SG_PLUS
         if number[1] = SPLUS then
            number = remove(number, 1)
         end if
      end if

      pos = find('.', number)
      if pos then
         number = remove(number, pos)
         big[EXPONENT] = pos - 2
      else
         big[EXPONENT] = length(number) - 1
      end if
      big[DIGITS] = number - '0'
   end if

   return normalize(big)
end function
--


--
-- Compara dos bigatoms
-- (NO_VALUE es menor que cualquier número)
--
export function ba_compare(bigatom A, bigatom B)
   integer cmp = compare(A, B)

   if A[SIGN] = SG_MINUS and B[SIGN] = SG_MINUS then
      cmp = -cmp
   end if

   return cmp
end function
--


-- -------------------------------------------------------------------------------------------


--
-- ###############################
-- ###   FUNCIONES DE SALIDA   ###
-- ###############################
--

--
-- Devuelve la representación completa de un bigatom en una cadena
--
export function ba_sprint(bigatom N)
   sequence str = to_string(N)
   if str[1] = SZERO then
      str = remove(str, 1)
   end if

   return str
end function
--


--
-- Escribe la representación completa de un bigatom en un archivo
-- (1 = STDOUT, 2 = STDERR)
export procedure ba_print(integer file, bigatom N)
   puts(file, ba_sprint(N))
end procedure
--


--
-- Devuelve una cadena con la representación con formato de un bigatom
-- Nota: Necesita una revisión a fondo, está hecha sobre la marcha pero es funcional.
--       Sería bueno integrarla con otros tipos, substituir los bigatoms y pasarle todo
--       a printf o sprintf. En realidad esta función solo lee la cadena de formato y
--       llama a to_string() que es quien realiza la conversión del número y ésta ajusta
--       al tamaño, pone el relleno.
--
-- La cadena de formato está formada por tres partes: la cabecera (el texto que va delante
-- del número), el propio formato y la cola (el texto que se pondrá detrás). Todas opcionales.
--
-- Si la cadena es nula, se representará el número con todos sus dígitos. Si no, el formato se
-- se delimita con la pareja de caracteres '%' y 'B'. La parte anterior y posterior al formato
-- se añaden delante y detrás respectivamente. Si no contiene un formato %B, la cadena se coloca
-- delante del número.
--
-- Es similar a printf(), las diferencias son:
--   - El tamaño se refiere solo a la parte entera, no al total. Si delante del tamaño hay un
--     carácter que no sea un dígito entre 1 y 9 o el signo, se toma como carácter de relleno
--     para completar la parte entera.
--
--   - Si detrás del punto hay un cero, se pondrá el signo si el número no es cero, pero sí en
--     la representación.
--
--   - Si el formato finaliza en 'a' se rellenan con ceros las posiciones decimales no ocupadas.
--     Si finaliza en 'c' se utilizan espacios pero si va precedido de un carácter que no sea
--     un dígito, se utiliza éste para rellenar las posiciones decimales.
--
--     'a' y 'c' se usan principalmente para producir salidas encolumnadas.
--     'c' aumenta notablemente la legibilidad al poner solo los decimales necesarios. Si no se
--      especifica un tamaño 'c' es ignorado. No así 'a'.
--
--
-- fmt = '%' ['+'] [fchar] [size] [['.'] ['0'] [decs] ['a' | [char] 'c']] 'B'
--
--          [+] poner signo siempre (*)
--
--      [fchar] carácter de relleno (cualquiera excepto dígitos del 1 al 9, sí el 0)
--
--       [size] tamaño mínimo de la parte entera del número (si no cabe es ignorado)
--
--          [.] un punto, tan solo separa las partes de la cadena de formato
--
--          [0] poner signo cuando el número es menor que los límites de representación pero
--              no cero en la escala actual, resultando ceros con signo
--
--       [decs] el número "máximo" de decimales a mostrar (si es negativo, todos los decimales)
--
--  [a|[char]c] 'a': rellenar con ceros las posiciones decimales no ocupadas.
--              'c': usa el carácter que le precede (si no es un dígito), o un espacio como
--                   carácter de relleno para la parte decimal. Si no se ha especificado un
--                   tamaño 'c' es ignorado.
--
-- Cuando el carácter de relleno es '0' se reserva la primera posición para el signo.
--
-- (*) El signo de cero (un espacio) solo se muestra cuando el carácter de relleno es '0' y
--     el signo se coloca delante del relleno y no detrás como en los demás casos.
--
export function ba_sprintf(sequence fmt, bigatom N)

   if not length(fmt) then
      return ba_sprint(N)
   end if

   sequence head, tail
   integer fpos = find('%', fmt)
   if fpos then
      head = fmt[1..fpos-1]
      fmt  = remove(fmt, 1, fpos)
      fpos = find('B', fmt)
      if fpos then
         tail = fmt[fpos+1..$]
         fmt  = fmt[1..fpos-1]
      else
         tail = fmt
         fmt  = ""
      end if
   else
      head = fmt
      fmt  = ""
      tail = ""
   end if

   sequence ifmt, dfmt
   integer decs = -1, dp = find('.', fmt)
   if dp then
      ifmt = fmt[1..dp-1]
      dfmt = fmt[dp+1..$]
      decs = 0
   else
      ifmt = fmt
      dfmt = ""
   end if

   integer c, sg = 0, size = 0, fill = 0, fchar = SPACE
   while length(ifmt) do
      c = ifmt[1]
      if not sg and c = SPLUS then
         sg = 1
      elsif c < '1' or c > '9' then
         if fchar = SPACE then
            fchar = c
         end if
      else
         size = int_value(ifmt)
         exit
      end if
      ifmt = remove(ifmt, 1)
   end while

   integer zsign = 0, all = 0
   if length(dfmt) then
      c = dfmt[$]
      if c = 'a' then
         all = 1
      elsif c = 'c' then
         c = dfmt[$-1]
         if size then
            all = 2
            if c < '0' or c > '9' then
               DFCHAR = c
               dfmt = dfmt[1..$-1]
            else
               DFCHAR = SPACE
            end if
         end if
      end if
      if all then
         dfmt = dfmt[1..$-1]
      end if

      while length(dfmt) do
         c = dfmt[1]
         if not zsign and c = '0' then
            zsign = 1
         else
            decs = int_value(dfmt)
            exit
         end if
         dfmt = remove(dfmt, 1)
      end while
   end if

   if decs < 0 then
      decs = scale_of(N)
   end if

   sequence str = to_string(N, sg, decs, zsign, all)

   -- elimina el signo de cero si el relleno no es cero
   if str[1] = SZERO and fchar != '0' then
      if size then
         str[1] = fchar
      else   --  if not sg and fchar != SPACE then
         str = remove(str, 1)
      end if
   elsif str[1] = '<' and size and all then
      if all > 1 then   -- 'c'
         str = str & repeat(DFCHAR, decs + 1)
      else              -- 'a' (usa relleno parte entera)
         str = str & repeat(fchar, decs + 1)
      end if
   end if

   -- si no hay decimales, substituye el punto por el relleno
   integer len = find('.', str)
   if len and str[len + 1] = DFCHAR then
      str[len] = DFCHAR
   end if

   if len then
      decs = length(str) - len
      len -= 1
   else
      len = length(str)
      decs = 0
   end if

   -- relleno parte entera
   len = size - len - decs - (decs > 0)
   if len > 0 then
      if fchar = '0' and str[1] != '<' then
         str = str[1] & repeat(fchar, len) & str[2..$]
      elsif str[1] = '<' then
         str = repeat(DFCHAR, len) & str
      else
         str = repeat(fchar, len) & str
      end if
   end if

   return head & str & tail
end function
--


--
-- Escribe la representación con formato de un bigatom en un archivo
-- (1 = STDOUT, 2 = STDERR)
--
export procedure ba_printf(integer file, sequence fmt, bigatom N)
   puts(file, ba_sprintf(fmt, N))
end procedure
--

-- -------------------------------------------------------------------------------------------


-- ###################################
-- ###   OPERACIONES ARITMÉTICAS   ###
-- ###################################


--
-- Devuelve un bigatom con el resultado de la suma de dos números.
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_add(object A, object B)

   if not bigatom(A) then
      A = ba_new(A)
   end if
   if not bigatom(B) then
      B = ba_new(B)
   end if

   integer signA = A[SIGN], signB = B[SIGN]
   if signA = SG_ZERO then
      return B
   elsif signB = SG_ZERO then
      return A
   elsif signA = SG_NOVALUE or signB = SG_NOVALUE then
      return NO_VALUE
   end if

   integer sign = SG_PLUS
   if signA = signB then
      sign = signA
   elsif signA = SG_MINUS then
      A[SIGN] = SG_PLUS
      return ba_sub(B, A)
   else
      B[SIGN] = SG_PLUS
      return ba_sub(A, B)
   end if

   {A, B} = align(A, B)

   integer exponent = A[EXPONENT]
   sequence res = digits_add(A[DIGITS], B[DIGITS])
   if res[1] then
      exponent += 1
   end if
   res = {sign, exponent, res[2]}

   return  normalize(res)
end function
--


--
-- Devuelve un bigatom con el resultado de la resta de dos números.
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_sub(object A, object B)

   if not bigatom(A) then
      A = ba_new(A)
   end if
   if not bigatom(B) then
      B = ba_new(B)
   end if

   integer signA = A[SIGN], signB = B[SIGN]
   if signA = SG_ZERO then
      B[SIGN] = -signB
      return B
   elsif signB = SG_ZERO then
      return A
   elsif signA = SG_NOVALUE or B[SIGN] = SG_NOVALUE then
      return NO_VALUE
   end if

   if signA = SG_MINUS or signB = SG_MINUS then
      B[SIGN] = -signB
      return ba_add(A, B)
   end if

   {A, B} = align(A, B)

   sequence res = digits_sub(A[DIGITS], B[DIGITS])
   integer sign = SG_PLUS
   if res[1] then
      sign = SG_MINUS
   end if
   res = normalize({sign, A[EXPONENT], res[2]})

   return  res
end function
--


--
-- Devuelve un bigatom con el resultado de la multiplicación de dos números.
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_multiply(object A, object B, integer round = 0)

   if not bigatom(A) then
      A = ba_new(A)
   end if
   if not bigatom(B) then
      B = ba_new(B)
   end if

   integer signA = A[SIGN], signB = B[SIGN]
   integer signR = signA * signB
   if signA = SG_ZERO or signB = SG_ZERO then
      return ZERO
   elsif signA = SG_NOVALUE or signB = SG_NOVALUE then
      return NO_VALUE
   end if

   sequence res, digsA = A[DIGITS], digsB = B[DIGITS]
   integer  expA = A[EXPONENT], expB = B[EXPONENT], expR = expA + expB
   if equal({1}, digsA) then       -- múltiplo de 10
      res = {signR, expR, digsB}
   elsif equal({1}, digsB) then    -- múltiplo de 10
      res = {signR, expR, digsA}
   else
      res = digits_multiply(digsA, digsB)
      res = {signR, expR + (res[1] > 0), res[2]}
   end if

   -- límite decimales
   integer ndecs
   if SC_MODE and res[EXPONENT] < 0 then
      ndecs = SCALE + 1
   else
      ndecs = res[EXPONENT] + SCALE + 2
   end if
   if ndecs > 0 then
      if round then
         res = round_digits(res, ndecs)
      else
         res[DIGITS] = remove(res[DIGITS], ndecs, length(res[DIGITS]))
         res = normalize(res)
      end if
   else
      res = ZERO
   end if

   return res
end function
--


--
-- Devuelve el resultado de la división entera de  dos números.
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
-- los números son truncados antes de hacer la división
--
export function ba_idivide(object A, object B)

   if not bigatom(A) then
      A = ba_trunc(ba_new(A))
   end if
   if not bigatom(B) then
      B = ba_trunc(ba_new(B))
   end if

   integer signA = A[SIGN], signB = B[SIGN]
   if signB = SG_ZERO or
      signB = SG_NOVALUE or
      signA = SG_NOVALUE
   then
      return NO_VALUE
   elsif signA = SG_ZERO then
      return ZERO
   end if

   sequence res, digsA
   sequence quotient = {}, partial = {SG_PLUS,-1,{}}   -- +0
   integer signR, expR = -1

   signR = A[SIGN] * B[SIGN]
   if equal({1}, B[DIGITS]) then    -- múltiplo de 10
      res = {signR, A[EXPONENT] - B[EXPONENT], A[DIGITS]}
   else
      A[SIGN] = SG_PLUS
      B[SIGN] = SG_PLUS

      A = expand(A)
      digsA = A[DIGITS]
      for i = 1 to length(digsA) do
         partial[DIGITS] &= digsA[i]
         partial[EXPONENT] += 1
         quotient &= 0
         expR += 1
         partial = normalize(partial)
         while compare(partial, B) != SG_MINUS do
            quotient[$] += 1
            partial = ba_sub(partial, B)
         end while
         partial[SIGN] = SG_PLUS    -- por si se ha puesto a cero
         partial = expand(partial)
      end for
      res = normalize({ signR, expR, quotient })
   end if

   return ba_floor(res)
end function
--


--
-- Devuelve un bigatom con el resultado de la división de dos números.
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_divide(object A, object B, integer round = 0)

   if not bigatom(A) then
      A = ba_new(A)
   end if
   if not bigatom(B) then
      B = ba_new(B)
   end if

   integer signA = A[SIGN], signB = B[SIGN]
   if signB = SG_ZERO or
      signB = SG_NOVALUE or
      signA = SG_NOVALUE
   then
      return NO_VALUE
   elsif signA = SG_ZERO then
      return ZERO
   end if

   sequence res
   integer decsB, ndecs
   integer expA = A[EXPONENT], expB = B[EXPONENT]
   if equal({1}, B[DIGITS]) then
      A[EXPONENT] -= expB
      A[SIGN] *= signB
      res = A
   else
      ndecs = scale_of(A)
      decsB = scale_of(B)
      if ndecs < decsB then
         ndecs = decsB
      end if

      integer mult = SCALE + 1  -- uno extra para redondeo
      if SC_MODE then
         if expA > expB then
         mult += expA - expB
         else
            mult += expB - expA
         end if
      end if
      A[EXPONENT] += ndecs + mult
      B[EXPONENT] += ndecs

      res = ba_idivide(A, B)
      res[EXPONENT] -= mult
   end if

   if SC_MODE and res[EXPONENT] < 0 then
      ndecs = SCALE + 1
   else
      ndecs = res[EXPONENT] + SCALE + 2
   end if
   if ndecs > 0 then
      if round then
         res = round_digits(res, ndecs)
      else
         res[DIGITS] = remove(res[DIGITS], ndecs, length(res[DIGITS]))
         res = normalize(res)
      end if
   else
      res = ZERO
   end if

   return res
end function
--


-- -------------------------------------------------------------------------------------------


-- #####################################
-- ###   RESTOS, REDONDEOS, VARIOS   ###
-- #####################################

--
-- Devuelve un bigatom con el resto de la división de dos números.
--
-- Admite atoms, bigatoms o cadenas que representen números
--
export function ba_remainder(object A, object B)

   if not bigatom(A) then
      A = ba_new(A)
   end if
   if not bigatom(B) then
      B = ba_new(B)
   end if

   if A[SIGN] = SG_NOVALUE or
      B[SIGN] = SG_NOVALUE or
      B[SIGN] = SG_ZERO
   then
      return NO_VALUE
   end if

   sequence res = ZERO
   integer comp = compare(ba_abs(A), ba_abs(B))
   if comp > 0 then
      res = ba_trunc(ba_divide(A, B))
      res = ba_multiply(res, B)
      res = ba_sub(A, res)
   elsif comp < 0 then
      res = A
   end if

   return res
end function
--


--
-- mi función de redondeo adaptada a bigatom
-- redondea igual lo números negativos que los positivos
--    -0.5 --> -1 en lugar de -0.5 --> 0
--
export function ba_round(object N, atom precision = 1, integer mode = 0)

   if not bigatom(N) then
      N = ba_new(N)
   end if

   if N[SIGN] = SG_NOVALUE then
      return NO_VALUE
   end if

   if not precision then
      precision = 1
   elsif precision < 0 then
      precision = -1 / precision
   end if

   sequence res
   if mode > 0 then
   -- n = ceil(n * precision) / precision
      res = ba_ceil(ba_multiply(N, precision))
   elsif mode < 0 then
   -- n = floor(n * precision) / precision
      res = ba_floor(ba_multiply(N, precision))
   else -- mode = 0
   -- n = sign * floor(sign * a * precision + 0.5) / precision
      res = ba_multiply(N, N[SIGN])
      res = ba_multiply(res, precision)
      res = ba_add(res, {SG_PLUS, -1, {5}})   -- 0.5
      res = ba_multiply(ba_floor(res), N[SIGN])
   end if
   res = ba_divide(res, precision, 1)

   return res
end function
--


--
-- devuelve el valor absoluto de un bigatom
--
export function ba_abs(bigatom N)
   if N[SIGN] = SG_MINUS then
      N[SIGN] = SG_PLUS
   end if

   return N
end function
--


--
-- devuelve la parte entera de un bigatom
--
export function ba_trunc(bigatom N)
   if N[SIGN] = SG_NOVALUE then
      return NO_VALUE
   elsif N[EXPONENT] < 0 then
      return ZERO
   end if

   N[DIGITS] = remove(N[DIGITS], N[EXPONENT] + 2, length(N[DIGITS]))

   return normalize(N)
end function
--


--
-- devuelve la parte decimal de un bigatom
--
export function ba_frac(bigatom N)
   if N[SIGN] = SG_NOVALUE then
      return NO_VALUE
   end if

   if N[EXPONENT] < 0 then
      return N
   end if
   N[DIGITS] = remove(N[DIGITS], 1, N[EXPONENT] + 1)
   N[EXPONENT] = -1

   return normalize(N)
end function
--


--
-- redondea al entero igual o el inmediato menor
--
export function ba_floor(bigatom N)

   if N[SIGN] = SG_NOVALUE then
      return NO_VALUE
   end if

   sequence I = ba_trunc(N)
   if N[SIGN] = SG_MINUS and compare(N, I) then
      I = ba_sub(I, ONE)
   end if

   return I
end function
--


--
-- redondea al entero igual o el inmediato mayor
-- (ceil(x) = -floor(-x) )
export function ba_ceil(bigatom N)

   if N[SIGN] = SG_NOVALUE then
      return NO_VALUE
   end if

   N[SIGN] = -N[SIGN]
   N = ba_floor(N)
   N[SIGN] = -N[SIGN]

   return N
end function
--


-- -------------------------------------------------------------------------------------------


-- ##################################################
-- ###   FUNCIONES EXPONENCIALES Y LOGARÍTMICAS   ###
-- ##################################################


--
-- adaptada de bigfixedmath.e (Lucius L. Hilley III) del Archivo de Euphoria
-- (más precisa que la interna de euphoria y la que mejor ha ido con bigatom)
--
constant LIMIT   = 1.1
--
override function log(atom x)

   if x <= 0 then       -- < 0 no es un real, es complejo (= log(x)+pi*i )
      return eu:log(x)  -- = 0 es indefinido (nan)
   end if

   integer inc
   atom res, limit, curr, prev
   if x > LIMIT then
      res = 0
      limit = bigatom:log(LIMIT)
      while x > LIMIT do
         x  /= LIMIT
         res += limit
      end while
      res += bigatom:log(x)

      return res
   end if

   prev = x
   x    = x - 1
   res  = x
   curr = x
   inc  = 1
   while res != prev do
      prev  = res
      curr *= x
      inc  += 1
      res  -= curr / inc
      curr *= x
      inc  += 1
      res  += curr/ inc
   end while

   return res
end function
--

export function logb(atom x, integer base = 10)
   if base > 0 then
      base = -base
   end if
   base = -floor(base)

   return log(x) / log(base)
end function
--


--
-- cálculo del número de Euler (e) y sus potencias
-- e = 2.718281828459045235360287471352662497757247093699959574
--       9669676277240766303535475945713821785251664274274663919...
-- adaptada de la librería de bc
--
-- override por si se ha incluido math.e
--
override function exp(atom x)

   integer neg = 0, mult = 0
   if x < 0 then
     neg = 1
     x   = -x
   end if

   while x > 1 do
     mult += 1
     x    /= 2
   end while

   integer inc = 2
   atom xpow = x, curr = x
   atom fact = 1, res = 1 + x
   while curr do
      xpow *= x
      fact *= inc
      curr  = xpow / fact
      res  += curr
      inc  += 1
   end while

   while mult do
      res  *= res
      mult -= 1
   end while

   if neg then
      res = 1 / res
   end if

   return res
end function
--


--
-- raíz cuadrada
--
override function sqrt(atom x)
   if x < 0 then
      return eu:sqrt(x)   -- no es un real, es imaginario (= sqrt(-x)i )
   elsif x = 0 then
      return 0
   elsif x = 1 then
      return 1
   end if

   atom res, res1
   if x < 1 then
      res = 1
   else
      res1 = -floor(-eu:log(x) * 0.5)
      res  = power(10, res1)
   end if

   while res != res1 do
      res1 = res
      res  = x / res
      res += res1
      res *= 0.5
   end while

   return res
end function
--


-- -------------------------------------------------------------------------------------------
--
-- -------------------------------------------------------------------------------------------



--
constant BLIMIT = {1,0,{1,1}}   -- 1.1
--
-- Devuelve un bigatom con el logaritmo natural
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
-- adaptada de bigfixedmath.e (Lucius L. Hilley III) del archivo de euphoria
-- (la que mejor ha rendido con bigatoms)
--
export function ba_log(object x, integer round = 0)

   if not bigatom(x) then
      x = ba_new(x)
   end if

   if x[SIGN] < SG_PLUS then
      return NO_VALUE
   end if

   sequence sc = scale(SCALE + 4, 0)

   sequence limit, ln_limit
   sequence res, curr, prev, inc
   if compare(ONE, x) = SG_PLUS then
      res = ba_log(ba_divide(ONE, x))
      res[SIGN] = -res[SIGN]
   elsif compare(x, BLIMIT) = SG_PLUS then
      res      = ZERO
      ln_limit = ba_log(BLIMIT)
      while compare(x, BLIMIT) = SG_PLUS do
         x   = ba_divide(x, BLIMIT)
         res = ba_add(res, ln_limit)
      end while
      res = ba_add(res, ba_log(x))
   else
      -- ln(x) = x - x^2/2 + x^3/3 - x^4/4 + x^5/5 - ...
      prev = x
      x    = ba_sub(x, ONE)
      curr = x
      res  = x
      inc  = {1, 0, {2}}  -- 2
      while compare(prev, res) do
         prev = res
         curr = ba_multiply(curr, x)
         res  = ba_sub(res, ba_divide(curr, inc))
         inc  = ba_add(inc, ONE)
         curr = ba_multiply(curr, x)
         res  = ba_add(res, ba_divide(curr, inc))
         inc  = ba_add(inc, ONE)
      end while
   end if

   scale(sc)

   if round then
      res = round_digits(res, res[EXPONENT] + 3 + SCALE)
   else
      res[DIGITS] = remove(res[DIGITS],
                             res[EXPONENT] + 2 + SCALE,
                             length(res[DIGITS]))
      res = normalize(res)
   end if

   return res
end function
--


--
-- Devuelve un bigatom con la potencia de e (e^x)
-- e = 2.71828182845904523536028747135266249775724709369995
--       95749669676277240766303535475945713821785251664274...
-- ba_exp(1) = e
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
-- adaptada de la librería de bc
--
export function ba_exp(object x, integer round = 0)

   if not bigatom(x) then
      x = ba_new(x)
   end if

   if equal(x, ONE) then
      return euler(SCALE, 1)  -- muchísimo más rápida
   end if

   integer neg = 0, mult = 0
   sequence sc, limit, inc = {1, 0, {2}}
   sequence xpow, curr, res, fact

   if x[SIGN] = SG_MINUS then
      neg = 1
      x[SIGN] = -x[SIGN]
   end if

   sc = scale(SCALE + 2, 0)

   while compare(x, ONE) = SG_PLUS do
      mult += 1
      x = ba_divide(x, {1, 0, {2}})
   end while

   curr = x
   xpow = x
   res  = ba_add(x, ONE)
   fact = ONE
   while curr[SIGN] do
      xpow = ba_multiply(xpow, x)
      fact = ba_multiply(fact, inc)
      curr = ba_divide(xpow, fact)
      res  = ba_add(res, curr)
      inc  = ba_add(inc, ONE)
   end while

   while mult do
      res   = ba_multiply(res, res)
      mult -= 1
   end while
   if neg then
      res = ba_divide(ONE, res)
   end if
   scale(sc)

   if round then
      res = round_digits(res, res[EXPONENT] + 2 + SCALE)
   else
      res[DIGITS] = remove(res[DIGITS],
                             res[EXPONENT] + 2 + SCALE,
                             length(res[DIGITS]))
      res = normalize(res)
   end if

   return res
end function
--


--
-- Devuelve un bigatom con la potencia de x (x^exp)
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_power(object x, object exponent, integer round = 0)
   sequence res

   if not bigatom(x) then
      if atom(x) and atom(exponent) then
         if floor(exponent) = exponent then
            res = intf_power(ba_new(x), exponent)
         elsif exponent > 0 then
            res = ba_exp(ba_multiply(ba_log(x), exponent), round)
         elsif exponent then
            res = ba_divide(ONE, ba_exp(ba_multiply(ba_log(x), -exponent), round))
         else
            res = ONE
         end if

         return res
      end if
      x = ba_new(x)
   end if

   if atom(exponent) and floor(exponent) = exponent then
      res = intf_power(x, exponent)

      return res
   end if

   if not bigatom(exponent) then
      exponent = ba_new(exponent)
   end if

   if equal(ba_floor(exponent), exponent) then
      res = intf_power(x, bigatom_to_atom(exponent))
   elsif exponent[SIGN] = SG_PLUS then
      res = ba_exp(ba_multiply(ba_log(x), exponent), round)
   elsif exponent[SIGN] = SG_MINUS then
      exponent[SIGN] = -exponent[SIGN]
      res = ba_divide(ONE, ba_exp(ba_multiply(ba_log(x), exponent)), round)
   else
      res = ONE
   end if

   return res
end function
--


--
-- Devuelve un bigatom con la raíz cuadrada de x
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_sqrt(object x, integer round = 0)

   if not bigatom(x) then
      x = ba_new(x)
   end if

   if x[SIGN] < 0 then
   -- no es un real, es imaginario (= sqrt(-x)*i )
      return NO_VALUE
   elsif x[SIGN] = SG_ZERO then
      return ZERO
   end if

   integer cmp = compare(x, ONE)
   if not compare(x, ONE) then
      return ONE
   end if

   sequence res, res1
   if cmp < 1 then
   -- si 0 < x < 1, inicialmente 1
      res = ONE
   else
   -- si x > 1, inicio en 10^(exp/2) = 10^(exp*0.5)
      res1 = ba_floor(ba_multiply(x[EXPONENT] + 1, {1,-1,{5}}))
      res  = ba_power({1,1,{1}}, res1)
   end if

   while compare(res, res1) do
      res1 = res
      res = ba_divide(x, res)
      res = ba_add(res, res1)
      res = ba_multiply(res, {1,-1,{5}})   -- 0.5
   end while

   if round then
      res = round_digits(res, res[EXPONENT] + 2 + SCALE)
   else
      res[DIGITS] = remove(res[DIGITS],
                           res[EXPONENT] + 2 + SCALE,
                           length(res[DIGITS]))
      res = normalize(res)
   end if

   return res
end function
--


--
-- Devuelve un bigatom con la raíz n de x (x^(1/n))
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_root(object x, object exponent, integer round = 0)
   sequence res

   sequence sc = scale(SCALE + 1)
   if not bigatom(x) then
      if atom(x) and atom(exponent) then
         if x < 0 or not exponent then
            res =  NO_VALUE
         elsif exponent > 0 then
            if exponent = 2 then
               res = ba_sqrt(x)  -- muy habitual y mucho más rápida
            else
               res = ba_power(x, 1 / exponent)
            end if
         else
            res = ba_divide(ONE, ba_power(x, 1 / -exponent))
         end if
         scale(sc)

         return ba_multiply(res, 1, round)
      end if
      x = ba_new(x)
   end if

   if not bigatom(exponent) then
      if atom(exponent) and exponent = 2 then
         res = ba_sqrt(x)
         scale(sc)

         return ba_multiply(res, 1, round)
      end if
      exponent = ba_new(exponent)
   end if

   if x[SIGN] != SG_PLUS or
      exponent[SIGN] = SG_NOVALUE or
      not exponent[SIGN]
   then
      res = NO_VALUE
   elsif exponent[SIGN] = SG_PLUS then
      res = ba_power(x, ba_divide(ONE, exponent))
   else
      exponent[SIGN] = -exponent[SIGN]
      res = ba_divide(ONE, ba_power(x, ba_divide(ONE, exponent)))
   end if
   scale(sc)

   return ba_multiply(res, 1, round)
end function
--


--
-- Devuelve un bigatom con el logaritmo en base b de x
--
-- Admite atoms, bigatoms y representación de números en una cadena
--
export function ba_logb(object x, object base = 10, integer round = 0)
   if not bigatom(x) then
      x = ba_new(x)
   end if
   if not bigatom(base) then
      base = ba_new(base)
   end if
   base[SIGN] = SG_PLUS
   base = ba_ceil(base)

   return ba_divide(ba_log(x), ba_log(base), round)
end function
--


-- -------------------------------------------------------------------------------------------


-- ###############################################################
-- ###   FUNCIONES DE CONVERSION Y  MANIPULACIÓN DE BIGATOMS   ###
-- ###############################################################


--
-- Devuelve una cadena con la representación del número formateado.
--
-- sign:  si es distinto de cero también muestra el signo en los números positivos.
--        si es cero el signo solo se muestra en los números negativos.
--
-- decs:  es el número de decimales a mostrar (con el último decimal redondeado).
--        Si es negativo se muestra el número con todos sus decimales.
--        Si no, se rellena de ceros u otro carácter hasta completar la cantidad de
--        decimales pedida o se recorta y redondea al último decimal.
--
-- zsign: poner signo al cero cuando el número es menor que los límites de representación
--        pero su valor no es cero en la escala actual.
--
-- all:   indica si deben mostrarse también decimales en todos los números. No se muestran
--        si es cero. Si es uno le añade el punto decimal y los ceros necesarios hasta
--        completar el número de decimales y, si es mayor que uno se utiliza el carácter
--        DFCHAR en lugar de ceros.
--
function to_string(sequence N, integer sign = 0, integer decs = -1,
                   integer zsign = 0, integer all = 0)

   if N[SIGN] = SG_NOVALUE then
      return "<no_value>"
   end if

   N = expand(N)

   integer sg, exponent
   sequence digits
   {sg, exponent, digits} = N

   -- ajustar al número de decimales solicitado
   integer decsN = length(digits) - exponent - 1
   if decs < 0 then
      decs = decsN
   end if

   integer last = exponent + 1 + decs
   if decs < decsN then
      if digits[last + 1] >= 5 then
         -- redondea el último decimal
         digits = digits_add(digits[1..last], {1})
         exponent += digits[1]
         digits = digits[2]
      else
         digits = digits[1..last]
      end if
   elsif all or decsN then
      if all = 1 then
         digits &= repeat(0, decs - decsN)
      elsif all then
         digits &= repeat(DFCHAR - '0', decs - decsN)
      end if
   end if

   if compare(digits, repeat(0, length(digits) - exponent)) then
      zsign = 0
   end if

   digits += '0'

   -- poner punto decimal
   integer dp = exponent + 2
   if dp <= length(digits) then
      digits = insert(digits, '.', dp)
   end if

   -- poner signo
   if sg = SG_MINUS then
      digits = prepend(digits, SMINUS)
   elsif sg and (sign or zsign) then
      digits = prepend(digits, SPLUS)
   else   --  if sign then
      digits = prepend(digits, SZERO)
   end if

   return digits
end function
--


--
-- Devuelve un atom con el valor de un bigatom
--
export function bigatom_to_atom(bigatom N)
   atom val = 0
   integer sg = N[SIGN]
   sequence digits = N[DIGITS]
   if length(digits) then
      atom div = 10
      val = digits[1]
      for i = 2 to length(digits) do
         val += digits[i] / div
         div *= 10
      end for
      val *= sg * power(10, N[EXPONENT])
   end if

   return val
end function
--


--
-- Devuelve un bigatom con el valor de un entero (integer)
--
function int_to_bigatom(integer i)
   integer sign = SG_ZERO
   if i > 0 then
      sign = SG_PLUS
   else
      sign = SG_MINUS
   end if
   sequence s = sprintf("%+d", i) - '0'
   s = remove(s, 1)

   return normalize({sign, length(s) - 1, s})
end function
--


--
-- Devuelve un atom con el valor entero de una cadena que contiene la
-- representación de un número.
--
-- La conversión se detiene al encontrar el primer carácter no válido.
-- Permite el uso de espacios y subrayados como separadores de grupos.
-- Si la cadena contiene un número con decimales devuelve el valor de
-- su parte entera, el punto o la coma detienen la conversion.
--
export function int_value(sequence str)
   if not length(str) then
      return 0
   end if

   integer digit
   integer sign = (str[1] = SPLUS) - (str[1] = SMINUS)
   if sign then
      str = remove(str, 1)
   else
      sign = 1
   end if

   atom ival = 0
   for i = 1 to length(str) do
      if not integer(str[i]) then
         exit
      end if
      digit = str[i]
      if digit >= '0' and digit <= '9' then
         ival = ival * 10 + (digit - '0')
      elsif digit != UNDERLINE and digit != SPACE then
         exit
      end if
   end for

   return sign * ival
end function
--


--
-- Normaliza un bigatom.
-- Reduce un bigatom a su forma más corta eliminando
-- los ceros superfluos y ajusta el exponente.
--
export function normalize(sequence N)
   if N[SIGN] < SG_MINUS or N[SIGN] > SG_PLUS then
      return NO_VALUE
   end if

   sequence mantissa = N[DIGITS]
   integer first = 1, last = length(mantissa)
   while last and mantissa[last] = 0 do
      last -= 1
   end while
   while first <= last and mantissa[first] = 0 do
      first += 1
   end while

   if first > last or not N[SIGN] then
      N = ZERO
   else
      N[DIGITS]  = mantissa[first..last]
      N[EXPONENT] -= first - 1
   end if

   return N
end function
--


--
-- expande un bigatom con todos sus dígitos
--
function expand(sequence N)
   sequence digits = N[DIGITS]
   integer exponent = N[EXPONENT]
   integer len = length(digits) - 1

   if exponent < 0 then
      digits = repeat(0, -exponent) & digits
      N[EXPONENT] = 0
   elsif exponent > len then
      digits &= repeat(0, exponent - len)
   end if
   N[DIGITS] = digits

   return N
end function
--


--
-- devuelve ambos números con las sequences de dígitos alineadas e igualadas
-- en tamaño añadiendo los ceros necesarios tanto a derecha como a izquierda
-- de los dígitos de ambos números.
-- Igualando los exponentes y con el mismo número de dígitos en ambas mantisas.
--
function align(sequence A, sequence B)
   integer expA, expB, offset
   sequence digsA, digsB

   {?, expA, digsA} = A
   {?, expB, digsB} = B

   -- poner ceros a la izquierda
   offset = expA - expB
   if offset > 0 then
      digsB = repeat(0, offset) & digsB
      expB += offset
   else
      digsA = repeat(0, -offset) & digsA
      expA -= offset
   end if

   -- poner ceros a la derecha
   offset = length(digsA) - length(digsB)
   if offset > 0 then
      digsB &= repeat(0, offset)
   else
      digsA &= repeat(0, -offset)
   end if

   -- elimina ceros sobrantes
   integer last = length(digsA)
   while last and digsA[last] = 0 and digsB[last] = 0 do
      last -= 1
   end while
   A[DIGITS] = digsA[1..last]
   B[DIGITS] = digsB[1..last]
   A[EXPONENT] = expA
   B[EXPONENT] = expB

   return { A, B }
end function
--


--
-- suma dos sequences de dígitos
-- devuelve una sequence de dos elementos: { acarreo, resultado }
-- si acarreo no es cero indica que se ha producido un exceso en la suma
-- y por ello la longitud del resultado se ha incrementado en 1.
--
function digits_add(sequence a, sequence b)
   integer len = length(a), lenb = length(b)
   if len < lenb then
      a = repeat(0, lenb - len) & a
      len = lenb
   elsif lenb < len then
      b = repeat(0, len - lenb) & b
   end if

   sequence result = a + b

   integer carry = 0, digit
   for i = len to 1 by -1 do
      digit = result[i] + carry
      carry = digit > 9
      if carry then
         digit -= 10
      end if
      result[i] = digit
   end for
   if carry then
      result = prepend(result, carry)
   end if

   return { carry, result }
end function
--


--
-- resta dos sequences de dígitos
-- devuelve una sequence de dos elementos: { negativo, resultado }
-- si negativo no es cero indica que el resultado es negativo.
--
function digits_sub(sequence a, sequence b)

   integer len = length(a), lenb = length(b)
   if len < lenb then
      a = repeat(0, lenb - len) & a
      len = lenb
   elsif lenb < len then
      b = repeat(0, len - lenb) & b
   end if

   sequence result = a - b

   integer digit, neg = 0
   for i = len to 1 by -1 do
      digit = result[i] - neg
      neg   = digit < 0
      if neg then
         digit += 10
      end if
      result[i] = digit
   end for
   -- si es negativo está en complemento a 10
   if neg then
      result = 9 - result              -- complemento a 9
      result = digits_add(result, {1}) -- más 1
      result = result[2]
   end if

   return { neg, result }
end function
--


--
-- multiplica dos sequences de dígitos
-- devuelve una sequence de dos elementos: { acarreo, resultado }
-- si acarreo no es cero indica que se ha producido un exceso en la suma
-- y la longitud del resultado se ha incrementado en 1.
--
function digits_multiply(sequence a, sequence b)
   integer lena = length(a), lenb = length(b)
   if lena < lenb then
      {a, lena, b, lenb} = {b, lenb, a, lena}
   end if

   sequence partial, result = {}
   integer digit, carry
   for i = lenb to 1 by -1 do
      carry = 0
      partial = a * b[i] & repeat(0, lenb - i)
      for j = lena to 1 by -1 do
         digit = partial[j] + carry
         carry = 0
         if digit > 9 then
            carry = floor(digit / 10)
            digit = remainder(digit, 10)
         end if
         partial[j] = digit
      end for
      if carry then
         partial = prepend(partial, carry)
      end if
      result  = digits_add(result, partial)
      carry  += result[1]
      result  = result[2]
   end for

   return { carry, result }
end function
--


--
-- redondea el dígito dignum de la mantisa de un bigatom
-- devuelve un bigatom con el número redondeado con un digito menos.
--
function round_digits(sequence N, integer dignum)
   sequence res = N[DIGITS]

   if dignum < 1 or dignum > length(res) then
      return N
   end if

   res = digits_add(res[1..dignum], {5})
   N[EXPONENT] += res[1]
   N[DIGITS] = res[2][1..$-1]

   return normalize(N)
end function
--


--
-- eleva un bigatom a una potencia entera (uno a uno)
--
function ipower(sequence A, integer exponent)
   sequence res
   if exponent = 0 then
      return ONE
   else
      res = A
      for i = 2 to exponent do
         res = ba_multiply(A, res)
      end for
   end if

   return res
end function
--


--
-- eleva un bigatom a una potencia entera (descompone en factores primos)
-- seguro que hay mejores modos de hacerlo, pero al menos da el resultado
-- exacto con exponentes enteros
--
function intf_power(sequence A, integer exponent)
   if exponent = 0 then
      return  ONE   -- por conveniencia 0^0 = 1 (el término neutro de la multiplicación)
   end if

   sequence res
   if A[SIGN] = SG_ZERO then
      if exponent < 0 then
         res = NO_VALUE
      else
         res = ZERO
      end if
   elsif exponent < 0 then
      res = ba_divide(ONE, intf_power(A, -exponent))
      if res[SIGN] then
         res[SIGN] = A[SIGN]
      end if
   elsif equal({1}, A[DIGITS]) then
      res = A
      res[EXPONENT] *= exponent
   else
      sequence factors = get_factors(exponent)
      res = A
      for i = 1 to length(factors) do
         res = ipower(res, factors[i])
      end for
   end if

   return res
end function
--


--
-- descompone un entero en sus factores primos
--
function get_factors(integer n)
   sequence flist = {}

   while not remainder(n, 2) do
      flist &= 2
      n /= 2
   end while

   integer factor = 3
   while n >= factor do
      while remainder(n, factor) do
         factor += 2
      end while
      flist &= factor
      n /= factor
   end while
   if not length(flist) then
      flist = {n}
   end if

   return flist
end function
--


--
-- Calcula el número de Euler (e) con precisión arbitraria
-- comprobado hasta un millón de decimales, donde miré por internet un diff y solo
-- difería el último decimal, aunque tarda un poquillo... jejeje echará humo...
-- según la máquina claro.
--
--    decs:   número de decimales de e (>=0)
--    output: si es cero devuelve una cadena, si es
--            distinto de cero, devuelve un bigatom
--
--  Algoritmo original escrito en Algol por Serge Batalov
--  Reescrito en Modula-2 y modificado por Eugene Nalimov y Pavel Zemtsov
--  Adaptado a Euphoria por Carlos J. Gómez Andreu (cargoan)
--
constant N = 4,     -- grupo
	     P = 10000  -- base (10^N)
--
export function euler(integer decs, integer output = 0)
   if decs < 0 then
      decs = 0
   end if
   integer d = decs
   if remainder(decs, N) then      -- al múltiplo de N decimales
      d += N - remainder(decs, N)
   end if

   atom size = floor(d / N) + 1
   sequence x = repeat(0, size)
   sequence y = repeat(0, size)
   sequence res = {}

   if d > 0 then
-- ------------------------------------------------------------------
-- de esto ni idea, adaptado literalmente de un ejemplo de XDS,
-- compilador de Modula y Oberon para OS/2
-- ------------------------------------------------------------------
      integer e
      atom  k,l,c
      e = size
      y[1] = P
      l = 0
      c = 1
      for i = 1 to e do
         while 1 do
            c += 1
            for j = i to e do
               l = y[j] + l * P
               y[j] = floor(l / c)
               x[j] += y[j]
               l -= c * y[j]
            end for
            if y[i] < c then
               exit
            end if
            l = 0
         end while
         l = y[i]
      end for

      l = 0
      for i = e to 1 by -1 do
         k = x[i] + l
         l = floor(k / P)
         x[i] = k - l * P
      end for
-- ------------------------------------------------------------------
--  lo que sé es que aquí x es un array de enteros con los decimales
--  en base 10^N, que son N dígitos en base 10 por cada elemento
-- ------------------------------------------------------------------
      sequence s
      integer bdigit
      for i = 1 to e do
         s = repeat(0, N)
         bdigit = x[i]
         size = N
         while size do
            s[size] = remainder(bdigit, 10)
            bdigit = floor(bdigit / 10)
            size -= 1
         end while
         res &= s
      end for
   end if

   if output then
      res = {1, 0, 2 & res[1..decs]}
   elsif decs then
      res = "2." & res[1..decs] + '0'
   else
      res = "2"
   end if

   return res
end function
--

-- -------------------------------------------------------------------------------------------
-- end bigatom.e
-- -------------------------------------------------------------------------------------------

ifdef PRUEBA then
   puts(1, "\nSabías que...\n\n")
   integer n = 2, decs = 100
   scale(decs)
   bigatom ba = ba_sqrt(n)
   printf(1, "La raíz cuadrada de %d con %d decimales es:\n\t%s\n\n",
          { n, decs, ba_sprintf("%.100aB", ba) })
   ba = ba_root(n, 3)
   ba_printf(1, "Y que su raíz cúbica es:\n\t%.100aB\n\n", ba)
   ba = ba_log(n)
   ba_printf(1, "Su logaritmo natural es:\n\t%.100aB\n\n", ba)
   ba = ba_logb(n, 10)
   ba_printf(1, "Y su logaritmo decimal es:\n\t%.100aB\n\n", ba)
   decs = 843
   scale(decs)
   ba = ba_exp("1")
   printf(1, "Y que 'e' con %d decimales es:\n\t%s\n", { decs, ba_sprint(ba) })
end ifdef