Ha' centenas de fórmulas para calcular o pi. O primeiro problema que se põe é que as máquinas de calcular e os computadores trabalham com poucos digitos decimais (cerca de 10). A primeira vez que eu calculei muitos digitos de pi usei a técnica de fazer contas à mão (excepto que fiz um programa de computador para fazer operações com milhares de digitos em vez de 10 digitos). Além disso, como usei a série mais lenta a convergir (pi/4 = 1 - 1/3 + 1/5 - ...), só consegui calcular uns 100 digitos. O problema desta técnica é que não é muito prático armazenar numeros de 1 milhão de digitos. Então, as técnicas actuais armazenam apenas os restos das operaçõees da mesma maneira que quando fazemos uma divisão à mão usamos os restos para calcular os novos digitos.
O programa está escrito em javascript, isto significa que para ver a listagem do programa basta ver o texto fonte desta página html (view source).
O método usado foi o cálculo da série: pi/2 = 1 + (1/3) + (1/3)*(2/5) + (1/3)*(2/5)*(3/7) + Na verdade o algoritmo calcula x=pi/10 para eliminar a necessidade de lidar com o ponto decimal (virgula).
O algoritmo trabalha da seguinte maneira:
Seja xx a soma dos primeiros T termos. Seja s[0] = 1/R, s[1] = (1/R)*(1/3), s[2] = (1/R)*(1/3)*(2/5), ... até ao numero de termos - onde R é um inteiro positivo divisivel por 5 (para não serem necessárias fracções).
O algoritmo imprimirá D digitos de cada vez ao multiplicar repetidamente xx por BASE=10^D (imprimindo e subtraindo a parte inteira). A parte fraccional será armazenada como uma série na qual o termo i é um inteiro multiplo de s[i]. Mais especificamente, antes de cada iteração, digamos iteração i, xx*BASE^i será igual a (digitos já impressos) + coeff[0]*s[0] + coeff[1]*s[1] + coeff[2]*s[2] + ... com cada coeff[i] > 0.
Na inicialização do programa (antes da iteração 0) precisamos de ter xx = (1/5) + (1/5)*(1/3) + (1/5)*(1/3)*(2/5) + (1/5)*(1/3)*(2/5)*(3/7) + ... (R/5)*s[0] + (R/5)*s[1] + (R/5)*s[2] + ... para deixar cada coeff[i] igual a R/5.
Em cada iteração, multiplica-se cada coeff por BASE. Já que em geral o resultado será maior que BASE, transporta-se o mais possivel do resultado para o termo anterior. Já que os termos na série original têm tamanhos diferentes, a quantidade transportada (carry) tem que ser multiplicada pela razão dos dois termos antes de ser adicionada ao termo anterior. Para tornar isto possivel em inteiros, transporta-se uma quantidade divisivel pelo denominador da razão (2*i+1 para o termo i) e deixa-se o resto como o novo coeficiente. A maneira natural de fazer isto é começar no ultimo termo e continuar para a esquerda. Quando chegamos ao termo 0, calculamos (BASE*coeff[0]+carry) e calculamos o quociente e resto quando dividido por R. O quociente é enviado para o output como a parte inteira.
Note-se que depois de cada iteração, coeff[0] está entre 0 e R-1 e coeff[i] está entre 0 e 2*i. Para isto funcionar a "parte fraccional", coeff[0]*s[0] + coeff[1]*s[1] + coeff[2]*s[2] + ... tem que ser menor que um, mas na verdade quando atribuimos a cada coeff[i] o valor maximo e calculamos o resultado ficamos com 1 + 1/R. O resultado no fim de cada iteração pode ser qualquer numero entre 0 e BASE+(BASE-1)/R. Aqui o resultado é filtrado por uma rotina de correcção.
Note-se que este algoritmo é a generalização do algoritmo para converter uma fracção duma base para outra base, a diferença é que em vez da representação ser da forma x = c[0] + c[1]/Q + c[2]/(Q*Q) + ... é da forma x = c[0]*s[0] + c[1]*s[1] + c[2]*s[2] + ...
O número de iterações determinará o número de digitos resultantes.
O número de termos abandonados depois de cada iteração deve ser ajustado de modo a que o resultado dos cálculos não seja afectado. O número de termos deve ser suficientemente grande. D deve ser suficientemente pequeno para que a multiplicação de numeros de tamanho BASE = 10^D não causar overflow. Quando se usa inteiros de 32bit deve-se udar D=4 e BASE=10000. A análise de erro mostra que para BASE=10000 podemos usar DROP=14 e o número de termos até DROP*numiters+log2(DROP*numiters)+2. Resumindo, DROP deve ser o mínimo inteiro de modo que 2^DROP>BASE. Seja p=BASE/(2^DROP) e q o menor inteiro que 2^q>1/(1-p). Então o número de termos deve ser DROP*numiters+log2(DROP*numiters)+q.
Para ver a análise de erro, inicializa-se cada coeff[term] com o seu maximo valor, 2*term. A série modificada convergirá para 0.4 - assim será fácil confirmar o resultado.