Calculando Pi na Plaquinha (Uno, PiPico, ESP32, BlackPill, etc)

Uma coisa que está na minha lista faz algum tempo é fazer testes de velocidade com as várias plaquinhas que eu tenho. Mais pela brincadeira, pois estes testes normalmente não representam direito o desempenho em aplicações reais. E aí me surgiu a ideia de fazer a comparação calculando os dígitos do Pi.

Chovendo no molhado, o Pi é a razão entre a circunferência de um círculo e o seu diâmetro. É um número "irracional", que não pode ser colocado na forma de fração e na forma decimal tem um número infinito de dígitos. Calcular os dígitos do Pi tem ocupado pessoas desde muito tempo atrás (a wikipedia tem uma cronologia).

A partir de 1400 sabe-se que o Pi pode ser calculado através da soma de séries infinitas. Com o advento do computador foram desenvolvidos algorítimos que facilitam e aceleram o cálculo. Na minha busca eu rapidamente achei o seguinte programa (de autoria de Dik T Winter):

view plaincopy to clipboardprint?

1. int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;  
2. for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,  
3. f[b]=d%--g,d/=g--,--b;d*=b);}  

int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,
f[b]=d%--g,d/=g--,--b;d*=b);}

Você talvez encontre uma versão com pequenas diferenças, esta acima calcula os primeiros 800 dígitos do Pi. Este programa está, propositalmente, compactado e ofuscado. Para os meus testes eu fiz uma versão um pouco mais legível e adaptada para compilação na IDE do Arduino:

view plaincopy to clipboardprint?

1. /* 
2.  * Calculo dos dígito de Pi 
3.  * Usando o algoritimo Spigot de Rabinowitz e Wagon 
4.  * Adaptação da versão compacta e obsfucada escrita por Dik T. Winter: 
5.  *  
6.  * int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5; 
7.  * for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a, 
8.  * f[b]=d%--g,d/=g--,--b;d*=b);} 
9.  *  
10.  * Comentários nas declarações copiados de 
11.  * https://stackoverflow.com/questions/4084571/implementing-the-spigot-algorithm-for-%CF%80-pi 
12.  *  
13.  * Daniel Quadros junho/2021 
14.  */  
15.    
16. #define NDIGITS 10000          //max digits to compute  
17. #define LEN (NDIGITS/4+1)*14   //nec. array length  
18.    
19. int32_t a = 10000;             //new base, 4 decimal digits  
20. int32_t b = 0;                 //nominator prev. base  
21. int32_t c = LEN;               //index  
22. int32_t d = 0;                 //accumulator and carry  
23. int32_t e = 0;                 //save previous 4 digits  
24. int32_t f[LEN+1];              //array of 4-digit-decimals  
25. int32_t g = 0;                 //denom previous base  
26.   
27. void setup() {  
28.   delay(1000);  
29.   Serial.begin(115200);  
30. }  
31.   
32. // Não usado  
33. void loop() {  
34.   while (Serial.read() != '\n') {  
35.     delay(100);   
36.   }  
37.   Serial.println("PI = ");  
38.   long inicio = millis();  
39.   calcula();  
40.   long tempo = millis() - inicio;  
41.   Serial.println();  
42.   Serial.print ("Tempo = ");  
43.   Serial.print (tempo);  
44.   Serial.println ("ms");  
45. }  
46.   
47. // Cálculo dos dígitos do Pi  
48. void calcula() {  
49.   char dig[5] = "0000"; // para fazer o print  
50.   int n = 0;            // para mudar de linha a cada 100 dígitos  
51.   
52.   c = LEN;  
53.   for(b = 0; b < c; b++) {  
54.     f[b] = a/5;  
55.   }  
56.   
57.   e = 0;  
58.   for (; c > 0; c -= 14) {  
59.       d = 0;  
60.       g = c*2;  
61.       b = c;  
62.       for (;;) {  
63.           d += f[b]*a;  
64.           f[b] = d % --g;  
65.           d /= g--;  
66.           if (--b == 0) {  
67.               break;  
68.           }  
69.           d *= b;  
70.       }  
71.       uint16_t val = e+d/a;  
72.       dig[0] = (val / 1000) + '0';   
73.       dig[1] = ((val / 100) % 10) + '0';   
74.       dig[2] = ((val / 10) % 10) + '0';   
75.       dig[3] = (val % 10) + '0';   
76.       Serial.print (dig);  
77.       n += 4;  
78.       if (n == 100) {  
79.         Serial.println();  
80.         n = 0;  
81.       }  
82.       e = d % a;  
83.   }  
84. }  

/*
 * Calculo dos dígito de Pi
 * Usando o algoritimo Spigot de Rabinowitz e Wagon
 * Adaptação da versão compacta e obsfucada escrita por Dik T. Winter:
 * 
 * int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;
 * for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,
 * f[b]=d%--g,d/=g--,--b;d*=b);}
 * 
 * Comentários nas declarações copiados de
 * https://stackoverflow.com/questions/4084571/implementing-the-spigot-algorithm-for-%CF%80-pi
 * 
 * Daniel Quadros junho/2021
 */
 
#define NDIGITS 10000          //max digits to compute
#define LEN (NDIGITS/4+1)*14   //nec. array length
 
int32_t a = 10000;             //new base, 4 decimal digits
int32_t b = 0;                 //nominator prev. base
int32_t c = LEN;               //index
int32_t d = 0;                 //accumulator and carry
int32_t e = 0;                 //save previous 4 digits
int32_t f[LEN+1];              //array of 4-digit-decimals
int32_t g = 0;                 //denom previous base

void setup() {
  delay(1000);
  Serial.begin(115200);
}

// Não usado
void loop() {
  while (Serial.read() != '\n') {
    delay(100); 
  }
  Serial.println("PI = ");
  long inicio = millis();
  calcula();
  long tempo = millis() - inicio;
  Serial.println();
  Serial.print ("Tempo = ");
  Serial.print (tempo);
  Serial.println ("ms");
}

// Cálculo dos dígitos do Pi
void calcula() {
  char dig[5] = "0000"; // para fazer o print
  int n = 0;            // para mudar de linha a cada 100 dígitos

  c = LEN;
  for(b = 0; b < c; b++) {
    f[b] = a/5;
  }

  e = 0;
  for (; c > 0; c -= 14) {
      d = 0;
      g = c*2;
      b = c;
      for (;;) {
          d += f[b]*a;
          f[b] = d % --g;
          d /= g--;
          if (--b == 0) {
              break;
          }
          d *= b;
      }
      uint16_t val = e+d/a;
      dig[0] = (val / 1000) + '0'; 
      dig[1] = ((val / 100) % 10) + '0'; 
      dig[2] = ((val / 10) % 10) + '0'; 
      dig[3] = (val % 10) + '0'; 
      Serial.print (dig);
      n += 4;
      if (n == 100) {
        Serial.println();
        n = 0;
      }
      e = d % a;
  }
}

Mas o que está sendo feito, e porque? Não me arrisco a dizer que entendi completamente, mas aqui vão algumas informações e referências. O algoritmo usado pertence a uma classe chamada "spigot" (torneira): os dígitos vão sendo gerados na ordem e não são usados diretamente no cálculo dos dígitos seguintes. Este algorítimo específico foi criado em 1995 por  Stan Wagon and Stanley Rabinowitz.

Algumas características interessantes deste algorítimo:

• Requer somente soma, subtração, multiplicação e divisão de inteiros
• Os dígitos são gerados de 4 em 4
• A quantidade de dígitos a calcular precisa ser definida a priori
• Requer um vetor com 14 inteiros para cada conjunto de 4 dígitos

Este último ponto é importante quando falamos em microcontroladores com pouca memória. Para rodar no Arduino Uno (que tem modestos 2K) eu reduzi o número de dígitos para 200 e usei inteiros sem sinal de 16 bits no vetor. Aliás, não encontrei nenhum comentário sobre o tamanho dos inteiros utilizados no algorítimo, pelas minhas experiências dá para calcular 200 dígitos usando inteiro sem sinal de 16 bits no vetor (e inteiros de 32 bits nas demais variáveis). Usando inteiros de 32 bits em todas as variáveis calculou corretamente 20000 dígitos.

Algumas referências sobre o algorítimo:

Explicação (incompleta) e a versão ofuscada: https://crypto.stanford.edu/pbc/notes/pi/code.html

Artigo que inclui uma explicação mais completa: https://www.cs.ox.ac.uk/people/jeremy.gibbons/publications/spigot.pdf

Fonte para os comentários nas variáveis: https://stackoverflow.com/questions/4084571/implementing-the-spigot-algorithm-for-%cf%80-pi/25837097#25837097

Eu rodei o programa (acertando o número de dígitos conforme a memória disponível) em algumas placas da minha coleção: Arduino Uno, Raspberry Pi Pico, ESP32, Black Pill (STM32) e Seeduino XIAO. O resultado esta aqui, no meu canal do YouTube.