Ile można szybko pomnożyć?

12

W związku z ostatnim uderzeniem w Python , oto próba pokazania mocnych stron Pythona. Twoim wyzwaniem jest napisanie programu, który oblicza silnię tak dużej liczby, jak to możliwe, w ciągu 10 sekund.n

Twój wynik będzie (highest n for your program on your machine)/(highest n for my program on your machine)

Zasady

  • Musisz obliczyć dokładne rozwiązanie liczb całkowitych. Ponieważ silnia byłaby znacznie wyższa niż to, co może zmieścić się w 64-bitowej liczbie całkowitej bez znaku, możesz użyć ciągów, jeśli twój język nie obsługuje dużych liczb całkowitych
  • Standardowe luki są zabronione. W szczególności nie można używać żadnych zasobów zewnętrznych.
  • Tylko część obliczeniowa (obejmuje czas na wszelkie obejścia wykorzystujące ciągi) dodaje do całkowitego czasu, który powinien być średnio poniżej 10 sekund.
  • Tylko programy jednowątkowe.
  • Musisz przechowywać dane wyjściowe w formie łatwej do wydrukowania (ponieważ drukowanie wymaga czasu) (patrz mój program poniżej), ciąg znaków, zmienna, tablica znaków itp.

EDYTOWAĆ:

  • Twój program musi dawać prawidłowe dane wyjściowe dla wszystkich n:1 <= n <= (your highest n)

EDYCJA 2:


Mój 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()
answer = factorial( 90000 )
end = time.clock()

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

Najwyższy wynik wygrywa. Dla przypomnienia, mój kod może zarządzać n = 90000w około 9.89sekundach na Pentium 4 3,0 GHz


EDYCJA: Czy każdy może dodać wynik, a nie tylko najwyższą liczbę n . Po prostu najwyższy nnie ma żadnego znaczenia, ponieważ zależy od twojego sprzętu. W przeciwnym razie nie można mieć obiektywnego kryterium wygranej. Odpowiedź ali0sha robi to poprawnie.


Mamy zwycięzcę. Nie zaakceptowałem odpowiedzi java /codegolf//a/26974/8766, ponieważ to rodzaj spódnic blisko http://meta.codegolf.stackexchange.com/a/1080/8766

użytkownik80551
źródło
1
Możesz użyć operator.mulzamiast funkcji lambda
gnibbler
1
Nieco zaskoczony, że to działa, ale przy założeniu, że poprawnie przeczytałem zasady, to rozwiązanie MATLAB byłoby fajne: factorial(Inf)zwraca Infw ułamku sekundy.
Dennis Jaheruddin
1
@Doorknob Pasuje do standardowych luk.
Justin
1
@DennisJaheruddin, określenie „Inf” jako „dokładnego rozwiązania liczb całkowitych” jest trochę skomplikowane.
tobyink
1
@Quincunx Nie, każdy język jest dozwolony.
user80551

Odpowiedzi:

7

C ++ z GMP, wynik = 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;
}
Keith Randall
źródło
8

Python 2.7

42,575 = (6,812 000/160 000) ok


Kod:

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

Wynik:

10.0069999695
10.0729999542
10.0360000134
9.98699998856

Jak to działa:

Większe zwielokrotnienia zajmują więcej czasu, dlatego chcemy wykonać jak najwięcej małych zwielokrotnień. Jest to szczególnie prawdziwe w Pythonie, gdzie dla liczb mniejszych niż 2^64używamy arytmetyki sprzętowej, a powyżej używamy oprogramowania. Tak więc, m(L)zaczynamy od listy L; jeśli ma nieparzystą długość, usuwamy jedną liczbę z rozważań, aby była jeszcze większa. Następnie mnożymy element 1przez element -2, element za 3pomocą -4itd., Aby

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]

Takie podejście zapewnia, że ​​używamy arytmetyki sprzętowej tak długo, jak to możliwe, po czym przełączamy się na wydajną bibliotekę arytmetyczną gmc.

W fac2, przyjmujemy również bardziej klasyczne podejście dziel i zwyciężaj, w którym dzielimy każdą wielokrotność 2 i przesuwamy je pod koniec bitów, aby uzyskać trywialny wzrost wydajności. Umieściłem go tutaj, ponieważ zwykle jest o około 0,5% szybszy niż fac1.

Wersja w golfa fac1(bo mogę), 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]
Alexander-Brett
źródło
1
Jeśli zaplecze GMP zawiera funkcję przesunięcia bitów, możesz zmniejszyć liczby, dzieląc każdą liczbę na liście przez 2, aż będzie parzysta, a następnie wykonując jedną zmianę na końcu.
Peter Taylor
Skąd masz gmpy2? $ python Python 2.7.3 (domyślnie, 27 lutego 2014, 19:58:35) [GCC 4.6.3] na linux2 Wpisz „pomoc”, „prawo autorskie”, „kredyty” lub „licencja”, aby uzyskać więcej informacji. >>> z gmpy2 import mpz Traceback (ostatnie połączenie ostatnio): Plik „<stdin>”, wiersz 1, w <module> ImportError: Brak modułu o nazwie gmpy2 >>>
user80551
@ user80551: code.google.com/p/gmpy (najlepszy wynik wyszukiwania Google) ma instalatory dla wielu różnych platform.
Alexander-Brett
W przypadku wersji golfowej nie mógłbyś while len(L): ...tego zrobić while len(L)>1: ...?
user80551
Nie: funkcja wewnątrz tej pętli nigdy nie zajmie listy poniżej długości 1 i tak czy inaczej potrzebujemy pierwszego elementu!
Alexander-Brett
2

Java - 125,15 (21 400 000 / 171,000)

Również bezwstydnie skopiowane z repozytorium Github Petera Luschnego (dzięki @ semi-extrinsic) i licencjonowane na licencji MIT, wykorzystuje algorytm „pierwiastkowania zagnieżdżonego do kwadratu”, jak zaproponowali Albert Schönhage i in. (zgodnie ze stroną opisu algorytmów czynnikowych Luschnego ).

Lekko dostosowałem algorytm, aby używał BigInteger Javy i nie używał tabeli odnośników dla n <20.

Skompilowany z gcj, który używa GMP do implementacji BigInteger i działał na Linuksie 3.12.4 (Gentoo), na Core i7 4700MQ przy 2,40 GHz

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);
    }
}
14mRh4X0r
źródło
Skompilowano zgcj -O3 --main=PrimeSieveFactorialSchoenhage PrimeSieveFactorialSchoenhage.java -o pf_nest_square_fact
14mRh4X0r
1

Python 3, n = 100000

Wystarczyła prosta zmiana algorytmu, aby podnieść kod przykładowy o 10000.

import time

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

start = time.clock()
answer = factorial(100000)
end = time.clock()

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

Oczywiście nie jest to najbardziej kreatywna odpowiedź, ale tak naprawdę jest tylko jeden sposób, aby zrobić czynnik.

Klamka
źródło
Proszę podać wynik, zobacz moją edycję. Guz będzie prawdopodobnie dlatego, że urządzenie jest lepsze niż moje.
user80551
1

Perl + C, n = około 3 milionów

Tutaj używam biblioteki Math :: BigInt :: GMP dostępnej na CPAN, która zapewnia ogromne przyspieszenie podstawowych obiektów Perla Math :: BigInt.

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 $answer;
say "Time: ", $end - $start, " sec";

Pamiętaj, że mój komputer jest prawdopodobnie nieco wolniejszy od twojego. Używając twojego oryginalnego skryptu Python, mogę obliczyć tylko factorial(40000)w 10 sekund; factorial(90000)zajmuje dużo dłużej. (Wcisnąłem Ctrl + C po minucie.) Na twoim sprzęcie, używając Math :: BigInt :: GMP, możesz być w stanie obliczyć silnię 5 milionów lub więcej w mniej niż 10 sekund.

Jedną z rzeczy, które możesz zauważyć, jest to, że chociaż silnia jest obliczana niewiarygodnie szybko, wydrukowanie wyniku jest bardzo wolne, zajmuje około trzy razy więcej niż pierwotne obliczenie. Wynika to z tego, że GMP wewnętrznie używa reprezentacji binarnej zamiast dziesiętnej, a wydrukowanie go wymaga konwersji binarnej na dziesiętną.

Tobyink
źródło
1
Myślę, że GMP liczy się jako zasób zewnętrzny. (Chociaż z pewnością sprawia to o wiele łatwiejsze niż wdrażanie faktoryzacji pierwszej i mnożenie Schönhage-Strassen od zera.)
r3mainer
3
Zakładałem, że „zasób zewnętrzny” odnosi się do wyszukiwania rozwiązań z wcześniej obliczonego zestawu odpowiedzi w bazie danych, serwisie internetowym itp.
tobyink
Squeamish: biblioteki zwykle nie liczą się jako zasoby zewnętrzne, chyba że mają funkcję, która wchodzi w zakres nudnej luki.
Alexander-Brett
1
Tobyink: czy możesz wyjaśnić, co robi twój program? wygląda na to, że używasz tylko wbudowanej funkcji (bfac?)
Alexander-Brett
Tak. Ta odpowiedź jest nieprawidłowa, ponieważ wykorzystuje metodę silni zMath::BigInt
14mRh4X0r
1

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)

Wykorzystuje względną łatwość mnożenia wielu grup małych liczb, a następnie mnożenia ich w porównaniu do dużej liczby mnożeń z udziałem dużej liczby.

Cthulhu
źródło
1

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

Nie jestem z tego zadowolony.

Czy C # jest tak wolny?

Ale i tak jest mój wpis.

static void Main(string[] args)
    {
        Console.WriteLine("Enter N for fatorial:");
        int n = Convert.ToInt32(Console.ReadLine());

        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);

    }

Gdy n = 77000 trzeba 00:00:09:8708952obliczyć.

Pracuję w trybie Release, poza Visual Studio, używając Core i3-2330M @ 2.2GHz.

Edycja: Ponieważ nie robię niczego inteligentnego, akceptuję ten wynik. Być może .NET Framework 4.5 dodaje trochę narzutu (lub BigInteger nie jest tak szybki).

Ricardo Pieper
źródło
Podaj wynik, a nie tylkon
użytkownik80551
1
Możesz użyć zero approached byoperatora, aby był ładniejszy (na przykład zacznij od, n = ... + 1a potem zrób while (0 <-- n) result *= n;)
Cthulhu
1
BigInteger dla .NET prawdopodobnie nie zaimplementował algorytmów do pomnażania większych liczb, takich jak Karatsuba lub Toom-3. Jeśli tak, to dobry przykład tego, jak Python jest szybszy.
kernigh
1

bc, wynik = 0,19

Co do cholery, oto mój pretendent do „Ile możesz powoli pomnożyć?”

bc jest „językiem kalkulatora o dowolnej precyzji”, ale niestety raczej powolnym:

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

W około 10 sekund na moim MacBooku Pro z połowy 2012 roku (2,3 GHz Intel Core i7) odpowiedź na python referencyjny może obliczyć 122000 !, ale ten skrypt bc może obliczyć tylko 23600 !.

Z drugiej strony 10000! zajmuje 1.5s ze skryptem referencyjnym Pythona, ale skrypt bc zajmuje 50s.

O jej.

Cyfrowa trauma
źródło
1
OpenBSD bc (1) jest szybszy. Twój program uzyskał wynik 0,29 = 28000/98000. Nie ma read(), więc pobiegłem time sed 's/read()/28000/' factorial.bc | bc.
kernigh
1

Bash: wynik = 0,001206 (181/150000)

Ukradłem funkcje matematyczne Rosettacode - długie mnożenie, którego nie analizowałem ani nie próbowałem zoptymalizować.
Możesz zmienić algorytm lub wypróbować inną metodę podziału ciągów.

#!/bin/bash


add() { # arbitrary-precision addition
  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"
    product="$(add "$product" "$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 
Emmanuel
źródło
1
Dodaj wynik, a nie tylko najwyższy n
user80551
@ user80551 gotowe
Emmanuel
1

Python 3, Advanced Algo autor: Peter Luschny: 8,25x (1 280 000/155 000)

Bezwstydnie skopiowane od Petera Luschnego,
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm ,
który udostępnia ten kod na licencji „Creative Commons Uznanie autorstwa-Na tych samych warunkach 3.0”.

Jest to właściwie dość zaawansowany algorytm, wykorzystujący coś, co nazywa się „zmiennym silnikiem” i listą liczb pierwszych. Podejrzewam, że mogłoby być jeszcze szybciej, gdyby podobało się wielu innym odpowiedziom i wykonał większość multiplikacji z 32-bitowymi liczbami całkowitymi.

#! /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()
answer = factorialPS(1280000) 
print(time.time()-start)
pół-zewnętrzny
źródło
1

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 {
            BigInteger mid = start.add(end.subtract(start).divide(Factorials.two));
            return Factorials.mul(start, mid).multiply(Factorials.mul(mid.add(BigInteger.ONE), end));
        }
    }

    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()+" @ "+
    }
}

BigIntegers są powolne.

Zalecenia dotyczące bibliotek liczb całkowitych Java o dowolnej precyzji? : P

cjfaure
źródło
Czy mogę ukraść twój kod, aby był wielowątkowy?
Simon Kuang
@SimonKuang Śmiało. : P Wpisy wielowątkowe nie są tu jednak dozwolone. Możesz także użyć bardziej wydajnej implementacji BigInteger.
cjfaure
Mergesort-y Nazywa się to „dziel i rządź”.
johnchen902
1

C ++ (specyficzny dla x86_64) - 3.0 (390000/130000)

(łatwo przenośny na x86-32, przenoszenie na inne architektury oznacza znaczną utratę prędkości)

Oto moja własna mikro-implementacja długiej arytmetyki.
Samo obliczenie zajmuje 10 sekund, a mimo że wydruk jest łatwy do wydrukowania (zobacz operator<<przeciążenie), wydrukowanie go zajmuje więcej czasu.

#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);
        *i = dadd(res, carry, carry);
        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;
}
mniip
źródło
Sprawdź swój wynik. Musisz uruchomić program Python 2.7 na swoim komputerze. Na moim komputerze skompilowałem twój program g++ -O2 factorial.cc -o factoriali jego wynik wynosi 3,90 = 382000 / 98000.
kernigh
Dziwne, mam 3,9, a ty 3.0 dla tego programu. Chyba twój szybszy komputer to kara. Być może twój program traci przewagę nad Pythonem w miarę rwzrostu. Jeśli tak, i możesz zrobić wyższą rw 10 sekund, twój wynik spada.
kernigh
0

Python 3: 280000/168000

Czas działania programu: pomiędzy 9.87585953253i 10.3046453994. Czas działania mojego programu: około 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)

Przeczytałem tę odpowiedź na cs.SE i postanowiłem spróbować ją zaimplementować w Pythonie. Jednak przypadkowo odkryłem, że n! = (⌊n / 2⌋)! * 2**(⌊n / 2⌋) * n!!(uwaga: !!jest podwójnym silnikiem ). Przekształciłem to w formę nierekurencyjną.

Komentarze pokazują moją próbę uniknięcia ponownego obliczenia podwójnej silni, ale odkryłem, że przechowywanie każdej wartości było zbyt kosztowne dla pamięci, co spowodowało, że mój komputer działał jeszcze wolniej. Mogę to poprawić, przechowując tylko to, co jest potrzebne.

O dziwo zaimplementowałem naiwne proste mnożenie w Pythonie 3 i robi to lepiej niż twój program: n = 169000w 10 sekund .:

def factorial(n):
    p=1
    for i in range(n):
        p*=i+1
    return p
Justin
źródło
0

Ruby 2.1

wynik = 1,80 = 176_000 / 98_000

EDYCJA: poprawiona z 1.35 = 132_000 / 98_000

Wziąłem pomysły z algorytmu czynnikowego GMP . Ten program używa standardowej biblioteki do generowania liczb pierwszych. Ruby to zły wybór, ponieważ mnożenie wydaje się wolniejsze w Rubim niż w Pythonie.

  1. Mój program w Ruby 2.1: wynik = 1,80 = 176_000 / 98_000
  2. Trywialny algorytm w Pythonie 2.7: wynik = 1 = 98_000 / 98_000
  3. Trywialny algorytm w Ruby 2.1: wynik = 0,878 = 86_000 / 98_000

Tak, mój plik binarny ruby 2.1.0p0 (2013-12-25 revision 44422) [x86_64-openbsd]linków przeciwko GMP. Ruby 2.1 dodał funkcję do korzystania z GMP do dużego mnożenia, ale nadal wydaje się wolniejszy niż 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"
kernigh
źródło
0

Julia - wynik = 15,194

Wykorzystując dokładnie to samo podejście, co w programie referencyjnym ...

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

Wykorzystuje więc redukcję, podstawową operację mnożenia binarnego i zakres (w tym przypadku użycie dużej (n) do wymuszenia wykonania obliczeń w BigInt zamiast Int64). Z tego rozumiem

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

Na moim komputerze, z programem referencyjnym działającym z wejściem 154000, otrzymuję Time: 10.041181 secdane wyjściowe (uruchom przy użyciu python ./test.py, gdzie test.py to nazwa pliku zawierającego kod referencyjny)

Glen O
źródło
0

tcl, 13757

Moja odpowiedź to sprawdzić limity tcl.

Pierwszy wiersz służy tylko do ustawienia parametru wejściowego:

set n 13757

Pozostałe to sam algorytm:

set r 2
for {set i 3} {$i <= $n} {incr i} {set r [expr {$r*$i}]}   
puts $r

Przetestowałem swój kod na http://rextester.com/live/WEL36956 ; Jeśli powiększę n, dostanę SIGKILL; may n może stać się większy na lokalnym tłumaczu tcl, którego nie mam.

sergiol
źródło