¿Cuál es la forma más rápida de obtener el valor de π?
Estoy buscando la manera más rápida de obtener el valor de π, como desafío personal. Más específicamente, estoy usando formas que no implican el uso #definede constantes comoM_PI o codificar el número.
El siguiente programa prueba las diversas formas que conozco. La versión ensamblada en línea es, en teoría, la opción más rápida, aunque claramente no es portátil. Lo he incluido como punto de referencia para compararlo con las otras versiones. En mis pruebas, con funciones integradas, la 4 * atan(1)versión es más rápida en GCC 4.2, porque se pliega automáticamente atan(1)en una constante. Con -fno-builtinlo especificado, elatan2(0, -1) versión es más rápida.
Aquí está el programa de prueba principal ( pitimes.c):
#include <math.h>
#include <stdio.h>
#include <time.h>
#define ITERS 10000000
#define TESTWITH(x) { \
diff = 0.0; \
time1 = clock(); \
for (i = 0; i < ITERS; ++i) \
diff += (x) - M_PI; \
time2 = clock(); \
printf("%s\t=> %e, time => %f\n", #x, diff, diffclock(time2, time1)); \
}
static inline double
diffclock(clock_t time1, clock_t time0)
{
return (double) (time1 - time0) / CLOCKS_PER_SEC;
}
int
main()
{
int i;
clock_t time1, time2;
double diff;
/* Warmup. The atan2 case catches GCC's atan folding (which would
* optimise the ``4 * atan(1) - M_PI'' to a no-op), if -fno-builtin
* is not used. */
TESTWITH(4 * atan(1))
TESTWITH(4 * atan2(1, 1))
#if defined(__GNUC__) && (defined(__i386__) || defined(__amd64__))
extern double fldpi();
TESTWITH(fldpi())
#endif
/* Actual tests start here. */
TESTWITH(atan2(0, -1))
TESTWITH(acos(-1))
TESTWITH(2 * asin(1))
TESTWITH(4 * atan2(1, 1))
TESTWITH(4 * atan(1))
return 0;
}
Y el material de ensamblaje en línea ( fldpi.c) que solo funcionará para sistemas x86 y x64:
double
fldpi()
{
double pi;
asm("fldpi" : "=t" (pi));
return pi;
}
Y un script de compilación que crea todas las configuraciones que estoy probando ( build.sh):
#!/bin/sh
gcc -O3 -Wall -c -m32 -o fldpi-32.o fldpi.c
gcc -O3 -Wall -c -m64 -o fldpi-64.o fldpi.c
gcc -O3 -Wall -ffast-math -m32 -o pitimes1-32 pitimes.c fldpi-32.o
gcc -O3 -Wall -m32 -o pitimes2-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -fno-builtin -m32 -o pitimes3-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -ffast-math -m64 -o pitimes1-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -m64 -o pitimes2-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -fno-builtin -m64 -o pitimes3-64 pitimes.c fldpi-64.o -lm
Además de realizar pruebas entre varios indicadores del compilador (también he comparado 32 bits con 64 bits porque las optimizaciones son diferentes), también intenté cambiar el orden de las pruebas. Pero aún así, la atan2(0, -1)versión siempre sale victoriosa.
El método Monte Carlo , como se mencionó, aplica algunos conceptos excelentes pero, claramente, no es el más rápido, ni de lejos, ni de ninguna manera razonable. Además, todo depende del tipo de precisión que busques. El π más rápido que conozco es el que tiene los dígitos codificados. En cuanto a Pi y Pi[PDF] , hay muchas fórmulas.
Aquí hay un método que converge rápidamente: alrededor de 14 dígitos por iteración. PiFast , la aplicación más rápida actualmente, utiliza esta fórmula con la FFT. Simplemente escribiré la fórmula, ya que el código es sencillo. Esta fórmula casi fue encontrada por Ramanujan y descubierta por Chudnovsky . En realidad, así es como calculó varios miles de millones de dígitos del número, por lo que no es un método que deba ignorarse. La fórmula se desbordará rápidamente y, dado que estamos dividiendo factoriales, sería ventajoso retrasar dichos cálculos para eliminar términos.
dónde,
A continuación se muestra el algoritmo Brent-Salamina . Wikipedia menciona que cuando a y b están "lo suficientemente cerca", entonces (a + b)² / 4t será una aproximación de π. No estoy seguro de lo que significa "lo suficientemente cerca", pero según mis pruebas, una iteración obtuvo 2 dígitos, dos obtuvieron 7 y tres tuvieron 15; por supuesto, esto es con dobles, por lo que podría tener un error basado en su representación y el verdadero cálculo podría ser más preciso.
let pi_2 iters =
let rec loop_ a b t p i =
if i = 0 then a,b,t,p
else
let a_n = (a +. b) /. 2.0
and b_n = sqrt (a*.b)
and p_n = 2.0 *. p in
let t_n = t -. (p *. (a -. a_n) *. (a -. a_n)) in
loop_ a_n b_n t_n p_n (i - 1)
in
let a,b,t,p = loop_ (1.0) (1.0 /. (sqrt 2.0)) (1.0/.4.0) (1.0) iters in
(a +. b) *. (a +. b) /. (4.0 *. t)
Por último, ¿qué tal un pi golf (800 dígitos)? ¡160 caracteres!
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);}
Realmente me gusta este programa, porque aproxima π mirando su propia área.
IOCCC 1988: westley.c
#define _ -F<00||--F-OO--; int F=00,OO=00;main(){F_OO();printf("%1.3f\n",4.*-F/OO/OO);}F_OO() { _-_-_-_ _-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_-_-_-_-_ _-_-_-_-_-_-_-_ _-_-_-_ }
