Estoy tratando de escribir un programa que calcule dígitos decimales de π a 1000 dígitos o más.

Para practicar la programación de bajo nivel por diversión, el programa final se escribirá en ensamblador, en una CPU de 8 bits que no tiene multiplicaciones ni divisiones, y solo realiza sumas de 16 bits. Para facilitar la implementación, es deseable poder usar solo operaciones con enteros sin signo de 16 bits y usar un algoritmo iterativo. La velocidad no es una preocupación importante. Y la multiplicación y división rápidas está más allá del alcance de esta pregunta, así que no consideres esos problemas también.

Antes de implementarlo en ensamblaje, todavía estoy tratando de encontrar un algoritmo utilizable en C en mi computadora de escritorio. Hasta ahora, encontré que la siguiente serie es razonablemente eficiente y relativamente fácil de implementar.

La fórmula se deriva de la serie de Leibniz utilizando una técnica de aceleración de convergencia. Para obtenerla, consulte Computing the Digits in π, de Carl D. Offner (https://cs.umb.edu/~offner/files/pi.pdf), página 19-26. La fórmula final se muestra en la página 26. La fórmula inicial que escribí tenía algunos errores tipográficos, actualice la página para ver la fórmula corregida. El término constante 2 como máximo, se explica en la página 54. El documento también describía un algoritmo iterativo avanzado, pero no lo usé aquí.

Si uno evalúa la serie usando muchos (por ejemplo, 5000) términos, es posible obtener miles de dígitos de π fácilmente, y encontré que esta serie también es fácil de evaluar iterativamente usando este algoritmo:

Algoritmo

1. Primero, reorganiza la fórmula para obtener sus términos constantes de una matriz.

2. Complete la matriz con 2 para comenzar la primera iteración, por lo tanto, la nueva fórmula se parece a la original.

3. Dejar carry = 0.

4. Empiece por el término mayor. Obtenga un término (2) de la matriz, multiplique el término por PRECISION para realizar una división de punto fijo contra 2 * i + 1y guarde el recordatorio como el nuevo término de la matriz. Luego suma el siguiente término. Ahora decrementar ipasar al siguiente término, repetir hasta i == 1. Finalmente agregue el término final x_0.

5. Debido a que se usa un entero de 16 bits, PRECISION es 10, por lo tanto se obtienen 2 dígitos decimales, pero solo el primer dígito es válido. Guarde el segundo dígito como acarreo. Mostrar el primer dígito más llevar.

6. x_0 es el entero 2, no se debe sumar para las iteraciones sucesivas, borrarlo.

7. Ir al paso 4 para calcular el siguiente dígito decimal, hasta que tengamos todos los dígitos que queremos.

Implementación 1

Traduciendo este algoritmo a C:

#include <stdio.h>
#include <stdint.h>

#define N 2160
#define PRECISION 10

uint16_t terms[N + 1] = {0};

int main(void)
{
    /* initialize the initial terms */
    for (size_t i = 0; i < N + 1; i++) {
        terms[i] = 2;
    }

    uint16_t carry = 0;
    for (size_t j = 0; j < N / 4; j++) {
        uint16_t numerator = 0;
        uint16_t denominator;
        uint16_t digit;

        for (size_t i = N; i > 0; i--) {
            numerator += terms[i] * PRECISION;
            denominator = 2 * i + 1;

            terms[i] = numerator % denominator;
            numerator /= denominator;
            numerator *= i;
        }
        numerator += terms[0] * PRECISION;
        digit = numerator / PRECISION + carry;
        carry = numerator % PRECISION;

        printf("%01u", digit);

        /* constant term 2, only needed for the first iteration. */
        terms[0] = 0;
    }
    putchar('\n');
}

El código puede calcular π con 31 dígitos decimales, hasta que comete un error.

31415926535897932384626433832794
10 <-- wrong

Algunas veces digit + carry es mayor que 9, por lo que necesita un acarreo adicional. Si tenemos mucha mala suerte, puede haber incluso un doble acarreo, un triple acarreo, etc. Usamos un ring-buffer para almacenar los últimos 4 dígitos. Si se detecta un acarreo adicional, generamos un retroceso para borrar el dígito anterior, realizamos un acarreo y los reimprimimos. Esta es solo una solución fea para la Prueba de concepto, que es irrelevante para mi pregunta sobre el desbordamiento, pero para completar, aquí está. Algo mejor sería implementado en el futuro.

Implementación 2 con acarreo repetido

#include <stdio.h>
#include <stdint.h>

#define N 2160
#define PRECISION 10

#define BUF_SIZE 4

uint16_t terms[N + 1] = {0};

int main(void)
{
    /* initialize the initial terms */
    for (size_t i = 0; i < N + 1; i++) {
        terms[i] = 2;
    }

    uint16_t carry = 0;
    uint16_t digit[BUF_SIZE];
    int8_t idx = 0;

    for (size_t j = 0; j < N / 4; j++) {
        uint16_t numerator = 0;
        uint16_t denominator;

        for (size_t i = N; i > 0; i--) {
            numerator += terms[i] * PRECISION;
            denominator = 2 * i + 1;

            terms[i] = numerator % denominator;
            numerator /= denominator;
            numerator *= i;
        }
        numerator += terms[0] * PRECISION;
        digit[idx] = numerator / PRECISION + carry;

        /* over 9, needs at least one carry op. */
        if (digit[idx] > 9) {
            for (int i = 1; i <= 4; i++) {
                if (i > 3) {
                    /* allow up to 3 consecutive carry ops */
                    fprintf(stderr, "ERROR: too many carry ops!\n");
                    return 1;
                }
                /* erase a digit */
                putchar('\b');

                /* carry */
                digit[idx] -= 10;
                idx--;
                if (idx < 0) {
                    idx = BUF_SIZE - 1;
                }
                digit[idx]++;            
                if (digit[idx] < 10) {
                    /* done! reprint the digits */
                    for (int j = 0; j <= i; j++) {
                        printf("%01u", digit[idx]);
                        idx++;
                        if (idx > BUF_SIZE - 1) {
                            idx = 0;
                        }
                    }
                    break;
                }
            }
        }
        else {
            printf("%01u", digit[idx]);
        }

        carry = numerator % PRECISION;
        terms[0] = 0;

        /* put an element to the ring buffer */
        idx++;
        if (idx > BUF_SIZE - 1) {
            idx = 0;
        }
    }
    putchar('\n');
}

Genial, ahora el programa puede calcular correctamente 534 dígitos de π, hasta que comete un error.

3141592653589793238462643383279502884
1971693993751058209749445923078164062
8620899862803482534211706798214808651
3282306647093844609550582231725359408
1284811174502841027019385211055596446
2294895493038196442881097566593344612
8475648233786783165271201909145648566
9234603486104543266482133936072602491
4127372458700660631558817488152092096
2829254091715364367892590360011330530
5488204665213841469519415116094330572
7036575959195309218611738193261179310
5118548074462379962749567351885752724
8912279381830119491298336733624406566
43086021394946395
22421 <-- wrong

Desbordamiento de enteros de 16 bits

Resulta que, durante el cálculo de los términos más grandes al principio, el término de error se vuelve bastante grande, ya que los divisores al principio están en el rango de ~4000. Al evaluar la serie, numerator en realidad comienza a desbordarse en la multiplicación inmediatamente.

El desbordamiento de enteros es insignificante a la hora de calcular los primeros 500 dígitos, pero va empeorando cada vez más, hasta dar un resultado incorrecto.

Cambiando uint16_t numerator = 0 para uint32_t numerator = 0 puede resolver este problema y calcular π a más de 1000 dígitos.

Sin embargo, como mencioné antes, mi plataforma de destino es una CPU de 8 bits y solo tiene operaciones de 16 bits. ¿Hay algún truco para resolver el problema de desbordamiento de enteros de 16 bits que estoy viendo aquí? usando solo uno o más uint16_t? Si no es posible evitar la aritmética de precisión múltiple, ¿cuál es el método más simple para implementarla aquí? Sé que de alguna manera necesito introducir una “palabra de extensión” adicional de 16 bits, pero no estoy seguro de cómo puedo implementarla.

Y gracias de antemano por su paciencia para comprender el largo contexto aquí.

Eche un vistazo al control de calidad relacionado:

• Desafío Baking-Pi – Entender y mejorar

está usando Wiki: Bailey-Borwein-Plouffe_formula que es más adecuado para la aritmética de enteros.

Sin embargo, el verdadero desafío sería:

• ¿Cómo convierto un número binario muy largo a decimal?

Como probablemente quieras imprimir el número en base dec…

Además, si necesita llevar un lenguaje de nivel más alto que asm, eche un vistazo a esto:

• No puedo hacer que el valor se propague a través del carry

Puede modificarlo para manejar tantos bits de acarreo como necesite (si aún es menor que el ancho de bits del tipo de datos).

[Edit1] Ejemplo de BBP en C++/VCL

Utilicé esta fórmula (tomada de la página Wiki vinculada anteriormente):

convertido a punto fijo…

//---------------------------------------------------------------------------
AnsiString str_hex2dec(const AnsiString &hex)
    {
    char c;
    AnsiString dec="",s;
    int i,j,l,ll,cy,val;
    int  i0,i1,i2,i3,sig;
    sig=+1; l=hex.Length();
    if (l) { c=hex[l]; if (c=='h') l--; if (c=='H') l--; }
    i0=0; i1=l; i2=0; i3=l;
    for (i=1;i<=l;i++)      // scan for parts of number
        {
        char c=hex[i];
        if (c=='-') sig=-sig;
        if ((c=='.')||(c==',')) i1=i-1;
        if ((c>='0')&&(c<='9')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        if ((c>='A')&&(c<='F')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        if ((c>='a')&&(c<='f')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        }

    l=0; s=""; if (i0) for (i=i0;i<=i1;i++)
        {
        c=hex[i];
             if ((c>='0')&&(c<='9')) c-='0';
        else if ((c>='A')&&(c<='F')) c-='A'-10;
        else if ((c>='a')&&(c<='f')) c-='A'-10;
        for (cy=c,j=1;j<=l;j++)
            {
            val=(s[j]<<4)+cy;
            s[j]=val%10;
            cy  =val/10;
            }
        while (cy>0)
            {
            l++;
            s+=char(cy%10);
            cy/=10;
            }
        }
    if (s!="")
        {
        for (j=1;j<=l;j++) { c=s[j]; if (c<10) c+='0'; else c+='A'-10; s[j]=c; }
        for (i=l,j=1;j<i;j++,i--) { c=s[i]; s[i]=s[j]; s[j]=c; }
        dec+=s;
        }
    if (dec=="") dec="0";
    if (sig<0) dec="-"+dec;

    if (i2)
        {
        dec+='.';
        s=hex.SubString(i2,i3-i2+1);
        l=s.Length();
        for (i=1;i<=l;i++)
            {
            c=s[i];
                 if ((c>='0')&&(c<='9')) c-='0';
            else if ((c>='A')&&(c<='F')) c-='A'-10;
            else if ((c>='a')&&(c<='f')) c-='A'-10;
            s[i]=c;
            }
        ll=((l*1234)>>10);  // num of decimals to compute
        for (cy=0,i=1;i<=ll;i++)
            {
            for (cy=0,j=l;j>=1;j--)
                {
                val=s[j];
                val*=10;
                val+=cy;
                s[j]=val&15;
                cy=val>>4;
                }
            dec+=char(cy+'0');
            for (;;)
                {
                if (!l) break;;
                if (s[l]) break;
                l--;
                }
            if (!l) break;;
            }
        }

    return dec;
    }
//---------------------------------------------------------------------------
AnsiString pi_BBP() // https://en.wikipedia.org/wiki/Bailey–Borwein–Plouffe_formula
    {
    const int N=100;        // 32*N bit uint arithmetics
    int sh;
    AnsiString s;
    uint<N> pi,a,b,k,k2,k3,k4;

    for (pi=0,sh=(N<<5)-8,k=0;sh>=0;k++,sh-=4)
        {
        k2=k*k;
        k3=k2*k;
        k4=k3*k;
        a =k2* 120;
        a+=k * 151;
        a+=     47;
        b =k4* 512;
        b+=k3*1024;
        b+=k2* 712;
        b+=k * 194;
        b+=     15;
        a<<=sh;
        pi+=a/b;
        }
    pi<<=4;
    s=pi.strhex();
    s=s.Insert(".",2);
    return str_hex2dec(s);
    }
//---------------------------------------------------------------------------

El código está usando VCL AnsiString que es una cadena autoasignable y mía uint<N> plantilla que es aritmética de enteros sin signo de 32*N ancho de bits basado en el mío ALU32. Como puede ver, solo necesita sumas y multiplicaciones de divisiones de enteros grandes para esto (todas las demás cosas se pueden hacer en enteros normales).

Aquí el resultado decádico frente a la referencia Pi de 1000 dígitos:

ref: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
BPP: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778048187

El valor bigint calculado se exporta a una cadena hexadecimal y luego se convierte a una base decádica usando str_hex2dec del enlace de arriba. El número de iteraciones depende del ancho de bits de destino.

El código aún no está optimizado…

¿Qué pasa con la implementación de la aritmética de 32 bits?

Para agregar, agregue las dos palabras de orden superior (16 bits), luego las dos palabras de orden inferior, pruebe el bit de desbordamiento y lleve al resultado de orden superior si es necesario.

Si puede predecir cuándo ocurrirá el desbordamiento, puede cambiar de aritmética de 16 a 32 bits cuando sea necesario.

La prueba del bit de desbordamiento no se puede hacer en C puro, requerirá algún ensamblaje en línea o una función intrínseca.

De lo contrario, puede inspirarse en esta respuesta: https://codereview.stackexchange.com/a/37178/39646

Hay un truco:

Considere usar una matriz para los numeradores y otra matriz para los denominadores. Cada posición representaría la cantidad de veces que ese número se multiplica para obtener el número real.

Un ejemplo:

(1 * 2 * 3 * 7 * 7) / (3 * 6 * 8)

Se representaría como:

num[] = {1, 1, 1, 0, 0, 0, 2};
denom[] = {0, 0, 1, 0, 0, 1, 0, 1};

Luego considere factorizar en números primos cada número antes de almacenarlo, para que tenga números más bajos. Ahora necesitará otra matriz para almacenar todos los números primos:

primes[] = {2, 3, 5, 7};
num[] = {1, 1, 0, 2};
denom[] = {4, 2, 0, 0};

Esto le permitirá almacenar números inimaginablemente grandes, pero tarde o temprano querrá volver a transformarlos en números, por lo que primero querrá simplificar esto. La forma de hacerlo es simplemente restar factors[i] += num[i] - denom[i] para cada campo en las matrices, para cada fracción en la serie. Querrá simplificar después de cada iteración, para minimizar el riesgo de desbordamiento.

factors[] = {-3, -1, 0, 2};

Cuando necesite el número, simplemente hágalo num *= pow(primes[i], factors[i]); si el factor es positivo, o num /= pow(primes, -factors[i]); si es negativo, para cada campo de las matrices. (No hacer nada si es 0.

num y denom son matrices temporales que se utilizan para almacenar una fracción, la matriz donde se almacena el resultado es factors. Recuerda memset las matrices temporales antes de cada uso.

Esta explicación es útil para cualquier fracción grande. Para adaptarlo a su problema específico, es posible que necesite usar una función de potencia entera y también multiplicar por 10^algo para convertir la parte decimal en una parte integral. Esa es tu misión, si la aceptas 🙂