Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

MEJORANDO LA PRECISION NUMERICA DE LOS RESULTADOS DE SIMULACIONES DIGITALES CONTINUAS USANDO ARITMETICA DE PRECISION ARBITRARIA Graciela Bueno1, Enrique Arjona1, Andrés Arjona2 RESUMEN En corridas de simulación y particularmente cuando se utilizan modelos digitales continuos, la propagación del error numérico resultante de la ejecución de procesos iterativos puede afectar considerablemente la exactitud de los resultados de la simulación. Algunos lenguajes de programación, como C y Pascal, atenúan la propagación del error a través de permitir el uso de un número fijo relativamente grande de dígitos significativos. Sin embargo, investigar si un número de dígitos significativos es adecuado para un modelo particular no es una tarea fácil para el usuario promedio. Además, si el número máximo de dígitos significativos permitidos por el lenguaje no es satisfactorio, no existe una solución directa para este problema. En este artículo, para atenuar la propagación del error se propone un conjunto de rutinas sencillas para aritmética de precisión arbitraria. La implementación de estas rutinas en programas existentes es directa y una 1

Especialidad de Estadística, Colegio de Postgraduados, Km. 36.5 Carr. México Texcoco, Montecillo 56230 Edo. de México 2 Universidad de las Américas A.C., Puebla Col. Roma, México D.F.

63

Graciela Bueno, Enrique Arjona y Andrés Arjona

vez que la implementación se concluye, probar que un número de dígitos significativos particular es adecuado para el modelo es muy fácil. Palabras clave: Simulación digital continua, Error de redondeo, Aritmética de precisión arbitraria INTRODUCCION En corridas de simulación, y particularmente al usar modelos digitales continuos, el cálculo de variables endógenas y de salida requiere de la retroalimentación de valores intermedios dentro de procesos iterativos que son ejecutados muchas veces. La propagación del error numérico resultante de dichos procesos puede afectar considerablemente la exactitud de los resultados de la simulación. Los errores numéricos pueden estar inherentes en los parámetros y en las relaciones aproximadas usadas en las ecuaciones de los modelos para representar sistemas reales. Estos errores son independientes de la precisión del lenguaje o la computadora en la que los modelos son implementados y la única manera de evitarlos es utilizando parámetros y ecuaciones más precisas. Además de los errores numéricos inherentes existen dos tipos de errores que son más bien independientes de la precisión de las ecuaciones usadas en los modelos. 64

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

El primer tipo de error es producido por la naturaleza aproximada de las rutinas de integración numéricas usadas para extrapolar valores simulados. Este tipo de error puede ser evitado o atenuado usando integración simbólica siempre que sea posible, o rutinas de integración elaboradas y tamaños de paso pequeños (Maron y López 1991). El segundo tipo de error, que es en el que estamos interesados, es producido por la precisión numérica fija usada para la representación interna de los números en una computadora digital. Este tipo de error es conocido como error de redondeo y se puede propagar. Algunos lenguajes de programación, como C y Pascal, atenúan la propagación del error a través de permitir el uso de un número fijo, relativamente largo, de dígitos significativos. Sin embargo, investigar si un número de dígitos significativos es adecuado para un modelo particular no es tarea fácil para el usuario promedio. Además, si el número máximo de dígitos significativos permitidos por el lenguaje no es satisfactorio, no existe una solución directa para este problema. La propagación del error debida a la precisión numérica fija no es exclusiva de los modelos continuos digitales y ocurre en muchos algoritmos matemáticos. Con el paso de los años se han propuesto diversas soluciones a este problema. Las propuestas incluyen el uso de aritmética entera de precisión ilimitada (Collins, 1966), aritmética 65

Graciela Bueno, Enrique Arjona y Andrés Arjona

racional (Knuth 1969) y aritmética entera de precisión arbitraria (Knuth 1969). Las maneras en las que estos sistemas aritméticos pueden ser implementados en una computadora, son usualmente considerados complejos para el usuario promedio, por lo que no se incluyen en libros de simulación, estadística u optimización. Por consiguiente, el uso de sistemas aritméticos especializados ha sido el dominio de paquetes muy especializados de programación. Otro hecho que posiblemente ha retrasado la generalización del uso de rutinas de atenuación del error de redondeo, es que los algoritmos modificados con esas rutinas consumen más tiempo que los algoritmos originales. El tiempo de cómputo era un recurso caro hasta la llegada de computadoras personales poderosas lo cual ocurrió hace apenas algunos años. Los sistemas computacionales algebraicos como Mathematica (Wolfram 1991), que permiten el uso de aritmética racional y

de

precisión arbitraria, parecen ser una solución a mano para evitar la propagación del error. Sin embargo, el uso de dichos sistemas para la implementación de modelos de simulación continuos no ha sido generalizado aún. Los usuarios usan mayormente otras herramientas como los lenguajes C, CSL, y SIMULINK (Bennet 1995). El uso de rutinas de precisión arbitraria no es necesariamente una tarea compleja. Su complejidad está relacionada íntimamente al tipo de 66

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

rutinas utilizadas y la manera en que son implementadas y sus implementaciones no siempre requieren de la ayuda de un especialista en ciencias de la computación. En este artículo, un conjunto de rutinas aritméticas de precisión arbitraria es propuesto. Estas rutinas pueden proveer una solución a la propagación del error porque el número de dígitos significativos usados para la representación numérica en un modelo puede incrementarse tanto como sea necesario. La implementación de estas rutinas en programas existentes es directa y requiere básicamente del cambio de expresiones aritméticas de notación infija a postfija. Una vez que la implementación de las rutinas en un modelo ha sido realizada, probar que un número particular de dígitos es adecuado para un modelo es muy fácil y consiste en comparar los resultados usando diferentes precisiones que difieren en un pequeño número de dígitos. El artículo está organizado de la siguiente manera: Primero, se presenta un ejemplo de propagación del error y se dan generalidades acerca de la aritmética racional y de precisión arbitraria. Después, las rutinas de precisión arbitraria propuestas son descritas acompañadas junto con el proceso para su implementación en programas existentes. Finalmente una discusión es dada.

67

Graciela Bueno, Enrique Arjona y Andrés Arjona

PROPAGACION DEL ERROR DE REDONDEO, ARITMETICA RACIONAL Y PRECISION ARBITRARIA Como se ha mencionado, el error de redondeo es producido por la precisión numérica fija usada para la representación interna de números en una computadora digital. Este ocurre cuando el número de dígitos significativos disponible para la representación interna de números en una computadora, es menor que el número de dígitos significativos obtenidos del resultado de alguna operación aritmética. Una vez que ha ocurrido un error de

redondeo, puede propagarse en operaciones

aritméticas subsecuentes. Por ejemplo, si tenemos una computadora que internamente permite sólo ocho dígitos significativos y realizamos la siguiente suma: 1234567 + .1234567 = 1234567.1234567 el resultado contiene catorce dígitos significativos, aún cuando cada uno de los números sumados contienen sólo siete dígitos significativos. En la computadora, el resultado será almacenado como 1234567.1 . Eso es, se han perdido seis dígitos significativos en esta operación. Para ver como este error de redondeo puede propagarse a más dígitos significativos, realizaremos dos operaciones más. Primero, con el resultado obtenido arriba realizamos una simple resta: 68

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

1234567.1 - 1234565.1 = 2 y finalmente, con este nuevo resultado, realizamos la siguiente multiplicación: 2 x 7000 = 14000 El resultado exacto que podríamos haber obtenido al realizar las operaciones descritas, con un lápiz y papel es 14164.1969. Lo que significa, que en las operaciones arriba mostradas, a pesar de que se habían perdido inicialmente sólo seis dígitos, el resultado final es truncado en siete dígitos significativos debido a la propagación del error. El ejemplo anterior involucra sólo números con un número finito de dígitos significativos. Sin embargo, en aritmética real (aritmética con números reales), se pueden obtener resultados con un número infinito de dígitos significativos. Por ejemplo: 10 : 3 = 3.333333... En estos casos, obviamente, el número de dígitos significativos perdidos en cualquier representación decimal interna es también infinita.

69

Graciela Bueno, Enrique Arjona y Andrés Arjona

Los números que tienen una representación infinita en un sistema numérico posicional como el decimal, pueden tener una representación finita en otro sistema numérico posicional y viceversa. Por ejemplo, si usamos el sistema ternario (base 3), el número mencionado tendrá la representación finita 10.1 porque en el sistema ternario: 10.1 = 1 x 3**1 + 0 x 3**0 + 1 x 3**(-1) = 3.33... en base 10 y el número 10.111... , que tiene una representación infinita en el sistema ternario, es equivalente al número 3.5 en el sistema decimal. En cualquier sistema numérico posicional, habrá un conjunto infinito de números que para su representación requieren de un número infinito de dígitos significativos. Por lo que el cambio de base no es una solución al problema de error por redondeo. Para superar el problema de propagación del error de redondeo se ha propuesto el uso de aritmética racional y de precisión arbitraria. El uso de aritmética racional en una computadora es soportada por el hecho de que todos los números decimales finitos y periódicos como 3.3333... y 23.67813813813... , son números racionales y tienen una representación finita como el cociente de dos enteros. La aritmética 70

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

racional siempre producirá números racionales y por lo tanto, susceptibles a representación finita. Por consiguiente, en el sistema de números racional, los únicos números reales que no tienen una representación finita y que son susceptibles al error de redondeo son los números irracionales como π, e, √2, etc. Las reglas para transformar números decimales finitos y periódicos a sus representaciones racionales equivalentes son muy sencillas. Un número decimal finito puede ser expresado como un número racional cuyo denominador es una potencia de 10, por ejemplo: 14164.1969 = 141641969/10000 Un número decimal periódico puede ser expresado como la suma infinita de números racionales, que en turno puede ser expresada como la suma de dos números racionales y por consiguiente como un número racional único, por ejemplo: 23.67813813... = 23.67 + .00813 + .00000813 + ... = 2367/100 + 813/10000 + 813/100000000 + ... = 2367/100 + 813/99900 = 2365446/99900

71

Graciela Bueno, Enrique Arjona y Andrés Arjona

Mas aún, los números racionales pueden ser reducidos a números racionales equivalentes usando el algoritmo de Euclídes para calcular máximos comunes divisores. Los números reducidos pueden contener menos dígitos que los números originales. Por ejemplo: 2365446/99900 = 394241/16650 A pesar de su simplicidad, la aritmética racional no es fácil de implementar en una computadora. Las representaciones de números racionales como cociente de dos enteros pueden involucrar números muy largos, por lo que realizar operaciones con ellos puede requerir precisión infinita de aritmética de enteros. Además, las rutinas para la reducción de números racionales calculados deben de ser muy eficientes para minimizar el tiempo de ejecución. La aritmética de precisión arbitraria permite elegir un número arbitrario fijo máximo de dígitos significativos para la representación interna de números en una computadora. Usualmente, el número de dígitos para la precisión elegida depende de la exactitud requerida en los resultados de una aplicación en particular. A mayor precisión seleccionada, mayor será el tiempo que las aplicaciones requerirán para su ejecución.

72

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

La precisión arbitraria es más fácil de implementar que la aritmética racional. Sin embargo, su utilidad para atenuar la propagación del error de redondeo es más débil que la utilidad de la aritmética racional. Con precisión arbitraria, todos los números decimales infinitos almacenados en una computadora, racionales e irracionales, tendrán error de redondeo y podrían causar la propagación del error. Más aún, números decimales finitos pueden también tener error de redondeo dependiendo de su número de dígitos significativos. Sin embargo, en comparación con la máxima precisión permitida en muchos lenguajes de programación, la cual es del orden de 20 dígitos significativos, el uso de precisión arbitraria de 50 o más dígitos, puede llegar a ser muy satisfactoria. La implementación de precisión arbitraria en una computadora puede realizarse de diversas formas. Desde unas muy eficientes, usando lenguaje máquina y optimizando la ejecución y la cantidad de memoria de la computadora. A unos muy simples usando un lenguaje de alto nivel y un byte para almacenar cada dígito decimal. Las implementaciones simples son más fáciles de rastrear y requieren de menos programación, pero usualmente son lentas. Implementaciones completas requieren de una librería extensa de funciones matemáticas. Las rutinas presentadas en la siguiente sección contienen características eficientes, pero pueden considerarse simples. 73

Graciela Bueno, Enrique Arjona y Andrés Arjona

UN CONJUNTO DE RUTINAS DE PRECISION ARBITRARIA Las rutinas aquí presentadas han sido implementadas como procedimientos en Pascal, pero pueden fácilmente ser cambiadas a funciones o subrutinas en otro lenguaje de alto nivel como C o FORTRAN. El almacenamiento interno de números y operaciones aritméticas es realizado en base 256, por lo que cada dígito requiere un byte. En los procedimientos, AP significa precisión arbitraria; EP precisión externa; e IP precisión interna. AP (un tipo); EP (una constante); e IP (también una constante) son declaradas globales en el programa principal. Las rutinas han sido clasificadas en cuatro grupos. El primer grupo consiste de los procedimientos READLONG y WRITELONG usadas para entrada y salida. El segundo grupo está integrado por las funciones GREATERLONG, LESSLONG, y EQULONG usadas para pruebas de relación. El tercer grupo se refiere a los procedimientos CONSTLONG, ADDLONG, SUBLONG, MULLONG, y DIVLONG usadas para operaciones aritméticas. El último grupo consiste de los procedimientos de soporte INITLONG, CONVLONG, NORMLONG y ROUNDLONG. Las rutinas en este último grupo no son llamadas por el usuario. A continuación se presentan las rutinas.

74

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

Grupo 1 READLONG PROCEDURE READLONG(VAR APVAR:AP); { Este procedimiento lee un número decimal de punto fijo de precisión arbitraria desde el } {teclado y almacena su representación de precisión arbitraria de base BASE en APVAR } var inputstr:string[EP+2]; begin readln(inputstr); CONVLONG(inputstr, APVAR); end; WRITELONG PROCEDURE WRITELONG(APVAR:AP); { Este procedimiento despliega en la pantalla la representación decimal de punto fijo de} {número de precisión arbitraria arbitraria en base BASE almacenada en APVAR } var appart,digit,ten,aptemp:AP; ndigits,i,exp:integer; outputstr:string[EP+2]; decdigit:string[1]; begin INITLONG(10,ten); outputstr:=''; exp:=APVAR[IP+1]-127; ndigits:=0; if APVAR[1]=0 then begin outputstr:='0'; writeln(outputstr); exit; end; if exp>0 then begin appart:=APVAR; appart[0]:=1; for i:=exp+1 to IP do appart[i]:=0; repeat DIVLONG(appart,ten,aptemp); for i:=aptemp[IP+1]-127+1 to IP do aptemp[i]:=0; MULLONG (aptemp,ten,aptemp); SUBLONG(appart,aptemp,digit); if digit[IP+1]>127 then begin str(digit[1]:1,decdigit); outputstr:=decdigit+outputstr; end 75

Graciela Bueno, Enrique Arjona y Andrés Arjona

else outputstr:='0'+outputstr; ndigits:=ndigits+1; DIVLONG(aptemp,ten,aptemp);appart:=aptemp; until (appart[IP+1]=127); exp:=APVAR[IP+1]-127; i:=exp+1; while (APVAR[i]=0) and (i0 then begin i:=exp+1; while (iOP2 regresa} {un valor verdadero(true) } var pos:byte; begin GREATERLONG:=true; if OP1[1]=0 then if (OP2[0]=1) or (OP2[1]=0) then begin GREATERLONG:=false; exit; end else exit; if OP2[1]=0 then if OP1[0]=0 then begin GREATERLONG:=false; exit; end else exit; if (OP1[0]=0) and (OP2[0]=1) then begin GREATERLONG:=false; exit; end; if (OP1[0]=1) and (OP2[0]=0) then exit; if OP1[0]=1 then begin if OP1[IP+1]OP2[IP+1] then exit; end; if OP1[0]=0 then begin if OP1[IP+1]>OP2[IP+1] then begin GREATERLONG:=false; exit; end else if OP1[IP+1]= length(inputstr) then exit; inputstr:=copy(inputstr,point+1,length(inputstr)-point); INITLONG(1,powers); DIVLONG(powers,ten,powers); for pos:=1 to length(inputstr) do if (inputstr[pos] in ['0'..'9']) then begin INITLONG(0,digit); digit[IP+1]:=128; digit[1]:=ord(inputstr[pos])-48; MULLONG(powers,digit,digit); ADDLONG(right,digit,right); DIVLONG(powers,ten,powers);end else begin writeln('arbitrary precision invalid input data'); halt; end; ADDLONG(LEFT,RIGHT,APVAR); APVAR[0]:=sign; end; NORMLONG PROCEDURE NORMLONG(VAR UNNORMED); { Este procedimiento normaliza un número doble de precisión arbitraria APVAR } var acum:array[0..2*IP+1] of longint absolute UNNORMED; i,j,limit,twoip:integer; begin twoip:=2*IP; limit:=acum[twoip+1]; if twoipIP then begin writeln('arbitrary precision overflow');halt;end; if acum[1]=0 then begin for i:=2 to twoip do if acum[i]0 then begin writeln('arbitrary precision underflow'); halt; end; acum[0]:=1; acum[twoip+1]:=127;end; end; TRUNCLONG PROCEDURE TRUNCLONG(VAR UNROUND; VAR APVAR:AP); { Este procedimiento trunca un número doble de precisión arbitraria UNROUND y } { almacena el resultado en el número de precisión arbitraria APVAR } var i,pos:integer; acum:array[0..2*IP+1] of longint absolute UNROUND; begin APVAR[IP+1]:=acum[2*IP+1]; if ROUND=0 then begin for pos:=0 to IP do APVAR[pos]:=acum[pos]; exit; end else begin APVAR[0]:=acum[0]; if acum[IP+1]>=trunc(BASE/2) then acum[IP]:=acum[IP]+1; for pos:=IP downto 2 do if acum[pos]=BASE then begin APVAR[pos]:=0; acum[pos-1]:=acum[pos-1]+1; end else APVAR[pos]:=acum[pos]; if acum[1]=BASE then begin for pos:=IP downto 3 do APVAR[pos]:=APVAR[pos-1]; APVAR[1]:=1; APVAR[2]:=0; APVAR[IP+1]:=APVAR[IP+1]+1;end else APVAR[1]:=acum[1]; end; 85

Graciela Bueno, Enrique Arjona y Andrés Arjona

for i:=0 to IP*2+1 do acum[i]:=0; end;

IMPLEMENTACION DE PRECISION ARBITRARIA EN PROGRAMAS EXISTENTES La implementación de las rutinas mencionadas en programas existentes es directa. Toda la aritmética real en sus programas es transformada en aritmética de precisión arbitraria. Las reglas para la implementación son las siguientes: 1. Incluya en su programa principal (main) la declaración global de BASE, ROUND, EP, IP y AP. BASE es la base utilizada internamente por las rutinas y puede ser cualquier número desde 10 a 256. Las bases que son múltiplos de 10 usualmente se desempeñan más rápidamente que otras bases. EP es la precisión deseada (máximo número de dígitos significativos). El máximo valor permitido por EP en las rutinas dadas es 256, sin embargo, cuando los números son introducidos por el teclado, el número máximo para EP será 130. IP es la precisión interna equivalente en bytes de EP y depende de la base usada, para base 10 IP=EP; para base 256 IP=1/2EP. Para otras bases use valores intermedios para IP y revise la adecuidad. ROUND es utilizado para indicar si se desea redondeo interno, un valor de 1 significa un redondeo al siguiente dígito 86

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

significativo, un valor de 0 significa truncar sin redondear. AP es la declaración del tipo de precisión arbitraria. Las declaraciones a incluir en el programa principal, antes de la primera declaración BEGIN son similares a las siguientes: CONST EP=100; IP=50; BASE=250; ROUND=1; TYPE AP=ARRAY[0..IP+1] OF BYTE; 2. Incluya en su programa principal y también antes de la primera declaración begin, los encabezados de todas las rutinas con la opción forward, y después del grupo de encabezados, el código de todas las rutinas. Por ejemplo, el encabezado del procedimiento READLONG es: PROCEDURE READLONG(VAR APVAR:AP); forward; 3. Incluya en su programa principal la declaración de las variables de precisión arbitraria que serán utilizadas para almacenar resultados intermedios de operaciones aritméticas (paso 5). Usualmente, dos variables serán suficientes, pero use más si así lo desea. Las variables temporales podrán ser declaradas de la siguiente manera: VAR TEMP1, TEMP2: AP;

87

Graciela Bueno, Enrique Arjona y Andrés Arjona

si su programa ya usa variables con esos nombres, use otros nombres para las variables temporales. 4. Cambie todas las variables de tipo real utilizadas en su programa al tipo AP. Por ejemplo, si en su programa existe una declaración como la siguiente: VAR A,B,C: REAL; cámbiela a: VAR A,B,C: AP; 5. Cambie todas las declaraciones de lectura (escritura) de datos reales en

su

programa

a

llamadas

al

procedimiento

READLONG

(WRITELONG). Estos procedimientos sólo leen (escriben) una variable a la vez, por lo que una declaración puede requerir más de una llamada a READLONG (WRITELONG). Por ejemplo, si en su programa hay una declaración como la siguiente (A y B son variables reales): READLN (A,B); cámbiela a: READLONG(A); READLONG(B); 6. Cambie todas las expresiones aritméticas que involucren datos reales en su programa a llamadas a los procedimientos aritméticos 88

Comunicaciones en Socioeconomía, Estadística e Informática 2002, Vol. 6 Núm. 1, pp 63-91

ADDLONG, SUBLONG, MULLONG y DIVLONG. Note que esos procedimientos realizan operaciones sólo entre dos números, por lo que puede ser necesario reorganizar una expresión en su programa en varias llamadas a los procedimientos mencionados. Los resultados intermedios son almacenados en variables temporales de precisión arbitraria (véase paso 2). Si una expresión contiene constantes reales, éstas deberán ser transformadas a precisión arbitraria utilizando el procedimiento CONSTLONG. Por ejemplo, si su programa contiene una expresión como la siguiente (A, B, C y D son variables reales): A:= (B+2.0)/(C*D) cámbiela a: CONSTLONG(2.0,TEMP1); ADDLONG(B,TEMP1,TEMP1); MULTLONG(C,D,TEMP2); DIVLONG(TEMP1,TEMP2,A) 7. Cambie todas las expresiones relacionales que involucren datos reales en su programa, a llamadas a los procedimientos GREATERLONG, LESSLONG, y EQULONG. Por ejemplo, si su programa contiene la siguiente declaración (A, B y C son variables reales): IF A