Comparación de velocidad con Project Euler: C vs Python vs Erlang vs Haskell
performance (18)
Tomé el Problema # 12 del Proyecto Euler como un ejercicio de programación y comparé mis implementaciones (seguramente no óptimas) en C, Python, Erlang y Haskell. Para obtener un mayor tiempo de ejecución, busco el primer número de triángulo con más de 1000 divisores en lugar de 500 como se indica en el problema original.
El resultado es el siguiente:
DO:
lorenzo@enzo:~/erlang$ gcc -lm -o euler12.bin euler12.c
lorenzo@enzo:~/erlang$ time ./euler12.bin
842161320
real 0m11.074s
user 0m11.070s
sys 0m0.000s
Pitón:
lorenzo@enzo:~/erlang$ time ./euler12.py
842161320
real 1m16.632s
user 1m16.370s
sys 0m0.250s
Python con PyPy:
lorenzo@enzo:~/Downloads/pypy-c-jit-43780-b590cf6de419-linux64/bin$ time ./pypy /home/lorenzo/erlang/euler12.py
842161320
real 0m13.082s
user 0m13.050s
sys 0m0.020s
Erlang:
lorenzo@enzo:~/erlang$ erlc euler12.erl
lorenzo@enzo:~/erlang$ time erl -s euler12 solve
Erlang R13B03 (erts-5.7.4) [source] [64-bit] [smp:4:4] [rq:4] [async-threads:0] [hipe] [kernel-poll:false]
Eshell V5.7.4 (abort with ^G)
1> 842161320
real 0m48.259s
user 0m48.070s
sys 0m0.020s
Haskell:
lorenzo@enzo:~/erlang$ ghc euler12.hs -o euler12.hsx
[1 of 1] Compiling Main ( euler12.hs, euler12.o )
Linking euler12.hsx ...
lorenzo@enzo:~/erlang$ time ./euler12.hsx
842161320
real 2m37.326s
user 2m37.240s
sys 0m0.080s
Resumen:
- C: 100%
- Python: 692% (118% con PyPy)
- Erlang: 436% (135% gracias a RichardC)
- Haskell: 1421%
Supongo que C tiene una gran ventaja, ya que utiliza largos para los cálculos y no enteros de longitud arbitraria como los otros tres. Además, no es necesario cargar primero un tiempo de ejecución (¿los otros?).
Pregunta 1: ¿Erlang, Python y Haskell pierden velocidad debido al uso de enteros de longitud arbitraria o no, siempre que los valores sean menores que MAXINT
?
Pregunta 2: ¿Por qué Haskell es tan lento? ¿Hay un indicador de compilador que apague los frenos o es mi implementación? (Lo último es bastante probable ya que Haskell es un libro con siete sellos para mí).
Pregunta 3: ¿Me puede ofrecer algunos consejos sobre cómo optimizar estas implementaciones sin cambiar la forma en que determiné los factores? Optimización de cualquier manera: mejor, más rápido, más "nativo" al idioma.
EDITAR:
Pregunta 4: ¿Mis implementaciones funcionales permiten LCO (optimización de la última llamada, también conocida como eliminación de recursión de cola) y por lo tanto evitan agregar marcos innecesarios a la pila de llamadas?
Realmente intenté implementar el mismo algoritmo lo más similar posible en los cuatro idiomas, aunque debo admitir que mi conocimiento de Haskell y Erlang es muy limitado.
Códigos fuente utilizados:
#include <stdio.h>
#include <math.h>
int factorCount (long n)
{
double square = sqrt (n);
int isquare = (int) square;
int count = isquare == square ? -1 : 0;
long candidate;
for (candidate = 1; candidate <= isquare; candidate ++)
if (0 == n % candidate) count += 2;
return count;
}
int main ()
{
long triangle = 1;
int index = 1;
while (factorCount (triangle) < 1001)
{
index ++;
triangle += index;
}
printf ("%ld/n", triangle);
}
#! /usr/bin/env python3.2
import math
def factorCount (n):
square = math.sqrt (n)
isquare = int (square)
count = -1 if isquare == square else 0
for candidate in range (1, isquare + 1):
if not n % candidate: count += 2
return count
triangle = 1
index = 1
while factorCount (triangle) < 1001:
index += 1
triangle += index
print (triangle)
-module (euler12).
-compile (export_all).
factorCount (Number) -> factorCount (Number, math:sqrt (Number), 1, 0).
factorCount (_, Sqrt, Candidate, Count) when Candidate > Sqrt -> Count;
factorCount (_, Sqrt, Candidate, Count) when Candidate == Sqrt -> Count + 1;
factorCount (Number, Sqrt, Candidate, Count) ->
case Number rem Candidate of
0 -> factorCount (Number, Sqrt, Candidate + 1, Count + 2);
_ -> factorCount (Number, Sqrt, Candidate + 1, Count)
end.
nextTriangle (Index, Triangle) ->
Count = factorCount (Triangle),
if
Count > 1000 -> Triangle;
true -> nextTriangle (Index + 1, Triangle + Index + 1)
end.
solve () ->
io:format ("~p~n", [nextTriangle (1, 1) ] ),
halt (0).
factorCount number = factorCount'' number isquare 1 0 - (fromEnum $ square == fromIntegral isquare)
where square = sqrt $ fromIntegral number
isquare = floor square
factorCount'' number sqrt candidate count
| fromIntegral candidate > sqrt = count
| number `mod` candidate == 0 = factorCount'' number sqrt (candidate + 1) (count + 2)
| otherwise = factorCount'' number sqrt (candidate + 1) count
nextTriangle index triangle
| factorCount triangle > 1000 = triangle
| otherwise = nextTriangle (index + 1) (triangle + index + 1)
main = print $ nextTriangle 1 1
Pregunta 1: ¿Erlang, python y haskell pierden velocidad debido al uso de enteros de longitud arbitraria o no, siempre que los valores sean menores que MAXINT?
Esto es poco probable. No puedo decir mucho sobre Erlang y Haskell (bueno, quizás un poco sobre Haskell a continuación) pero puedo señalar muchos otros cuellos de botella en Python. Cada vez que el programa intenta ejecutar una operación con algunos valores en Python, debe verificar si los valores son del tipo adecuado y cuesta un poco de tiempo. Su función factorCount
simplemente asigna una lista con range (1, isquare + 1)
varias veces, y la asignación en memoria de estilo, malloc
-styled es mucho más lenta que iterar en un rango con un contador como lo hace en C. Notablemente, el factorCount()
Se llama varias veces y así se asigna una gran cantidad de listas. Además, no olvidemos que Python se interpreta y el intérprete de CPython no tiene un gran enfoque en la optimización.
EDIT : oh, bueno, observo que está utilizando Python 3, por lo que range()
no devuelve una lista, sino un generador. En este caso, mi punto sobre la asignación de listas es medio incorrecto: la función simplemente asigna objetos de range
, que sin embargo son ineficientes pero no tan ineficientes como asignar una lista con una gran cantidad de elementos.
Pregunta 2: ¿Por qué haskell es tan lento? ¿Hay un indicador de compilador que apague los frenos o es mi implementación? (Lo último es bastante probable, ya que haskell es un libro con siete sellos para mí).
¿Estás usando Hugs ? Los abrazos son un intérprete considerablemente lento. Si lo estás usando, tal vez puedas obtener un mejor momento con GHC , pero solo estoy pensando en la hipotesis, el tipo de cosas que un buen compilador de Haskell hace bajo el capó es bastante fascinante y mucho más allá de mi comprensión :)
Pregunta 3: ¿Me puede ofrecer algunos consejos sobre cómo optimizar estas implementaciones sin cambiar la forma en que determiné los factores? Optimización de cualquier manera: mejor, más rápido, más "nativo" al idioma.
Yo diría que estás jugando un juego sin gracia. Lo mejor de saber varios idiomas es usarlos de la manera más diferente posible :) Pero estoy divagando, simplemente no tengo ninguna recomendación para este punto. Lo siento, espero que alguien pueda ayudarte en este caso :)
Pregunta 4: ¿Mis implementaciones funcionales permiten LCO y, por lo tanto, evitan agregar marcos innecesarios a la pila de llamadas?
Por lo que recuerdo, solo debe asegurarse de que su llamada recursiva sea el último comando antes de devolver un valor. En otras palabras, una función como la que se muestra a continuación podría usar dicha optimización:
def factorial(n, acc=1):
if n > 1:
acc = acc * n
n = n - 1
return factorial(n, acc)
else:
return acc
Sin embargo, no tendría tal optimización si su función fuera como la de abajo, porque hay una operación (multiplicación) después de la llamada recursiva:
def factorial2(n):
if n > 1:
f = factorial2(n-1)
return f*n
else:
return 1
Separé las operaciones en algunas variables locales para dejar claro qué operaciones se ejecutan. Sin embargo, lo más habitual es ver estas funciones como se muestra a continuación, pero son equivalentes para el punto que estoy señalando:
def factorial(n, acc=1):
if n > 1:
return factorial(n-1, acc*n)
else:
return acc
def factorial2(n):
if n > 1:
return n*factorial(n-1)
else:
return 1
Tenga en cuenta que le corresponde al compilador / intérprete decidir si hará la recursión de la cola. Por ejemplo, el intérprete de Python no lo hace si recuerdo bien (en mi ejemplo usé Python solo por su sintaxis fluida). De todos modos, si encuentras cosas extrañas como funciones factoriales con dos parámetros (y uno de los parámetros tiene nombres como acc
, accumulator
, etc.), ahora sabes por qué la gente lo hace :)
Pregunta 3: ¿Me puede ofrecer algunos consejos sobre cómo optimizar estas implementaciones sin cambiar la forma en que determiné los factores? Optimización de cualquier manera: mejor, más rápido, más "nativo" al idioma.
La implementación de C es subóptima (como lo sugiere Thomas M. DuBuisson), la versión usa enteros de 64 bits (es decir, tipo de datos largo ). Investigaré la lista de ensamblajes más adelante, pero con una conjetura informada, hay algunos accesos de memoria en el código compilado, lo que hace que el uso de enteros de 64 bits sea significativamente más lento. Es ese o el código generado (ya que el hecho de que pueda incluir menos ints de 64 bits en un registro SSE o redondear un doble a un entero de 64 bits es más lento).
Aquí está el código modificado (simplemente reemplace por largo con int y explícitamente en línea factorCount, aunque no creo que esto sea necesario con gcc -O3):
#include <stdio.h>
#include <math.h>
static inline int factorCount(int n)
{
double square = sqrt (n);
int isquare = (int)square;
int count = isquare == square ? -1 : 0;
int candidate;
for (candidate = 1; candidate <= isquare; candidate ++)
if (0 == n % candidate) count += 2;
return count;
}
int main ()
{
int triangle = 1;
int index = 1;
while (factorCount (triangle) < 1001)
{
index++;
triangle += index;
}
printf ("%d/n", triangle);
}
Corriendo + cronometrando da:
$ gcc -O3 -lm -o euler12 euler12.c; time ./euler12
842161320
./euler12 2.95s user 0.00s system 99% cpu 2.956 total
Para referencia, la implementación de haskell por Thomas en la respuesta anterior da:
$ ghc -O2 -fllvm -fforce-recomp euler12.hs; time ./euler12 [9:40]
[1 of 1] Compiling Main ( euler12.hs, euler12.o )
Linking euler12 ...
842161320
./euler12 9.43s user 0.13s system 99% cpu 9.602 total
Conclusión: no quitarle nada a ghc, es un gran compilador, pero gcc normalmente genera un código más rápido.
Pregunta 1: ¿Erlang, Python y Haskell pierden velocidad debido al uso de enteros de longitud arbitraria o no, siempre que los valores sean menores que MAXINT?
La pregunta uno puede ser respondida negativamente para Erlang. La última pregunta se responde utilizando Erlang adecuadamente, como en:
http://bredsaal.dk/learning-erlang-using-projecteuler-net
Ya que es más rápido que su ejemplo inicial de C, supongo que hay numerosos problemas, ya que otros ya lo han cubierto en detalle.
Este módulo de Erlang se ejecuta en una netbook barata en aproximadamente 5 segundos ... Utiliza el modelo de hilos de red en erlang y, como tal, demuestra cómo aprovechar el modelo de eventos. Podría ser distribuido en muchos nodos. Y es rapido No es mi código.
-module(p12dist).
-author("Jannich Brendle, [email protected], http://blog.bredsaal.dk").
-compile(export_all).
server() ->
server(1).
server(Number) ->
receive {getwork, Worker_PID} -> Worker_PID ! {work,Number,Number+100},
server(Number+101);
{result,T} -> io:format("The result is: /~w./~n", [T]);
_ -> server(Number)
end.
worker(Server_PID) ->
Server_PID ! {getwork, self()},
receive {work,Start,End} -> solve(Start,End,Server_PID)
end,
worker(Server_PID).
start() ->
Server_PID = spawn(p12dist, server, []),
spawn(p12dist, worker, [Server_PID]),
spawn(p12dist, worker, [Server_PID]),
spawn(p12dist, worker, [Server_PID]),
spawn(p12dist, worker, [Server_PID]).
solve(N,End,_) when N =:= End -> no_solution;
solve(N,End,Server_PID) ->
T=round(N*(N+1)/2),
case (divisor(T,round(math:sqrt(T))) > 500) of
true ->
Server_PID ! {result,T};
false ->
solve(N+1,End,Server_PID)
end.
divisors(N) ->
divisor(N,round(math:sqrt(N))).
divisor(_,0) -> 1;
divisor(N,I) ->
case (N rem I) =:= 0 of
true ->
2+divisor(N,I-1);
false ->
divisor(N,I-1)
end.
La siguiente prueba se realizó en una: CPU Intel (R) Atom (TM) N270 a 1,60 GHz
~$ time erl -noshell -s p12dist start
The result is: 76576500.
^C
BREAK: (a)bort (c)ontinue (p)roc info (i)nfo (l)oaded
(v)ersion (k)ill (D)b-tables (d)istribution
a
real 0m5.510s
user 0m5.836s
sys 0m0.152s
Con Haskell, realmente no necesitas pensar explícitamente en las recursiones.
factorCount number = foldr factorCount'' 0 [1..isquare] -
(fromEnum $ square == fromIntegral isquare)
where
square = sqrt $ fromIntegral number
isquare = floor square
factorCount'' candidate
| number `rem` candidate == 0 = (2 +)
| otherwise = id
triangles :: [Int]
triangles = scanl1 (+) [1,2..]
main = print . head $ dropWhile ((< 1001) . factorCount) triangles
En el código anterior, he reemplazado las recursiones explícitas en la respuesta de @Thomas con operaciones de lista comunes. El código sigue haciendo exactamente lo mismo sin que nos preocupe la recursión de la cola. Funciona (~ 7.49s ) aproximadamente un 6% más lento que la versión en la respuesta de @Thomas (~ 7.04s ) en mi máquina con GHC 7.6.2, mientras que la versión C de @Raedwulf ejecuta ~ 3.15s . Parece que GHC ha mejorado a lo largo del año.
PD. Sé que es una pregunta antigua, y me topé con ella a partir de búsquedas en Google (olvidé lo que estaba buscando, ahora ...). Solo quería hacer un comentario sobre la pregunta sobre LCO y expresar mis sentimientos sobre Haskell en general. Quería comentar sobre la respuesta principal, pero los comentarios no permiten bloques de código.
Echa un vistazo a este blog . Durante el año pasado, más o menos, ha solucionado algunos de los problemas del Proyecto Euler en Haskell y Python, y en general ha encontrado que Haskell es mucho más rápido. Creo que entre esos idiomas tiene más que ver con tu fluidez y estilo de codificación.
Cuando se trata de la velocidad de Python, ¡estás usando la implementación incorrecta! Prueba PyPy , y para cosas como esta, encontrarás que es mucho más rápido.
En lo que respecta a la optimización de Python, además de usar PyPy (para aceleraciones bastante impresionantes con cero cambios en su código), puede usar la cadena de herramientas de traducción de PyPy para compilar una versión compatible con RPython, o Cython para construir un módulo de extensión, ambos que son más rápidas que la versión C en mis pruebas, con el módulo Cython casi el doble de rápido . Para referencia, también incluyo los resultados de referencia de C y PyPy:
C (compilado con gcc -O3 -lm
)
% time ./euler12-c
842161320
./euler12-c 11.95s
user 0.00s
system 99%
cpu 11.959 total
PyPy 1.5
% time pypy euler12.py
842161320
pypy euler12.py
16.44s user
0.01s system
99% cpu 16.449 total
RPython (usando la última revisión de PyPy, c2f583445aee
)
% time ./euler12-rpython-c
842161320
./euler12-rpy-c
10.54s user 0.00s
system 99%
cpu 10.540 total
Cython 0.15
% time python euler12-cython.py
842161320
python euler12-cython.py
6.27s user 0.00s
system 99%
cpu 6.274 total
La versión RPython tiene un par de cambios clave. Para traducir en un programa independiente, debe definir su target
, que en este caso es la función main
. Se espera que acepte sys.argv
como único argumento, y se requiere que devuelva un int. Puedes traducirlo usando translate.py, % translate.py euler12-rpython.py
que se traduce a C y lo compila por ti.
# euler12-rpython.py
import math, sys
def factorCount(n):
square = math.sqrt(n)
isquare = int(square)
count = -1 if isquare == square else 0
for candidate in xrange(1, isquare + 1):
if not n % candidate: count += 2
return count
def main(argv):
triangle = 1
index = 1
while factorCount(triangle) < 1001:
index += 1
triangle += index
print triangle
return 0
if __name__ == ''__main__'':
main(sys.argv)
def target(*args):
return main, None
La versión de Cython fue reescrita como un módulo de extensión _euler12.pyx
, el cual importo y llamo desde un archivo normal de Python. El _euler12.pyx
es esencialmente el mismo que su versión, con algunas declaraciones de tipo estático adicionales. Setup.py tiene la plantilla normal para construir la extensión, utilizando python setup.py build_ext --inplace
.
# _euler12.pyx
from libc.math cimport sqrt
cdef int factorCount(int n):
cdef int candidate, isquare, count
cdef double square
square = sqrt(n)
isquare = int(square)
count = -1 if isquare == square else 0
for candidate in range(1, isquare + 1):
if not n % candidate: count += 2
return count
cpdef main():
cdef int triangle = 1, index = 1
while factorCount(triangle) < 1001:
index += 1
triangle += index
print triangle
# euler12-cython.py
import _euler12
_euler12.main()
# setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
ext_modules = [Extension("_euler12", ["_euler12.pyx"])]
setup(
name = ''Euler12-Cython'',
cmdclass = {''build_ext'': build_ext},
ext_modules = ext_modules
)
Honestamente, tengo muy poca experiencia con RPython o Cython, y me sorprendió gratamente los resultados. Si está utilizando CPython, escribir los bits de código que requieren un uso intensivo de la CPU en un módulo de extensión Cython parece ser una forma realmente fácil de optimizar su programa.
Hay algunos problemas con la implementación de Erlang. Como punto de referencia para lo siguiente, mi tiempo de ejecución medido para su programa Erlang no modificado fue de 47.6 segundos, en comparación con 12.7 segundos para el código C.
Lo primero que debe hacer si desea ejecutar código de Erlang intensivo computacionalmente es usar código nativo. La erlc +native euler12
con erlc +native euler12
el tiempo a 41.3 segundos. Sin embargo, esta es una aceleración mucho menor (solo un 15%) de lo que se espera de la compilación nativa en este tipo de código, y el problema es su uso de -compile(export_all)
. Esto es útil para la experimentación, pero el hecho de que todas las funciones sean potencialmente accesibles desde el exterior hace que el compilador nativo sea muy conservador. (El emulador BEAM normal no se ve tan afectado). Reemplazar esta declaración con -export([solve/0]).
proporciona una aceleración mucho mejor: 31,5 segundos (casi el 35% de la línea de base).
Pero el código en sí tiene un problema: para cada iteración en el bucle factorCount, realiza esta prueba:
factorCount (_, Sqrt, Candidate, Count) when Candidate == Sqrt -> Count + 1;
El código C no hace esto. En general, puede ser difícil hacer una comparación equitativa entre diferentes implementaciones del mismo código, y en particular si el algoritmo es numérico, porque necesita estar seguro de que realmente están haciendo lo mismo. Un ligero error de redondeo en una implementación debido a algún encasillado en algún lugar puede hacer que haga muchas más iteraciones que la otra, aunque ambos alcancen el mismo resultado.
Para eliminar esta posible fuente de error (y deshacerme de la prueba adicional en cada iteración), reescribí la función factorCount de la siguiente manera, modelada en detalle en el código C:
factorCount (N) ->
Sqrt = math:sqrt (N),
ISqrt = trunc(Sqrt),
if ISqrt == Sqrt -> factorCount (N, ISqrt, 1, -1);
true -> factorCount (N, ISqrt, 1, 0)
end.
factorCount (_N, ISqrt, Candidate, Count) when Candidate > ISqrt -> Count;
factorCount ( N, ISqrt, Candidate, Count) ->
case N rem Candidate of
0 -> factorCount (N, ISqrt, Candidate + 1, Count + 2);
_ -> factorCount (N, ISqrt, Candidate + 1, Count)
end.
Esta reescritura, no export_all
y compilación nativa, me dio el siguiente tiempo de ejecución:
$ erlc +native euler12.erl
$ time erl -noshell -s euler12 solve
842161320
real 0m19.468s
user 0m19.450s
sys 0m0.010s
que no es tan malo en comparación con el código C:
$ time ./a.out
842161320
real 0m12.755s
user 0m12.730s
sys 0m0.020s
considerando que Erlang no está en absoluto orientado a escribir código numérico, ser solo un 50% más lento que C en un programa como este es bastante bueno.
Finalmente, con respecto a sus preguntas:
Pregunta 1: ¿Erlang, python y haskell pierden velocidad debido al uso de enteros de longitud arbitraria o no, siempre que los valores sean menores que MAXINT?
Sí, un poco. En Erlang, no hay manera de decir "usar aritmética de 32/64 bits con envolvente", así que a menos que el compilador pueda probar algunos límites en sus enteros (y generalmente no puede), debe verificar todos los cálculos para ver si pueden caber en una sola palabra etiquetada o si tiene que convertirlos en bignums asignados por el montón. Incluso si no se utilizan bignums en la práctica en el tiempo de ejecución, estos controles deberán realizarse. Por otro lado, eso significa que sabes que el algoritmo nunca fallará debido a un envolvente de enteros inesperado si de repente le das entradas más grandes que antes.
Pregunta 4: ¿Mis implementaciones funcionales permiten LCO y, por lo tanto, evitan agregar marcos innecesarios a la pila de llamadas?
Sí, su código Erlang es correcto con respecto a la optimización de la última llamada.
Mirando tu implementación de Erlang. La sincronización ha incluido el inicio de toda la máquina virtual, ejecutar su programa y detener la máquina virtual. Estoy bastante seguro de que configurar y detener el erlang vm lleva algún tiempo.
Si la sincronización se realizó dentro de la propia máquina virtual de erlang, los resultados serían diferentes, ya que en ese caso tendríamos la hora real solo para el programa en cuestión. De lo contrario, creo que el tiempo total tomado por el proceso de inicio y carga de Erlang Vm más el de detenerlo (como lo puso en su programa) se incluye en el tiempo total que el método que está utilizando para calcular El programa está dando salida. Considere usar el tiempo de erlang en sí mismo que usamos cuando queremos cronometrar nuestros programas dentro de la máquina virtual timer:tc/1 or timer:tc/2 or timer:tc/3
. De esta manera, los resultados de erlang excluirán el tiempo necesario para iniciar y detener / detener / detener la máquina virtual. Ese es mi razonamiento allí, piénselo y luego vuelva a probar su punto de referencia.
De hecho, sugiero que intentemos programar el programa (para los idiomas que tienen un tiempo de ejecución), dentro del tiempo de ejecución de esos idiomas para obtener un valor preciso. Por ejemplo, C no tiene gastos generales de iniciar y apagar un sistema de tiempo de ejecución como Erlang, Python y Haskell (98% seguro de esto, estoy en la corrección). Entonces (basado en este razonamiento) concluyo diciendo que este punto de referencia no fue lo suficientemente preciso / justo para los idiomas que se ejecutan en la parte superior de un sistema de tiempo de ejecución. Vamos a hacerlo de nuevo con estos cambios.
EDITAR: además, incluso si todos los idiomas tuvieran sistemas de tiempo de ejecución, la sobrecarga de iniciar cada uno y detenerlo diferiría. por lo que sugiero que hagamos tiempo desde dentro de los sistemas de tiempo de ejecución (para los idiomas para los que se aplica). Se sabe que la máquina virtual de Erlang tiene una sobrecarga considerable al inicio.
Solo por diversión. La siguiente es una implementación de Haskell más ''nativa'':
import Control.Applicative
import Control.Monad
import Data.Either
import Math.NumberTheory.Powers.Squares
isInt :: RealFrac c => c -> Bool
isInt = (==) <$> id <*> fromInteger . round
intSqrt :: (Integral a) => a -> Int
--intSqrt = fromIntegral . floor . sqrt . fromIntegral
intSqrt = fromIntegral . integerSquareRoot''
factorize :: Int -> [Int]
factorize 1 = []
factorize n = first : factorize (quot n first)
where first = (!! 0) $ [a | a <- [2..intSqrt n], rem n a == 0] ++ [n]
factorize2 :: Int -> [(Int,Int)]
factorize2 = foldl (/ls@((val,freq):xs) y -> if val == y then (val,freq+1):xs else (y,1):ls) [(0,0)] . factorize
numDivisors :: Int -> Int
numDivisors = foldl (/acc (_,y) -> acc * (y+1)) 1 <$> factorize2
nextTriangleNumber :: (Int,Int) -> (Int,Int)
nextTriangleNumber (n,acc) = (n+1,acc+n+1)
forward :: Int -> (Int, Int) -> Either (Int, Int) (Int, Int)
forward k val@(n,acc) = if numDivisors acc > k then Left val else Right (nextTriangleNumber val)
problem12 :: Int -> (Int, Int)
problem12 n = (!!0) . lefts . scanl (>>=) (forward n (1,1)) . repeat . forward $ n
main = do
let (n,val) = problem12 1000
print val
Al usar ghc -O3
, esto se ejecuta de manera consistente en 0.55-0.58 segundos en mi máquina (1.73GHz Core i7).
Una función factorCount más eficiente para la versión C:
int factorCount (int n)
{
int count = 1;
int candidate,tmpCount;
while (n % 2 == 0) {
count++;
n /= 2;
}
for (candidate = 3; candidate < n && candidate * candidate < n; candidate += 2)
if (n % candidate == 0) {
tmpCount = 1;
do {
tmpCount++;
n /= candidate;
} while (n % candidate == 0);
count*=tmpCount;
}
if (n > 1)
count *= 2;
return count;
}
Cambiando largos a ints en main, usando gcc -O3 -lm
, esto se ejecuta constantemente en 0.31-0.35 segundos.
Se puede hacer que ambos se ejecuten incluso más rápido si se aprovecha el hecho de que el número n del triángulo = n * (n + 1) / 2, yn y (n + 1) tienen factorizaciones primas completamente diferentes, por lo que el número de factores De cada mitad se puede multiplicar para encontrar el número de factores del conjunto. El seguimiento:
int main ()
{
int triangle = 0,count1,count2 = 1;
do {
count1 = count2;
count2 = ++triangle % 2 == 0 ? factorCount(triangle+1) : factorCount((triangle+1)/2);
} while (count1*count2 < 1001);
printf ("%lld/n", ((long long)triangle)*(triangle+1)/2);
}
reducirá el tiempo de ejecución del código c a 0.17-0.19 segundos, y puede manejar búsquedas mucho más grandes; más de 10000 factores toman aproximadamente 43 segundos en mi máquina. Dejo una aceleración de haskell similar al lector interesado.
Su implementación de Haskell podría acelerarse enormemente mediante el uso de algunas funciones de los paquetes de Haskell. En este caso, utilicé primos, que se instalan simplemente con ''cabal install primes'';)
import Data.Numbers.Primes
import Data.List
triangleNumbers = scanl1 (+) [1..]
nDivisors n = product $ map ((+1) . length) (group (primeFactors n))
answer = head $ filter ((> 500) . nDivisors) triangleNumbers
main :: IO ()
main = putStrLn $ "First triangle number to have over 500 divisors: " ++ (show answer)
Tiempos:
Su programa original:
PS> measure-command { bin/012_slow.exe }
TotalSeconds : 16.3807409
TotalMilliseconds : 16380.7409
Implementación mejorada
PS> measure-command { bin/012.exe }
TotalSeconds : 0.0383436
TotalMilliseconds : 38.3436
Como puede ver, este se ejecuta en 38 milisegundos en la misma máquina donde corrió el suyo en 16 segundos :)
Comandos de compilación:
ghc -O2 012.hs -o bin/012.exe
ghc -O2 012_slow.hs -o bin/012_slow.exe
Usando GHC 7.0.3
, gcc 4.4.6
, Linux 2.6.29
en una máquina x86_64 Core2 Duo (2.5GHz), compilando usando ghc -O2 -fllvm -fforce-recomp
para Haskell y gcc -O3 -lm
para C.
- Su rutina de C se ejecuta en 8.4 segundos (más rápido que su ejecución, probablemente debido a
-O3
) - La solución Haskell se ejecuta en 36 segundos (debido a la
-O2
) - El código de su
factorCount''
no se ha escrito explícitamente y el valorfactorCount''
esInteger
(¡gracias a Daniel por corregir mi diagnóstico erróneo aquí!). Dar una firma de tipo explícita (que es una práctica estándar de todos modos) usandoInt
y el tiempo cambia a 11.1 segundos - en
factorCount''
has llamado innecesariamente desdefromIntegral
. Sin embargo, una solución no produce cambios (el compilador es inteligente, por suerte para usted). - Usaste
mod
donderem
es más rápido y suficiente. Esto cambia el tiempo a 8,5 segundos . -
factorCount''
aplica constantemente dos argumentos adicionales que nunca cambian (number
,sqrt
). Una transformación de trabajador / envoltorio nos da:
$ time ./so
842161320
real 0m7.954s
user 0m7.944s
sys 0m0.004s
Así es, 7,95 segundos . Consistentemente medio segundo más rápido que la solución C. Sin el indicador -fllvm
todavía 8.182 seconds
, por lo que el backend NCG también está funcionando bien en este caso.
Conclusión: Haskell es impresionante.
Código resultante
factorCount number = factorCount'' number isquare 1 0 - (fromEnum $ square == fromIntegral isquare)
where square = sqrt $ fromIntegral number
isquare = floor square
factorCount'' :: Int -> Int -> Int -> Int -> Int
factorCount'' number sqrt candidate0 count0 = go candidate0 count0
where
go candidate count
| candidate > sqrt = count
| number `rem` candidate == 0 = go (candidate + 1) (count + 2)
| otherwise = go (candidate + 1) count
nextTriangle index triangle
| factorCount triangle > 1000 = triangle
| otherwise = nextTriangle (index + 1) (triangle + index + 1)
main = print $ nextTriangle 1 1
EDIT: Así que ahora que hemos explorado eso, vamos a abordar las preguntas
Pregunta 1: ¿Erlang, python y haskell pierden velocidad debido al uso de enteros de longitud arbitraria o no, siempre que los valores sean menores que MAXINT?
En Haskell, usar Integer
es más lento que Int
pero cuánto más lento depende de los cálculos realizados. Por suerte (para máquinas de 64 bits) Int
es suficiente. Por el bien de la portabilidad, probablemente debería volver a escribir mi código para usar Int64
o Word64
(C no es el único idioma con long
).
Pregunta 2: ¿Por qué haskell es tan lento? ¿Hay un indicador de compilador que apague los frenos o es mi implementación? (Lo último es bastante probable, ya que haskell es un libro con siete sellos para mí).
Pregunta 3: ¿Me puede ofrecer algunos consejos sobre cómo optimizar estas implementaciones sin cambiar la forma en que determiné los factores? Optimización de cualquier manera: mejor, más rápido, más "nativo" al idioma.
Eso fue lo que respondí arriba. La respuesta fue
- 0) Utilice la optimización a través de
-O2
- 1) Utilice tipos rápidos (notablemente: unbox-able) cuando sea posible
- 2)
rem
nomod
(una optimización frecuentemente olvidada) y - 3) transformación de trabajador / envoltorio (quizás la optimización más común).
Pregunta 4: ¿Mis implementaciones funcionales permiten LCO y, por lo tanto, evitan agregar marcos innecesarios a la pila de llamadas?
Sí, ese no era el problema. Buen trabajo y me alegra que hayas considerado esto.
C ++ 11, <20ms para mí - Ejecútalo aquí
Entiendo que desea consejos para ayudarlo a mejorar su conocimiento específico del idioma, pero como se ha cubierto bien aquí, pensé que agregaría un contexto para las personas que pueden haber visto el comentario matemático en su pregunta, etc. El código era mucho más lento.
Esta respuesta es principalmente para proporcionar un contexto para ayudar a las personas a evaluar el código en su pregunta / otras respuestas más fácilmente.
Este código utiliza solo un par de optimizaciones (feas), no relacionadas con el lenguaje utilizado, basadas en:
- cada número de traingle es de la forma n (n + 1) / 2
- n y n + 1 son coprime
- El número de divisores es una función multiplicativa.
#include <iostream>
#include <cmath>
#include <tuple>
#include <chrono>
using namespace std;
// Calculates the divisors of an integer by determining its prime factorisation.
int get_divisors(long long n)
{
int divisors_count = 1;
for(long long i = 2;
i <= sqrt(n);
/* empty */)
{
int divisions = 0;
while(n % i == 0)
{
n /= i;
divisions++;
}
divisors_count *= (divisions + 1);
//here, we try to iterate more efficiently by skipping
//obvious non-primes like 4, 6, etc
if(i == 2)
i++;
else
i += 2;
}
if(n != 1) //n is a prime
return divisors_count * 2;
else
return divisors_count;
}
long long euler12()
{
//n and n + 1
long long n, n_p_1;
n = 1; n_p_1 = 2;
// divisors_x will store either the divisors of x or x/2
// (the later iff x is divisible by two)
long long divisors_n = 1;
long long divisors_n_p_1 = 2;
for(;;)
{
/* This loop has been unwound, so two iterations are completed at a time
* n and n + 1 have no prime factors in common and therefore we can
* calculate their divisors separately
*/
long long total_divisors; //the divisors of the triangle number
// n(n+1)/2
//the first (unwound) iteration
divisors_n_p_1 = get_divisors(n_p_1 / 2); //here n+1 is even and we
total_divisors =
divisors_n
* divisors_n_p_1;
if(total_divisors > 1000)
break;
//move n and n+1 forward
n = n_p_1;
n_p_1 = n + 1;
//fix the divisors
divisors_n = divisors_n_p_1;
divisors_n_p_1 = get_divisors(n_p_1); //n_p_1 is now odd!
//now the second (unwound) iteration
total_divisors =
divisors_n
* divisors_n_p_1;
if(total_divisors > 1000)
break;
//move n and n+1 forward
n = n_p_1;
n_p_1 = n + 1;
//fix the divisors
divisors_n = divisors_n_p_1;
divisors_n_p_1 = get_divisors(n_p_1 / 2); //n_p_1 is now even!
}
return (n * n_p_1) / 2;
}
int main()
{
for(int i = 0; i < 1000; i++)
{
using namespace std::chrono;
auto start = high_resolution_clock::now();
auto result = euler12();
auto end = high_resolution_clock::now();
double time_elapsed = duration_cast<milliseconds>(end - start).count();
cout << result << " " << time_elapsed << ''/n'';
}
return 0;
}
Eso toma alrededor de 19 ms en promedio para mi computadora de escritorio y 80 ms para mi computadora portátil, muy lejos de la mayoría de los otros códigos que he visto aquí. Y hay, sin duda, muchas optimizaciones aún disponibles.
Cambio: case (divisor(T,round(math:sqrt(T))) > 500) of
A: case (divisor(T,round(math:sqrt(T))) > 1000) of
Esto producirá la respuesta correcta para el ejemplo de proceso múltiple de Erlang.
Intentando ir
package main
import "fmt"
import "math"
func main() {
var n, m, c int
for i := 1; ; i++ {
n, m, c = i * (i + 1) / 2, int(math.Sqrt(float64(n))), 0
for f := 1; f < m; f++ {
if n % f == 0 { c++ }
}
c *= 2
if m * m == n { c ++ }
if c > 1001 {
fmt.Println(n)
break
}
}
}
Yo obtengo:
versión c original: 9.1690 100%
go: 8.2520 111%
Pero usando:
package main
import (
"math"
"fmt"
)
// Sieve of Eratosthenes
func PrimesBelow(limit int) []int {
switch {
case limit < 2:
return []int{}
case limit == 2:
return []int{2}
}
sievebound := (limit - 1) / 2
sieve := make([]bool, sievebound+1)
crosslimit := int(math.Sqrt(float64(limit))-1) / 2
for i := 1; i <= crosslimit; i++ {
if !sieve[i] {
for j := 2 * i * (i + 1); j <= sievebound; j += 2*i + 1 {
sieve[j] = true
}
}
}
plimit := int(1.3*float64(limit)) / int(math.Log(float64(limit)))
primes := make([]int, plimit)
p := 1
primes[0] = 2
for i := 1; i <= sievebound; i++ {
if !sieve[i] {
primes[p] = 2*i + 1
p++
if p >= plimit {
break
}
}
}
last := len(primes) - 1
for i := last; i > 0; i-- {
if primes[i] != 0 {
break
}
last = i
}
return primes[0:last]
}
func main() {
fmt.Println(p12())
}
// Requires PrimesBelow from utils.go
func p12() int {
n, dn, cnt := 3, 2, 0
primearray := PrimesBelow(1000000)
for cnt <= 1001 {
n++
n1 := n
if n1%2 == 0 {
n1 /= 2
}
dn1 := 1
for i := 0; i < len(primearray); i++ {
if primearray[i]*primearray[i] > n1 {
dn1 *= 2
break
}
exponent := 1
for n1%primearray[i] == 0 {
exponent++
n1 /= primearray[i]
}
if exponent > 1 {
dn1 *= exponent
}
if n1 == 1 {
break
}
}
cnt = dn * dn1
dn = dn1
}
return n * (n - 1) / 2
}
Yo obtengo:
versión c original: 9.1690 100%
versión c de thaumkid: 0.1060 8650%
first go version: 8.2520 111%
second go version: 0.0230 39865%
También probé Python3.6 y pypy3.3-5.5-alpha:
versión c original: 8.629 100%
versión c de thaumkid: 0.109 7916%
Python3.6: 54.795 16%
pypy3.3-5.5-alpha: 13.291 65%
y luego con el siguiente código que tengo:
versión c original: 8.629 100%
versión c de thaumkid: 0.109 8650%
Python3.6: 1.489 580%
pypy3.3-5.5-alpha: 0.582 1483%
def D(N):
if N == 1: return 1
sqrtN = int(N ** 0.5)
nf = 1
for d in range(2, sqrtN + 1):
if N % d == 0:
nf = nf + 1
return 2 * nf - (1 if sqrtN**2 == N else 0)
L = 1000
Dt, n = 0, 0
while Dt <= L:
t = n * (n + 1) // 2
Dt = D(n/2)*D(n+1) if n%2 == 0 else D(n)*D((n+1)/2)
n = n + 1
print (t)
Algunos números más y explicaciones para la versión C. Al parecer nadie lo hizo durante todos esos años. Recuerde aumentar esta respuesta para que todos puedan ver y aprender.
Paso uno: Benchmark de los programas del autor.
Especificaciones de la computadora portátil:
- CPU i3 M380 (931 MHz - modo de ahorro máximo de batería)
- Memoria de 4GB
- Win7 64 bits
- Microsoft Visual Studio 2012 Ultimate
- Cygwin con gcc 4.9.3
- Python 2.7.10
Comandos:
compiling on VS x64 command prompt > `for /f %f in (''dir /b *.c'') do cl /O2 /Ot /Ox %f -o %f_x64_vs2012.exe`
compiling on cygwin with gcc x64 > `for f in ./*.c; do gcc -m64 -O3 $f -o ${f}_x64_gcc.exe ; done`
time (unix tools) using cygwin > `for f in ./*.exe; do echo "----------"; echo $f ; time $f ; done`
.
----------
$ time python ./original.py
real 2m17.748s
user 2m15.783s
sys 0m0.093s
----------
$ time ./original_x86_vs2012.exe
real 0m8.377s
user 0m0.015s
sys 0m0.000s
----------
$ time ./original_x64_vs2012.exe
real 0m8.408s
user 0m0.000s
sys 0m0.015s
----------
$ time ./original_x64_gcc.exe
real 0m20.951s
user 0m20.732s
sys 0m0.030s
Los nombres de archivo son: integertype_architecture_compiler.exe
- integertype es el mismo que el programa original por ahora (más sobre esto más adelante)
- La arquitectura es x86 o x64 dependiendo de la configuración del compilador.
- compilador es gcc o vs2012
Paso dos: Investigar, mejorar y comparar nuevamente
VS es 250% más rápido que gcc. Los dos compiladores deberían dar una velocidad similar. Obviamente, algo está mal con el código o las opciones del compilador. ¡A investigar!
El primer punto de interés son los tipos enteros. Las conversiones pueden ser costosas y la consistencia es importante para una mejor generación de código y optimizaciones. Todos los enteros deben ser del mismo tipo.
Es un lío mixto de int
y en long
este momento. Vamos a mejorar eso. ¿Qué tipo de uso? El más rápido. Tienes que compararlos a todos!
----------
$ time ./int_x86_vs2012.exe
real 0m8.440s
user 0m0.016s
sys 0m0.015s
----------
$ time ./int_x64_vs2012.exe
real 0m8.408s
user 0m0.016s
sys 0m0.015s
----------
$ time ./int32_x86_vs2012.exe
real 0m8.408s
user 0m0.000s
sys 0m0.015s
----------
$ time ./int32_x64_vs2012.exe
real 0m8.362s
user 0m0.000s
sys 0m0.015s
----------
$ time ./int64_x86_vs2012.exe
real 0m18.112s
user 0m0.000s
sys 0m0.015s
----------
$ time ./int64_x64_vs2012.exe
real 0m18.611s
user 0m0.000s
sys 0m0.015s
----------
$ time ./long_x86_vs2012.exe
real 0m8.393s
user 0m0.015s
sys 0m0.000s
----------
$ time ./long_x64_vs2012.exe
real 0m8.440s
user 0m0.000s
sys 0m0.015s
----------
$ time ./uint32_x86_vs2012.exe
real 0m8.362s
user 0m0.000s
sys 0m0.015s
----------
$ time ./uint32_x64_vs2012.exe
real 0m8.393s
user 0m0.015s
sys 0m0.015s
----------
$ time ./uint64_x86_vs2012.exe
real 0m15.428s
user 0m0.000s
sys 0m0.015s
----------
$ time ./uint64_x64_vs2012.exe
real 0m15.725s
user 0m0.015s
sys 0m0.015s
----------
$ time ./int_x64_gcc.exe
real 0m8.531s
user 0m8.329s
sys 0m0.015s
----------
$ time ./int32_x64_gcc.exe
real 0m8.471s
user 0m8.345s
sys 0m0.000s
----------
$ time ./int64_x64_gcc.exe
real 0m20.264s
user 0m20.186s
sys 0m0.015s
----------
$ time ./long_x64_gcc.exe
real 0m20.935s
user 0m20.809s
sys 0m0.015s
----------
$ time ./uint32_x64_gcc.exe
real 0m8.393s
user 0m8.346s
sys 0m0.015s
----------
$ time ./uint64_x64_gcc.exe
real 0m16.973s
user 0m16.879s
sys 0m0.030s
Los tipos enteros son int
long
int32_t
uint32_t
int64_t
y uint64_t
desde#include <stdint.h>
Hay LOTES de tipos de enteros en C, más algunos con signo / sin signo para jugar, más la opción de compilar como x86 o x64 (no debe confundirse con el tamaño del entero real). Eso es un montón de versiones para compilar y ejecutar ^^
Paso Tres: Entendiendo los Números
Conclusiones definitivas:
- Los enteros de 32 bits son ~ 200% más rápidos que los equivalentes de 64 bits
- los enteros sin signo de 64 bits son 25% más rápidos que los de 64 bits con signo (Desafortunadamente, no tengo una explicación para eso)
Pregunta del truco: "¿Cuáles son los tamaños de int y long en C?"
La respuesta correcta es: ¡ El tamaño de int y long en C no está bien definido!
De la especificación de C:
Int tiene al menos 32 bits de
longitud es al menos un Int
Desde la página del manual de gcc (banderas -m32 y -m64):
El entorno de 32 bits establece int, long y puntero a 32 bits y genera código que se ejecuta en cualquier sistema i386.
El entorno de 64 bits establece int en 32 bits y largo y puntero en 64 bits y genera código para la arquitectura x86-64 de AMD.
De la documentación de MSDN (Rangos de tipos de datos) https://msdn.microsoft.com/en-us/library/s3f49ktz%28v=vs.110%29.aspx :
int, 4 bytes, también conocido como firmado
largo, 4 bytes, también conocido como largo int y firmado largo int
Para concluir esto: Lecciones aprendidas
Los enteros de 32 bits son más rápidos que los enteros de 64 bits.
Los tipos de enteros estándar no están bien definidos en C ni en C ++, varían según los compiladores y las arquitecturas. Cuando necesite coherencia y previsibilidad, use la
uint32_t
familia de enteros de#include <stdint.h>
.Problemas de velocidad resueltos. ¡Todos los demás idiomas están atrasados cientos de por ciento, C y C ++ ganan nuevamente! Siempre lo hacen La próxima mejora será multiproceso usando OpenMP: D
Modifiqué la versión de "Jannich Brendle" a 1000 en lugar de 500. Y listé el resultado de euler12.bin, euler12.erl, p12dist.erl. Ambos códigos erl usan ''+ nativo'' para compilar.
zhengs-MacBook-Pro:workspace zhengzhibin$ time erl -noshell -s p12dist start
The result is: 842161320.
real 0m3.879s
user 0m14.553s
sys 0m0.314s
zhengs-MacBook-Pro:workspace zhengzhibin$ time erl -noshell -s euler12 solve
842161320
real 0m10.125s
user 0m10.078s
sys 0m0.046s
zhengs-MacBook-Pro:workspace zhengzhibin$ time ./euler12.bin
842161320
real 0m5.370s
user 0m5.328s
sys 0m0.004s
zhengs-MacBook-Pro:workspace zhengzhibin$
Supuse que el número de factores solo es grande si los números involucrados tienen muchos factores pequeños. Así que utilicé el excelente algoritmo de thaumkid, pero primero usé una aproximación al factor de conteo que nunca es demasiado pequeño. Es bastante simple: compruebe los factores primos hasta 29, luego verifique el número restante y calcule un límite superior para el número de factores. Use esto para calcular un límite superior para el número de factores, y si ese número es lo suficientemente alto, calcule el número exacto de factores.
El código a continuación no necesita esta suposición para la corrección, sino para ser rápido. Parece funcionar; solo uno de cada 100,000 números da una estimación que es lo suficientemente alta como para requerir una verificación completa.
Aquí está el código:
// Return at least the number of factors of n.
static uint64_t approxfactorcount (uint64_t n)
{
uint64_t count = 1, add;
#define CHECK(d) /
do { /
if (n % d == 0) { /
add = count; /
do { n /= d; count += add; } /
while (n % d == 0); /
} /
} while (0)
CHECK ( 2); CHECK ( 3); CHECK ( 5); CHECK ( 7); CHECK (11); CHECK (13);
CHECK (17); CHECK (19); CHECK (23); CHECK (29);
if (n == 1) return count;
if (n < 1ull * 31 * 31) return count * 2;
if (n < 1ull * 31 * 31 * 37) return count * 4;
if (n < 1ull * 31 * 31 * 37 * 37) return count * 8;
if (n < 1ull * 31 * 31 * 37 * 37 * 41) return count * 16;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43) return count * 32;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43 * 47) return count * 64;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43 * 47 * 53) return count * 128;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43 * 47 * 53 * 59) return count * 256;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43 * 47 * 53 * 59 * 61) return count * 512;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43 * 47 * 53 * 59 * 61 * 67) return count * 1024;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43 * 47 * 53 * 59 * 61 * 67 * 71) return count * 2048;
if (n < 1ull * 31 * 31 * 37 * 37 * 41 * 43 * 47 * 53 * 59 * 61 * 67 * 71 * 73) return count * 4096;
return count * 1000000;
}
// Return the number of factors of n.
static uint64_t factorcount (uint64_t n)
{
uint64_t count = 1, add;
CHECK (2); CHECK (3);
uint64_t d = 5, inc = 2;
for (; d*d <= n; d += inc, inc = (6 - inc))
CHECK (d);
if (n > 1) count *= 2; // n must be a prime number
return count;
}
// Prints triangular numbers with record numbers of factors.
static void printrecordnumbers (uint64_t limit)
{
uint64_t record = 30000;
uint64_t count1, factor1;
uint64_t count2 = 1, factor2 = 1;
for (uint64_t n = 1; n <= limit; ++n)
{
factor1 = factor2;
count1 = count2;
factor2 = n + 1; if (factor2 % 2 == 0) factor2 /= 2;
count2 = approxfactorcount (factor2);
if (count1 * count2 > record)
{
uint64_t factors = factorcount (factor1) * factorcount (factor2);
if (factors > record)
{
printf ("%lluth triangular number = %llu has %llu factors/n", n, factor1 * factor2, factors);
record = factors;
}
}
}
}
Esto encuentra el 14.753.024 triángulos triangulares con 13824 factores en aproximadamente 0,7 segundos, el número triangular 879.207.615 con 61.440 factores en 34 segundos, el triángulo número 12.524.486.975 con 138.240 factores en 10 minutos 5 segundos y el triángulo 26.467.792.064 en 172 21 minutos y 25 segundos (2.4GHz Core2 Duo), por lo que este código toma solo 116 ciclos de procesador por número en promedio. El último número triangular en sí es mayor que 2 ^ 68, por lo que
#include <stdio.h>
#include <math.h>
int factorCount (long n)
{
double square = sqrt (n);
int isquare = (int) square+1;
long candidate = 2;
int count = 1;
while(candidate <= isquare && candidate<=n){
int c = 1;
while (n % candidate == 0) {
c++;
n /= candidate;
}
count *= c;
candidate++;
}
return count;
}
int main ()
{
long triangle = 1;
int index = 1;
while (factorCount (triangle) < 1001)
{
index ++;
triangle += index;
}
printf ("%ld/n", triangle);
}
gcc -lm -Ofast euler.c
tiempo ./a.out
2.79s usuario 0.00s sistema 99% cpu 2.794 total