How much can you quickly multiply?

With the recent python bashing , here's an attempt to show python's strengths. Your challenge is to write a program that calculates the factorial of as high a number `n` as possible within 10 seconds.

Your score will be `(highest n for your program on your machine)/(highest n for my program on your machine)`

Rules

• You must calculate an exact integer solution. Since the factorial would be much higher than what can fit in a 64 bit unsigned integer, you can use strings if your language does not support large integers
• Standard loopholes are forbidden. Particularly, you cannot use any external resources.
• Only the calculation part(this includes time for any workarounds using strings) adds to the total time which should be under 10 seconds on average.
• You must store the output in an easily printable form (as printing takes time) (see my program below), string, variable, character array, etc.

EDIT:

• Your program must give the correct output for all `n`: `1 <= n <= (your highest n)`

EDIT2:

My program

``````from __future__ import print_function
import time

def factorial( n ):
return reduce( ( lambda x , y : x * y ) , xrange( 1 , n + 1 ) , 1 )

start = time.clock()
end = time.clock()

print ( "Time:" , end - start , "sec" )
``````

Highest score wins. For the record, my code can manage `n = 90000` in about `9.89` seconds on a Pentium 4 3.0 GHz

EDIT: Can everyone please add the score rather than just the highest n. Just the highest `n` has no meaning by itself as it depends on your hardware. It's impossible to have an objective winning criterion otherwise. ali0sha's anwer does this correctly.

We have a winner. I didn't accept the java answer http://codegolf.stackexchange.com/a/26974/8766 as it kind of skirts close to http://meta.codegolf.stackexchange.com/a/1080/8766

-
You can use `operator.mul` instead of the lambda function – gnibbler May 6 '14 at 11:57
Bit suprised this works, but assuming I read the rules correctly this MATLAB solution would be nice: `factorial(Inf)`, returns `Inf` in a fraction of a second. – Dennis Jaheruddin May 6 '14 at 12:40
@Doorknob That fits in the standard loopholes. – Justin May 6 '14 at 12:44
@DennisJaheruddin, it's a bit of a stretch to refer to "Inf" as an "exact integer solution". – tobyink May 6 '14 at 13:12
@Quincunx No, any language is allowed. – user80551 May 6 '14 at 17:22

C++ with GMP, score = 55.55 (10,000,000 / 180,000)

``````#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <vector>
#include <iostream>
#include <queue>
#include <gmpxx.h>

int main(int argc, char *argv[]) {
uint64_t n = atoi(argv[1]);

// Iterate through 1..n.  Strip off powers of 2.  Multiply
// remainders together into <= 64 bit chunks.
uint64_t twos = 0;
std::vector<uint64_t> terms;
uint64_t m = 1;
for(uint64_t i = 1; i <= n; i++) {
uint64_t j = __builtin_ctzll(i);
twos += j;
uint64_t k = i >> j;
if(__builtin_clzll(m) + __builtin_clzll(k) >= 64) {
m *= k;
} else {
terms.push_back(m);
m = k;
}
}
if(m != 1) terms.push_back(m);

// convert to gmp
// why isn't there a 64-bit constructor?
std::queue<mpz_class> gmpterms;
for(int i = 0; i < terms.size(); i++) {
mpz_class x = (uint32_t)(terms[i] >> 32);
x <<= 32;
x += (uint32_t)terms[i];
gmpterms.push(x);
}

// pop two from the bottom, multiply them, push on the end.
while(gmpterms.size() > 1) {
mpz_class a = gmpterms.front();
gmpterms.pop();
mpz_class b = gmpterms.front();
gmpterms.pop();
gmpterms.push(a * b);
}

mpz_class r = gmpterms.front();
r <<= twos;
//std::cout << r << std::endl;
}
``````
-

42.575 = ( 6,812,000 / 160,000 ) approx

Code:

``````import gmpy2

def fac1(n):
m=lambda(L):([]if len(L)%2==0 else[L.pop()])+map(lambda(l):l[0]*l[1],zip(L[1::2],L[-2::-2]))
L=map(gmpy2.mpz,xrange(1,n+1))
Number = (len(L)-1).bit_length()
while Number:Number-=1;L=m(L)
return L[0]

def fac2(n):
global E; E=0
def f(i):
global E; E+=i//2
return[]if i==1 else f(i//2)+range(3,i,2)+[[1,i][i%2]]
m=lambda(L):([]if len(L)%2==0 else[L.pop()])+map(lambda(l):l[0]*l[1],zip(L[1::2],L[-2::-2]))
L=map(gmpy2.mpz,f(n))
N=(len(L)-1).bit_length()
while N: N-=1;L=m(L)
return L[0]<<E
``````

Test:

``````import time

start = time.time()
baseline(160000)
print time.time()-start

start = time.time()
fac1(6811000)
print time.time()-start

start = time.time()
fac2(6812000)
print time.time()-start

start = time.time()
gmpy2.fac(26000000)
print time.time()-start
``````

Output:

``````10.0069999695
10.0729999542
10.0360000134
9.98699998856
``````

How it works:

Bigger multiplications take more time, thus we want to do as many small multiplications as possible. This is especially true in Python where for numbers less that `2^64` we use hardware arithmetic, and above that we use software. So, in `m(L)`, we start with a list `L`; if it's odd length we remove one number from consideration to make it even again. Then we multiply element `1` with element `-2`, element `3` with `-4`, etc, so that

``````m([1,2,3,4,5,6,7,8]) = [2*7, 4*5, 6*3, 8*1] = [14, 20, 18, 8]
m([10,12,6]) = [360,112]
m([120,6]) = [40320]
``````

This approach ensures we're using hardware arithmetic for as long as possible, following which we switch onto the efficient gmc arithmetic library.

In `fac2`, we take a more classic divide and conquer approach as well, where we split out every multiple of 2 and bitshift them at the end for a trivial performance boost. I've included it here because it's usually around 0.5% faster than `fac1`.

Golfed Version of `fac1` (because I can), 220B

``````import gmpy2
def f(n):
m=lambda(L):([]if len(L)%2==0 else[L.pop()])+map(lambda(l):l[0]*l[1],zip(L[1::2],L[-2::-2]))
L=map(gmpy2.mpz,xrange(1,n+1));N=(len(L)-1).bit_length()
while N:N-=1;L=m(L)
return L[0]
``````
-
If the GMP backend includes a bitshift function then you can keep the numbers even smaller by dividing each number in the list by 2 until it's even and then doing a single shift at the end. – Peter Taylor May 6 '14 at 15:58
Where did you get `gmpy2` from? \$ python Python 2.7.3 (default, Feb 27 2014, 19:58:35) [GCC 4.6.3] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> from gmpy2 import mpz Traceback (most recent call last): File "<stdin>", line 1, in <module> ImportError: No module named gmpy2 >>> – user80551 May 6 '14 at 17:57
@user80551: code.google.com/p/gmpy (the top google search result) has installers for many different platforms. – alexander-brett May 6 '14 at 19:35
For the golfed version, couldn't you do `while len(L): ...` instead of `while len(L)>1: ...`? – user80551 May 7 '14 at 9:10
No: the function inside that loop will never take the list below length 1 and anyhow we need the first element! – alexander-brett May 7 '14 at 9:12

Python 3, n = 100000

A simple algorithm change was all that was needed to bump the sample code up by 10000.

``````import time

def factorial(n):
result = 1
while n > 0:
result *= n
n = n - 1
return result

start = time.clock()
end = time.clock()

print("Time:", end - start, "sec")
``````

Obviously not the most creative answer, but there's really only one way to do a factorial....

-
Please give the score, see my edit. The bump would probably be because your machine is better than mine. – user80551 May 6 '14 at 17:46

Perl + C, n = about 3 million

Here I'm using the Math::BigInt::GMP library available on CPAN, which provides a massive speed boost for Perl's core Math::BigInt objects.

``````use v5.14;
use Time::HiRes 'time';
use Math::BigInt only => 'GMP';

sub factorial { Math::BigInt::->new(@_)->bfac }

my \$start  = time;
my \$answer = factorial( 3_000_000 );
my \$end    = time;

say "Time: ", \$end - \$start, " sec";
``````

Bear in mind that my computer is probably quite a bit slower than yours. Using your original Python script, I can only calculate `factorial(40000)` in 10 seconds; `factorial(90000)` takes a lot longer. (I hit Ctrl+C after a minute.) On your hardware, using Math::BigInt::GMP, you may well be able to calculate the factorial of 5 million or more in under 10 seconds.

One thing you may notice is that although the factorial is calculated incredibly quickly, printing out the result is very slow, taking about three times longer than the original calculation. This is because GMP internally uses a binary rather than decimal representation, and printing it out requires binary to decimal conversion.

-
I think GMP counts as an external resource. (Although it certainly makes things a hell of a lot easier than implementing prime factorization and Schönhage-Strassen multiplication from scratch.) – squeamish ossifrage May 6 '14 at 13:14
I was assuming that "external resource" referred to looking up solutions from a pre-computed set of answers in a database, or a web service, etc. – tobyink May 6 '14 at 13:29
Squeamish: libraries don't normally count as external resources unless they have a function which falls under the boring loopholes rule. – alexander-brett May 6 '14 at 17:01
Tobyink: can you explain what your program does? it looks like you're just using a built-in function (bfac?) – alexander-brett May 6 '14 at 17:27
Yup. This answer is invalid, since it uses the factorial method of `Math::BigInt` – 14mRh4X0r May 9 '14 at 10:34

Python 2.7
5.94 = 1'200'000/202'000

``````def fast_fac(n):
def prod(start, fin):
if fin - start <= 50:
return reduce(lambda x,y: x*y, xrange(start, fin+1), 1)
else:
mid = (start+fin) / 2
return prod(start, mid) * prod(mid+1, fin)
return prod(1, n)
``````

Makes use of relative ease of multiplication of many groups of small numbers and then multiplying them compared to large number of multiplyings involving huge number.

-

C#: 0,48 (77,000 / 160,000)

I'm not happy with this.

Is C# that slow?

But here is my entry anyway.

``````static void Main(string[] args)
{
Console.WriteLine("Enter N for fatorial:");

Stopwatch s = Stopwatch.StartNew();

BigInteger result = 1;
while (0 <-- n) result *= n;

s.Stop();

Console.WriteLine("Output: {0} ", result);

Console.WriteLine("Completed in {0}", s.Elapsed);

}
``````

When n = 77000 it takes `00:00:09:8708952` to calculate.

I'm running in Release mode, outside Visual Studio, using a Core i3-2330M @2.2GHz.

Edit: Since i'm not doing anything intelligent, I accept that result. Maybe the .NET Framework 4.5 is addind some overhead (or BigInteger isn't that fast).

-
Please give the score and not just `n` – user80551 May 6 '14 at 18:10
You could use `zero approached by` operator to make it prettier (like start with `n = ... + 1` then do `while (0 <-- n) result *= n;`) – Cthulhu May 6 '14 at 18:50
Thanks, @Cthulhu, it looks awesome. – Ricardo Pieper May 6 '14 at 19:26
BigInteger for .NET probably did not implement algorithms for multiplying larger numbers, like Karatsuba or Toom-3. If so, this is a good example of how Python is faster. – kernigh May 11 '14 at 20:20

bc, score=0.19

What the heck, here's my contender for "How much can you slowly multiply?"

bc is "An arbitrary precision calculator language", but unfortunately rather slow:

``````n=read()
for(f=i=1;i<=n;i++)f*=i
f
quit
``````

In about 10 seconds on my mid 2012 MacBook Pro (2.3 GHz Intel Core i7) the reference python answer can calculate 122000!, but this bc script can only calculate 23600!.

Conversely 10000! takes 1.5s with the python reference script, but the bc script takes 50s.

Oh dear.

-
OpenBSD bc(1) is faster. Your program scores 0.29 = 28000/98000. There is no `read()`, so I ran `time sed 's/read()/28000/' factorial.bc | bc`. – kernigh May 11 '14 at 15:33

C++ (x86_64-specific) - 3.0 (390000/130000)

(easily portable to x86-32, porting to other architectures implies a significant speed loss)

Here's my own micro-implementation of long arithmetic.
The calculation itself takes 10 seconds, and while the output is in easily printable form (see the `operator<<` overload), it takes some more time to print it.

``````#include <vector>
#include <iostream>
#include <stdint.h>
#include <ctime>

typedef uint64_t digit;
typedef std::vector<digit> number;

std::ostream &operator<<(std::ostream &s, const number &x)
{
std::vector<char> o;
size_t size = x.size() * 21;
o.resize(size);
size_t lud = 0;
for(number::const_reverse_iterator i = x.rbegin(), end = x.rend(); i != end; i++)
{
digit carry = 0;
int j;
for(j = 0; j <= lud || carry; j++)
{
digit r = o[j] * (1LL << 32) + carry;
o[j] = r % 10;
carry = r / 10;
}
lud = j;
carry = 0;
for(j = 0; j <= lud || carry; j++)
{
digit r = o[j] * (1LL << 32) + carry;
o[j] = r % 10;
carry = r / 10;
}
lud = j;
carry = *i;
for(j = 0; carry; j++)
{
digit r = o[j] + (carry % 10);
carry /= 10;
carry += r / 10;
o[j] = r % 10;
}
if(j > lud)
lud = j;
}
for(int j = lud; j--;)
s.put(o[j] + '0');
return s;
}

inline uint64_t dmul(uint64_t x, uint64_t y, uint64_t &carry)
{
asm("mulq %2" : "+a"(x), "=d"(carry) : "r"(y));
return x;
}
inline digit dadd(digit x, digit y, digit &carry)
{
asm("movq \$0, %1; addq %2, %0; adcq %1, %1" : "+r"(x), "=r"(carry), "+r"(y));
return x;
}

void multiply(number &x, digit y)
{
x.resize(x.size() + 2);
digit carry = 0;
for(number::iterator i = x.begin(), end = x.end(); i != end; i++)
{
digit nc, res = dmul(*i, y, nc);
carry += nc;
}
size_t sz = x.size();
for(number::const_reverse_iterator i = x.rbegin(), end = x.rend(); i != end; i++)
{
if(*i)
break;
sz--;
}
x.resize(sz);
}

int main()
{
const int r = 390000;
clock_t start = clock();
number n;
digit mult = 1;
n.push_back(1);
for(digit a = 2; a <= r; a++)
{
digit carry, m = dmul(mult, a, carry);
if(carry)
{
multiply(n, mult);
mult = a;
}
else
mult = m;
}
multiply(n, mult);
std::cout << "Took: " << (clock() - start)/((double)CLOCKS_PER_SEC) << std::endl;
std::cout << n << std::endl;
}
``````
-
Check your score. You need to run the question's Python 2.7 program on your computer. For my computer, I compiled your program with `g++ -O2 factorial.cc -o factorial` and it scores 3.90 = 382000 / 98000. – kernigh May 11 '14 at 20:11
Weird, I got 3.9 and you got 3.0 for this program. I guess your faster computer is a penalty. Perhaps your program loses its advantage over Python as `r` increases. If so, and you can do a higher `r` in 10 seconds, then your score goes down. – kernigh May 12 '14 at 21:30

Python 3: 280000 / 168000

Time running your program: between `9.87585953253` and `10.3046453994`. Time running my program: about `10.35296977897559`.

``````import time

def factorial(n):
f = 1
while n > 1:
hn = n >> 1
f = f * 2**hn * double_factorial(n) #dfl[hn + (n & 1) - 1]
n = hn
return f
def double_factorial(n):
#dfl = [1]
p = 1
l = 3
mh = n
while l <= n:
p *= l
l += 2
#dfl.append(p)
return p

start = time.clock()
factorial(280000)
end = time.clock()

print(end - start)
``````

I read this answer on cs.SE and decided to try to implement it in Python. However, I accidentally discovered that `n! = (⌊n / 2⌋)! * 2**(⌊n / 2⌋) * n!!` (note: `!!` is the double factorial). So I converted that to a non-recursive form.

The comments show my attempt to avoid recomputing the double factorial, but I discovered that storing every value was too memory-costly that it caused my computer to run even slower. I can improve this by only storing what is needed.

Strangely, I implemented the naive straight multiplication in Python 3, and it does better than your program: `n = 169000` in 10 seconds.:

``````def factorial(n):
p=1
for i in range(n):
p*=i+1
return p
``````
-

Bash: score = 0.001206 (181/150000)

I stole the math functions from Rosettacode - Long multiplication I didn't analyzed nor tried to optimize.
You are free to change the algorithm or to try a different strings split method.

``````#!/bin/bash

if (( \${#1} < \${#2} )); then
local a="\$2" b="\$1" sum= carry=0
else
local a="\$1" b="\$2" sum= carry=0
fi

while (( \${#a} )); do
local -i d1="\${a##\${a%?}}" d2="10#0\${b##\${b%?}}" s=carry+d1+d2
sum="\${s##\${s%?}}\$sum"
carry="10#0\${s%?}"
a="\${a%?}" b="\${b%?}"
done
echo "\$sum"
}

multiply() { # arbitrary-precision multiplication
if (( \${#1} < \${#2} )); then
local a="\$2" b="\$1" product=0
else
local a="\$1" b="\$2" product=0
fi

local zeroes=
while (( \${#b} )); do
local m1="\$a"
local m2="\${b##\${b%?}}"
local partial=\$zeroes
local -i carry=0
while (( \${#m1} )); do
local -i d="\${m1##\${m1%?}}"
m1="\${m1%?}"
local -i p=d*m2+carry
partial="\${p##\${p%?}}\$partial"
carry="10#0\${p%?}"
done
partial="\${carry#0}\$partial"
zeroes=0\$zeroes
b="\${b%?}"
done
echo "\$product"
}

# 'timerun' function
trap 'echo \$((i -1)) \$f; exit'  USR1
(sleep 9.9; kill -USR1 \$\$)&

declare -i i
f=1
for ((i=1; i< 10000 ; i++ ))   # 10000 is verry optimistic
do
f=\$(multiply \$f \$i)
done
``````
-
Please add the score and not just the highest n – user80551 May 7 '14 at 18:21
@user80551 it's done – Emmanuel May 7 '14 at 22:05

Python 3, advanced algo by Peter Luschny: 8.25x (1 280 000/155 000)

Shamelessly copied from Peter Luschny,
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm,

This is actually a quite advanced algorithm, using something called the "swinging factorial" and a list of primes. I suspect it could be even faster if it did like many of the other answers and performed most of the multiplications with 32 bit integers.

``````#! /usr/bin/python3
import time
import bisect

def Primes(n) :
primes = [2, 3]
lim, tog = n // 3, False
composite = [False for i in range(lim)]

d1 = 8; d2 = 8; p1 = 3; p2 = 7; s = 7; s2 = 3; m = -1

while s < lim :             # --  scan the sieve
m += 1                  # --  if a prime is found
if not composite[m] :   # --  cancel its multiples
inc = p1 + p2
for k in range(s,      lim, inc) : composite[k] = True
for k in range(s + s2, lim, inc) : composite[k] = True

tog = not tog
if tog: s += d2; d1 += 16; p1 += 2; p2 += 2; s2 = p2
else:   s += d1; d2 +=  8; p1 += 2; p2 += 6; s2 = p1

k, p, tog = 0, 5, False
while p <= n :
if not composite[k] : primes.append(p)
k += 1;
tog = not tog
p += 2 if tog else 4

return primes

def isqrt(x):
'''
Writing your own square root function
'''
if x < 0: raise ValueError('square root not defined for negative numbers')
n = int(x)
if n == 0: return 0
a, b = divmod(n.bit_length(), 2)
x = 2**(a + b)
while True:
y = (x + n // x) // 2
if y >= x: return x
x = y

def product(s, n, m):
if n > m: return 1
if n == m: return s[n]
k = (n + m) // 2
return product(s, n, k) * product(s, k + 1, m)

def factorialPS(n):

small_swing = [1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435,
109395,12155,230945,46189,969969,88179,2028117,676039,16900975,
1300075,35102025,5014575,145422675,9694845,300540195,300540195]

def swing(m, primes):
if m < 33: return small_swing[m]

s = bisect.bisect_left(primes, 1 + isqrt(m))
d = bisect.bisect_left(primes, 1 + m // 3)
e = bisect.bisect_left(primes, 1 + m // 2)
g = bisect.bisect_left(primes, 1 + m)

factors = primes[e:g]
factors += filter(lambda x: (m // x) & 1 == 1, primes[s:d])
for prime in primes[1:s]:
p, q = 1, m
while True:
q //= prime
if q == 0: break
if q & 1 == 1:
p *= prime
if p > 1: factors.append(p)

return product(factors, 0, len(factors) - 1)

def odd_factorial(n, primes):
if n < 2: return 1
return (odd_factorial(n // 2, primes)**2) * swing(n, primes)

def eval(n):
if n < 0:
raise ValueError('factorial not defined for negative numbers')

if n == 0: return 1
if n < 20: return product(range(2, n + 1), 0, n-2)

N, bits = n, n
while N != 0:
bits -= N & 1
N >>= 1

primes = Primes(n)
return odd_factorial(n, primes) * 2**bits

return eval(n)

start = time.time()
print(time.time()-start)
``````
-
+1 for extreme algo effort :) – alexander-brett May 12 '14 at 8:26

Java - 10.9

n = 885000

Mergesort-y.

``````import java.math.BigInteger;

public class Factorials {

public static BigInteger fac;

public static BigInteger two = BigInteger.valueOf(2);

static BigInteger mul(BigInteger start, BigInteger end) {
if(start.equals(end)) {
return start;
} else {
}
}

public static void main(String[] args) {
Factorials.fac = BigInteger.valueOf(Integer.parseInt(args[0]));
long t = System.nanoTime();
BigInteger result = mul(BigInteger.ONE, fac);
t = System.nanoTime() - t;
System.out.print(String.valueOf(((float) t) / 1000000000)); //result.toString()+" @ "+
}
}
``````

`BigInteger`s are slow.

Recommendations for arbitrary-precision high-speed Java integer libraries? :P

-
Can I steal your code to make it multithreaded? – Simon Kuang May 11 '14 at 6:47
@SimonKuang Go ahead. :P Multi-threaded entries aren't allowed here, though. Also, you might want to use a more efficient BigInteger implementation. – cjfaure May 11 '14 at 10:13

score = 1.80 = 176_000 / 98_000

EDIT: improved from 1.35 = 132_000 / 98_000

I took ideas from the GMP factorial algorithm. This program uses the standard library to generate prime numbers. Ruby is a bad choice because multiplication seems slower in Ruby than in Python.

1. My program in Ruby 2.1: score = 1.80 = 176_000 / 98_000
2. Trivial algorithm in Python 2.7: score = 1 = 98_000 / 98_000
3. Trivial algorithm in Ruby 2.1: score = 0.878 = 86_000 / 98_000

Yes, my binary of `ruby 2.1.0p0 (2013-12-25 revision 44422) [x86_64-openbsd]` links against GMP. Ruby 2.1 added a feature to use GMP for large multiplication, but it still seems slower than Python 2.7.

``````require 'benchmark'
require 'optparse'
require 'prime'

def factorial(n)
# calculate primes up to n, drop the 2
@odd_primes = Prime.each(n).drop(1)

# count prime factors of factorial(n)
@factors = Hash.new(0)
factorial_recurse(n)

shift = @factors.delete(2) || 0
@factors.inject(1) {|product, (base, exp)|
product * base**exp
} << shift
end

def factorial_recurse(n)
return if n < 2

# collect prime factors of 2 * 4 * 6 * .. * n
#  = (2 * 2 * 2 * .. * 2) * (1 * 2 * 3 * .. * exp)
#  = 2**exp * factorial(exp) where exp = floor(n/2)
exp = n >> 1
factorial_recurse(exp)
@factors[2] += exp

# collect prime factors 3 * 5 * 7 * ... * n
for prime in @odd_primes
break if prime > n
exp = 0
# count occurences of prime, prime**2, prime**3, .. n
prime_power = prime
until prime_power > n
# floor(n / prime_power) occurences in 1 * 2 * .. * n,
# but only ceil(count / 2) occurences in 3 * 5 * .. * n
@factors[prime] += (n / prime_power + 1) >> 1
prime_power *= prime
end
end
end

# usage: factorial.rb [-ct] [number]
cflag = tflag = false
OptionParser.new {|opts|
opts.on('-c', 'Check for bugs') { cflag = true }
opts.on('-t', 'Use trivial algorithm') { tflag = true }
opts.parse!
}
\$*[1] and fail 'too many arguments'
n = Integer(\$*[0] || 176_000)

if cflag
factorial(n) == (1..n).reduce(1, :*) or
fail "bad program: factorial(#{n}) is wrong"
puts "ok"
exit
end

# measure processor time to calculate factorial
f = nil
if tflag
time = Benchmark.measure { f = (1..n).reduce(1, :*) }
else
time = Benchmark.measure { f = factorial(n) }
end
puts f
puts "Time #{time.total} sec"
``````
-

Java - 125.15 (21,400,000 / 171,000)

Also shamelessly copied from Peter Luschny's Github repo (thanks @semi-extrinsic) and licensed under the MIT license, this uses the "prime factorization nested squaring" algorithm as proposed by Albert Schönhage et al. (according to Luschny's factorial algorithms description page).

I slightly adapted the algorithm to use Java's BigInteger and to not use a lookup table for n < 20.

Compiled with gcj, which uses GMP for its BigInteger implementation, and ran on Linux 3.12.4 (Gentoo), on a Core i7 4700MQ at 2.40GHz

``````import java.math.BigInteger;

public class PrimeSieveFactorialSchoenhage {

private static int[] primeList, multiList;

public static BigInteger factorial(int n) {
int log2n = 31 - Integer.numberOfLeadingZeros(n);
int piN = log2n < 2 ? 1 : 2 + (15 * n) / (8 * (log2n - 1));

primeList = new int[piN];
multiList = new int[piN];

int len = primeFactors(n);
return nestedSquare(len).shiftLeft(n - Integer.bitCount(n));
}

private static BigInteger nestedSquare(int len) {
if (len == 0) {
return BigInteger.ONE;
}

int i = 0, mult = multiList[0];

while (mult > 1) {
if ((mult & 1) == 1) { // is mult odd ?
primeList[len++] = primeList[i];
}

multiList[i++] = mult / 2;
mult = multiList[i];
}
BigInteger ns = nestedSquare(i);
if (len <= i) {
return ns.multiply(ns);
}

return product(primeList, i, len - i).multiply(ns.multiply(ns));
}

private static BigInteger product(int[] a, int start, int length) {
if (length == 0) {
return BigInteger.ONE;
}

int len = (length + 1) / 2;
long[] b = new long[len];

int i, j, k;

for (k = 0, i = start, j = start + length - 1; i < j; i++, k++, j--) {
b[k] = a[i] * (long) a[j];
}

if (i == j) {
b[k++] = a[j];
}

return recProduct(b, 0, k - 1);
}

private static BigInteger recProduct(long[] s, int n, int m) {
if (n > m) {
return BigInteger.ONE;
}
if (n == m) {
return BigInteger.valueOf(s[n]);
}
int k = (n + m) >> 1;
return recProduct(s, n, k).multiply(recProduct(s, k + 1, m));
}

private static int primeFactors(int n) {
int[] primes = new int[n < 17 ? 6 : (int) Math.floor(n / (Math.log(n) - 1.5))];
int numPrimes = makePrimeList(n, primes);

int maxBound = n / 2, count = 0;

int start = indexOf(primes, 2, 0, numPrimes - 1);
int end = indexOf(primes, n, start, numPrimes);

for (int i = start; i < end; i++) {
int prime = primes[i];
int m = prime > maxBound ? 1 : 0;

if (prime <= maxBound) {
int q = n;
while (q >= prime) {
m += q /= prime;
}
}

primeList[count] = prime;
multiList[count++] = m;
}
return count;
}

private static int indexOf(final int[] data, int value, int low, int high) {
while (low < high) {
int mid = (low + high) >>> 1;

if (data[mid] < value) {
low = mid + 1;
} else {
high = mid;
}
}

if (low >= data.length) {
return low;
}

if (data[low] == value) {
low++;
}

return low;
}

private static int makePrimeList(int n, int[] prime) {
boolean[] composite = new boolean[n / 3];

sieveOfEratosthenes(composite);

boolean toggle = false;
int p = 5, i = 0, j = 2;

prime[0] = 2;
prime[1] = 3;

while (p <= n) {
if (!composite[i++]) {
prime[j++] = p;
}
// -- never mind, it's ok.
p += (toggle = !toggle) ? 2 : 4;
}

return j; // number of primes
}

private static void sieveOfEratosthenes(final boolean[] composite) {
int d1 = 8;
int d2 = 8;
int p1 = 3;
int p2 = 7;
int s1 = 7;
int s2 = 3;
int n = 0;
int len = composite.length;
boolean toggle = false;

while (s1 < len) { // -- scan sieve
if (!composite[n++]) { // -- if a prime is found, cancel its multiples
int inc = p1 + p2;

for (int k = s1; k < len; k += inc) {
composite[k] = true;
}

for (int k = s1 + s2; k < len; k += inc) {
composite[k] = true;
}
}

if (toggle = !toggle) { // Never mind, it's ok.
s1 += d2;
d1 += 16;
p1 += 2;
p2 += 2;
s2 = p2;
} else {
s1 += d1;
d2 += 8;
p1 += 2;
p2 += 6;
s2 = p1;
}
}
}

public static void main(String[] args) {
int n = Integer.parseInt(args[0]);
long nanos = System.nanoTime();
BigInteger fact = factorial(n);
nanos = System.nanoTime() - nanos;
// Commented out because it takes ages to print
//System.out.println(fact);
System.out.println(nanos / 1e9);
}
}
``````
-
Compiled with `gcj -O3 --main=PrimeSieveFactorialSchoenhage PrimeSieveFactorialSchoenhage.java -o pf_nest_square_fact` – 14mRh4X0r May 12 '14 at 15:28

Julia - Score = 15.194

Utilising the exact same approach as that of the reference program... that is,

``````f(n)=reduce(*,1:big(n))
``````

So it uses reduce, the basic binary multiplication operation, and a range (in this case, using big(n) to force the calculation to be done in BigInt rather than Int64). From this, I get

``````julia> @time K=f(2340000);
elapsed time: 9.991324093 seconds (814552840 bytes allocated)
``````

On my computer, with reference program running with input of 154000, I get `Time: 10.041181 sec` output (run using `python ./test.py`, where test.py is the name of the file containing the reference code)

-

Mathematica 10^7

This is the plain vanilla, straight-out-of-the-box approach that shows typical performance for a particular processor. No points for creativity.

Ten million factorial takes a bit over 4 seconds on my MacBook (2.6GHz, i7 core). The result is roughly 65 million digits in length.

``````n=10000000! // Timing
``````

-
I think this counts as an external library. – Ypnypn May 6 '14 at 13:29
@Ypnypn if anything this is a built-in language feature. There's nothing "external" about it. I agree that using built-in features may be questionable but they haven't been ruled out here. – Martin Ender May 6 '14 at 13:32
@m.buettner I know it's not really an external library, but the point is that the answerer isn't doing the calculation and optimization himself, instead letting Mathmatica do the hard work, which has the same effect as using an external library. – Ypnypn May 6 '14 at 13:36
Since standard loopholes are mentioned in the OP, there is this: meta.codegolf.stackexchange.com/a/1078/14215 – Geobits May 6 '14 at 13:44
This is definitely a standard loophole, and also almost certainly uses multi-threading... – alexander-brett May 6 '14 at 15:26