﻿ FFF CompLang Shootout
```/// Author & Copyright for Scala and C# code: Peter Luschny
/// Author & Copyright for GO code: Sonia Keys
/// License: MIT licence (or at your option for Scala and C#)
 SCALA C# GO SWING SCALA SWING C# SWING GO `class Swing(n: Int) {` ` ` ` // OEIS: A056040 Swing numbers` ` def name() = "Swing"` ` ` ` private var primes: Primes = null` ` private var factors: Array[Long] = null` ` ` ` if (n >= Small.oddSwing.length)` ` {` ` primes = new Primes(n)` ` factors = new Array[Long](n)` ` }` ` ` ` def value(m: Int = n):BigInt =` ` {` ` return oddSwing(m) << MathFun.bitCount(m)` ` }` ` ` ` def oddSwing(k: Int): BigInt =` ` {` ` if (k > n)` ` {` ` throw new IllegalArgumentException(` ` "k has to be >= n, but was " + k)` ` }` ` ` ` if (k < Small.oddSwing.length)` ` return BigInt(Small.oddSwing(k));` ` ` ` val rootK = MathFun.floorSqrt(k)` ` var i = 0` ` ` ` primes.iterator(3, rootK, p =>` ` {` ` var q = k / p` ` while (q > 0)` ` {` ` if ((q & 1) == 1) {factors(i) = p; i += 1}` ` q /= p` ` }` ` })` ` ` ` primes.iterator(rootK + 1, k / 3, p =>` ` {` ` if (((k / p) & 1) == 1) {factors(i) = p; i += 1}` ` })` ` ` ` primes.iterator(k / 2 + 1, k, p =>` ` {` ` factors(i) = p; i += 1` ` })` ` ` ` return MathFun.product(factors,0,i)` ` }` `}` ` public class Swing` ` {` ` public string Name` ` {` ` get { return "Swing"; }` ` }` ` private uint n;` ` private Primes primes;` ` private Factors factors;` ` ` ` // OEIS: A056040 Swinging factorial` ` public Swing(uint n)` ` {` ` if (n >= Small.OddSwing.Length)` ` {` ` primes = new Primes(n);` ` factors = new Factors(n);` ` }` ` this.n = n;` ` }` ` ` ` public BigInt Value()` ` {` ` return OddSwing(n) << (int)MathFun.BitCount(n);` ` }` ` ` ` public BigInt OddSwing(uint k)` ` {` ` if (k > n)` ` {` ` throw new ArgumentOutOfRangeException(` ` "k has to be >= n, but was " + k);` ` }` ` if (k < Small.OddSwing.Length) return Small.OddSwing[k];` ` ` ` uint rootN = MathFun.FloorSqrt(k);` ` factors.Init();` ` ` ` factors.SetMax(rootN);` ` primes.Iterator(3, rootN, p =>` ` {` ` uint q = k;` ` while ((q /= p) > 0)` ` if ((q & 1) == 1) { factors.Add(p); }` ` });` ` ` ` factors.SetMax(k / 3);` ` primes.Iterator(rootN + 1, k / 3, p =>` ` {` ` if (((k / p) & 1) == 1) { factors.Add(p); }` ` });` ` ` ` factors.SetMax(k);` ` primes.Iterator(k / 2 + 1, k, p =>` ` {` ` factors.Add(p);` ` });` ` ` ` return factors.Product();` ` }` ` }` `type swing struct {` ` primes *primes` ` factors []uint64` `}` ` ` `func makeSwing(n uint) (ps *swing) {` ` ps = new(swing)` ` ps.primes = makePrimes(uint64(n))` ` ` ` if n >= uint(len(smallOddSwing)) {` ` ps.factors = make([]uint64, n)` ` }` ` ` ` return` `}` ` ` `func (ps *swing) swing(m uint) *big.Int {` ` if uint64(m) > ps.primes.sieveLen {` ` return nil` ` }` ` r := ps.oddSwing(m)` ` return r.Lsh(r, bitCount(uint64(m)))` `}` ` ` `func (ps *swing) oddSwing(k uint) *big.Int {` ` if k < uint(len(smallOddSwing)) {` ` return big.NewInt(smallOddSwing[k])` ` }` ` ` ` rootK := floorSqrt(k)` ` var i int` ` ` ` ps.primes.iterator(3, uint64(rootK), func(p uint64) {` ` q := uint64(k) / p` ` for q > 0 {` ` if q & 1 == 1 {` ` ps.factors[i] = p` ` i++` ` }` ` q /= p` ` }` ` })` ` ` ` ps.primes.iterator(uint64(rootK+1),uint64(k/3),func(p uint64){` ` if (uint64(k) / p & 1) == 1 {` ` ps.factors[i] = p` ` i++` ` }` ` })` ` ` ` ps.primes.iterator(uint64(k/2 + 1), uint64(k), func(p uint64){` ` ps.factors[i] = p` ` i++` ` })` ` ` ` return product(ps.factors[0:i])` `}` PRIME FACTORIAL SCALA PRIME FACTORIAL C# PRIME FACTORIAL GO `class PrimeFactorial(n: Int) {` ` ` ` def name() = "PrimeFactorial"` ` var swing: Swing = new Swing(n)` ` ` ` def value(k : Int = n): BigInt =` ` {` ` if (k > n)` ` {` ` throw new IllegalArgumentException(` ` "k has to be >= n, but was " + k)` ` }` ` ` ` val exp2 = k - MathFun.bitCount(k)` ` if (k < Small.oddFactorial.length)` ` return BigInt(Small.oddFactorial(k)) << exp2` ` ` ` return oddFactorial(k) << exp2` ` }` ` ` ` private def oddFactorial(n: Int): BigInt =` ` {` ` if (n < Small.oddFactorial.length)` ` {` ` return BigInt(Small.oddFactorial(n))` ` }` ` ` ` val of = oddFactorial(n / 2)` ` return (of * of) * swing.oddSwing(n)` ` }` `}` ` public class PrimeFactorial` ` {` ` public string Name` ` {` ` get { return "PrimeFactorial"; }` ` }` ` ` ` public PrimeFactorial() { }` ` ` ` private Swing swing;` ` ` ` // OEIS: A000142 Factorial` ` public BigInt Value(int n)` ` {` ` if (n < 0)` ` {` ` throw new ArgumentOutOfRangeException(` ` "Factorial: n has to be >= 0, but was " + n);` ` }` ` ` ` int exp2 = (int)(n - MathFun.BitCount((uint)n));` ` if (n < Small.OddFactorial.Length)` ` return new BigInt(Small.OddFactorial[n]) << exp2;` ` ` ` swing = new Swing((uint)n);` ` return OddFactorial((uint)n) << exp2;` ` }` ` ` ` private BigInt OddFactorial(uint n)` ` {` ` if (n < Small.OddFactorial.Length)` ` {` ` return Small.OddFactorial[n];` ` }` ` ` ` return BigInt.Pow(OddFactorial(n / 2), 2) * swing.OddSwing(n);` ` }` ` } ` `func (ps *swing) primeFactorial(k uint) (r *big.Int) {` ` if uint64(k) > ps.primes.sieveLen {` ` return nil` ` }` ` r = ps.oddFactorial(k)` ` return r.Lsh(r, k - bitCount(uint64(k)))` `}` ` ` `func (ps *swing) oddFactorial(n uint) *big.Int {` ` if n < uint(len(smallOddFactorial)) {` ` return big.NewInt(int64(smallOddFactorial[n]))` ` }` ` ` ` of := ps.oddFactorial(n / 2)` ` return of.Mul(of.Mul(of, of), ps.oddSwing(n))` `}` ` ` ```val start2 = System.nanoTime var ps = new PrimeFactorial(n) val v = ps.value(n) val stop2 = System.nanoTime val sec2 = DoubleFormat.fix2((stop2 - start2)/1e9) + " sec" n = 234567 -> 68,5 sec ``` ```var watch = new System.Diagnostics.Stopwatch(); watch.Reset(); watch.Start(); var pFact = new PrimeFactorial(); pFact.Value(n); watch.Stop(); Console.WriteLine(pFact.Name + " " + watch.Elapsed); n = 234567 -> PrimeFactorial 00:00:38.6056398 ``` `func main() {` ` const n = 234567` ` start := time.Nanoseconds()` ` fs := makeSwing(n)` ` a := fs.primeFactorial(n)` ` stop := time.Nanoseconds()` ` fmt.Printf("n = %d -> PrimeFactorial %.2f sec\n",` ` n, float64(stop - start) / 1e9)` ` ` ` dtrunc := int64(float64(a.BitLen()) * .30103) - 10` ` var first, rest big.Int` ` rest.Exp(first.SetInt64(10), rest.SetInt64(dtrunc), nil)` ` first.Quo(a, &rest)` ` fstr := first.String()` ` fmt.Printf("%d! begins %s... and has %d digits.\n",` ` n, fstr, int64(len(fstr)) + dtrunc)` `}` SWING FACTORIAL SCALA SWING FACTORIAL C# SWING FACTORIAL GO `class SwingFactorial(n: Int) {` ` ` ` def name() = "SwingFactorial"` ` ` ` def value(): BigInt =` ` {` ` if (n < 0)` ` {` ` throw new IllegalArgumentException(` ` "Factorial: n has to be >= 0, but was " + n)` ` }` ` ` ` ndiv2OddFact = BigInt(1)` ` ndiv4OddFact = ndiv2OddFact` ` ` ` return oddFactorial(n) << (n - MathFun.bitCount(n))` ` }` ` ` ` private def oddFactorial(n: Int): BigInt =` ` {` ` val oddFact =` ` if (n < Small.oddFactorial.length)` ` {` ` BigInt(Small.oddFactorial(n))` ` }` ` else` ` {` ` val of = oddFactorial(n / 2)` ` (of * of) * oddSwing(n)` ` }` ` ` ` ndiv4OddFact = ndiv2OddFact` ` ndiv2OddFact = oddFact` ` return oddFact` ` }` ` ` ` private def oddSwing(n: Int): BigInt =` ` {` ` if (n < Small.oddSwing.length)` ` {` ` return BigInt(Small.oddSwing(n))` ` }` ` ` ` val len = if ((n % 4) != 2) (n - 1) / 4 + 1 else (n - 1) / 4` ` val high = n - ((n + 1) & 1)` ` val ndiv4 = n / 4` ` val oddFact = if (ndiv4 < Small.oddFactorial.length)` ` BigInt(Small.oddFactorial(ndiv4)) else ndiv4OddFact` ` ` ` return product(high, len) / oddFact` ` }` ` ` ` private def product(m: Int, len: Int): BigInt =` ` {` ` if (len == 1) return BigInt(m) ` ` if (len == 2) {val M = m.toLong; return BigInt(M * (M - 2))}` ` ` ` val hlen = len >>> 1` ` return product(m - hlen * 2, len - hlen) * product(m, hlen)` ` }` ` ` ` private var ndiv4OddFact = BigInt(1)` ` private var ndiv2OddFact = BigInt(1)` `} ` ` public class SwingFactorial` ` {` ` public string Name` ` {` ` get { return "SwingFactorial"; }` ` }` ` ` ` public SwingFactorial() { }` ` ` ` private BigInt oddFactNdiv4, oddFactNdiv2;` ` private int SMALLSWING = Small.OddSwing.Length;` ` private int SMALLFACT = Small.OddFactorial.Length;` ` ` ` public BigInt Value(int n)` ` {` ` if (n < 0)` ` {` ` throw new ArgumentOutOfRangeException(` ` "Factorial: n has to be >= 0, but was " + n);` ` }` ` oddFactNdiv4 = oddFactNdiv2 = BigInt.One;` ` return OddFactorial(n) << (n-(int)MathFun.BitCount((uint)n));` ` }` ` ` ` private BigInt OddFactorial(int n)` ` {` ` BigInt oddFact;` ` ` ` if (n < SMALLFACT)` ` {` ` oddFact = Small.OddFactorial[n];` ` }` ` else` ` {` ` BigInt sqrOddFact = OddFactorial(n / 2);` ` oddFact = BigInt.Pow(sqrOddFact, 2) * OddSwing(n);` ` }` ` ` ` oddFactNdiv4 = oddFactNdiv2;` ` oddFactNdiv2 = oddFact;` ` return oddFact;` ` }` ` ` ` private BigInt OddSwing(int n)` ` {` ` if (n < SMALLSWING) return Small.OddSwing[n];` ` ` ` int len = (n - 1) / 4;` ` if ((n % 4) != 2) len++;` ` int high = n - ((n + 1) & 1);` ` int ndiv4 = n / 4;` ` BigInt oddFact = ndiv4 < SMALLFACT ? ` ` Small.OddFactorial[ndiv4] : oddFactNdiv4;` ` ` ` return Product(high, len) / oddFact;` ` }` ` ` ` private static BigInt Product(int m, int len)` ` {` ` if (len == 1) return new BigInt(m);` ` if (len == 2) return new BigInt((long)m * (m - 2));` ` ` ` int hlen = len >> 1;` ` return Product(m - hlen * 2, len - hlen) * Product(m, hlen);` ` }` ` }` `func swingFactorial(n uint) (r *big.Int) {` ` var oddFactNDiv2, oddFactNDiv4 big.Int` ` ` ` // closes on oddFactNDiv2, oddFactNDiv4` ` oddSwing := func(n uint) (r *big.Int) {` ` if n < uint(len(smallOddSwing)) {` ` return big.NewInt(smallOddSwing[n])` ` }` ` ` ` length := (n - 1) / 4` ` if n%4 != 2 {` ` length++` ` }` ` high := n - (n+1)&1` ` ndiv4 := n / 4` ` ` ` var oddFact big.Int` ` if ndiv4 < uint(len(smallOddFactorial)) {` ` oddFact.SetInt64(smallOddFactorial[ndiv4])` ` r = &oddFact` ` } else {` ` r = &oddFactNDiv4` ` }` ` ` ` return oddFact.Quo(oddProduct(high, length), r)` ` }` ` ` ` var oddFactorial func(uint) *big.Int` ` ` ` // closes on oddFactNDiv2, oddFactNDiv4, oddSwing, and itself` ` oddFactorial = func(n uint) (oddFact *big.Int) {` ` if n < uint(len(smallOddFactorial)) {` ` oddFact = big.NewInt(smallOddFactorial[n])` ` } else {` ` oddFact = oddFactorial(n / 2)` ` oddFact.Mul(oddFact.Mul(oddFact, oddFact), oddSwing(n))` ` }` ` ` ` oddFactNDiv4.Set(&oddFactNDiv2)` ` oddFactNDiv2.Set(oddFact)` ` return oddFact` ` }` ` ` ` oddFactNDiv2.SetInt64(1)` ` oddFactNDiv4.SetInt64(1)` ` r = oddFactorial(n)` ` return r.Lsh(r, n-bitCount(uint64(n)))` `}` ` ` `func oddProduct(m, length uint) *big.Int {` ` switch length {` ` case 1:` ` return big.NewInt(int64(m))` ` case 2:` ` var mb big.Int` ` mb.SetInt64(int64(m))` ` mb2 := big.NewInt(int64(m - 2))` ` return mb2.Mul(&mb, mb2)` ` }` ` hlen := length / 2` ` h := oddProduct(m-hlen*2, length-hlen)` ` return h.Mul(h, oddProduct(m, hlen))` `}` ```val start = System.nanoTime var fs = new SwingFactorial(n) val a = fs.value() val stop = System.nanoTime val sec = DoubleFormat.fix2((stop - start)/1e9) + " sec" n = 234567 -> 82,3 sec ``` ```watch.Reset(); watch.Start(); var sFact = new SwingFactorial(); sFact.Value(n); watch.Stop(); Console.WriteLine(sFact.Name + " " + watch.Elapsed); n = 234567 -> SwingFactorial 00:00:49.5329715 ``` `func testSwingFactorial() {` ` const n = 234567` ` start := time.Nanoseconds()` ` a := swingFactorial(n)` ` stop := time.Nanoseconds()` ` fmt.Printf("n = %d -> SwingFactorial %.2f sec\n",` ` n, float64(stop-start)/1e9)` ` ` ` dtrunc := int64(float64(a.BitLen())*.30103) - 10` ` var first, rest big.Int` ` rest.Exp(first.SetInt64(10), rest.SetInt64(dtrunc), nil)` ` first.Quo(a, &rest)` ` fstr := first.String()` ` fmt.Printf("%d! begins %s... and has %d digits.\n",` ` n, fstr, int64(len(fstr))+dtrunc)` `}` BINOMIAL SCALA BINOMIAL C# BINOMIAL GO `class Binomial(n: Int) {` ` ` ` val primes = new Primes(n)` ` var factors: Array[Long] = new Array[Long](n)` ` val rootN = MathFun.floorSqrt(n)` ` ` ` def value(j: Int): BigInt =` ` {` ` var k = j; var i = 0` ` if (0 > k || k > n) return BigInt(0)` ` if (k > n / 2) {k = n - k; }` ` if (k < 3)` ` {` ` if (k == 0) return BigInt(1)` ` if (k == 1) return BigInt(n)` ` val N = n.toLong` ` if (k == 2) return BigInt((N * (N - 1)) / 2)` ` }` ` ` ` primes.iterator(2, rootN, p =>` ` {` ` var r = 0; var N = n; var K = k` ` while (N > 0)` ` {` ` r = if ((N % p) < (K % p + r)) 1 else 0` ` if (r == 1) { factors(i) = p; i += 1 }` ` N /= p; K /= p` ` }` ` })` ` ` ` primes.iterator(rootN + 1, n / 2, p =>` ` {` ` if (n % p < k % p)` ` {` ` factors(i) = p; i += 1` ` }` ` })` ` ` ` primes.iterator(n - k + 1, n, p =>` ` {` ` factors(i) = p; i += 1` ` })` ` ` ` return MathFun.product(factors,0,i)` ` }` `}` ` public class Binomial` ` {` ` Primes prime;` ` Factors factors;` ` uint rootN;` ` uint n;` ` ` ` public Binomial(uint n)` ` {` ` prime = new Primes(n);` ` factors = new Factors(n);` ` rootN = MathFun.FloorSqrt(n);` ` this.n = n;` ` }` ` ` ` // OEIS: A007318 Pascal's triangle` ` public BigInt Value(uint k)` ` {` ` if (0 > k || k > n) return BigInt.Zero;` ` if (k > n / 2) { k = n - k; }` ` if (k < 3)` ` {` ` if (k == 0) return BigInt.One;` ` if (k == 1) return new BigInt(n);` ` if (k == 2) return new BigInt(((ulong)n * (n - 1)) / 2);` ` }` ` ` ` factors.Init();` ` ` ` factors.SetMax(rootN);` ` prime.Iterator(2, rootN, p =>` ` {` ` uint r = 0, N = n, K = k;` ` while (N > 0)` ` {` ` r = (N % p) < (K % p + r) ? 1u : 0u;` ` if (r == 1)` ` {` ` factors.Add(p);` ` }` ` N /= p; K /= p;` ` }` ` });` ` ` ` factors.SetMax(n / 2);` ` prime.Iterator(rootN + 1, n / 2, p =>` ` {` ` if (n % p < k % p)` ` {` ` factors.Add(p);` ` }` ` });` ` ` ` factors.SetMax(n);` ` prime.Iterator(n - k + 1, n, p =>` ` {` ` factors.Add(p);` ` });` ` ` ` return factors.Product();` ` }` ` } ` `func (p *primes) binomial(n, k uint) *big.Int {` ` if uint64(n) > p.sieveLen {` ` return nil` ` }` ` ` ` var r big.Int` ` if k > n {` ` return &r` ` }` ` ` ` if k > n/2 {` ` k = n - k` ` }` ` ` ` if k < 3 {` ` switch k {` ` case 0:` ` return r.SetInt64(1)` ` case 1:` ` return r.SetInt64(int64(n))` ` case 2:` ` var n1 big.Int` ` return r.Rsh(r.Mul(r.SetInt64(int64(n)), n1.SetInt64(int64(n-1))), 1)` ` }` ` }` ` ` ` var i int` ` rootN := uint64(floorSqrt(n))` ` factors := make([]uint64, n)` ` p.iterator(2, rootN, func(p uint64) {` ` var r, nn, kk uint64 = 0, uint64(n), uint64(k)` ` for nn > 0 {` ` if nn%p < kk%p+r {` ` r = 1` ` factors[i] = p` ` i++` ` } else {` ` r = 0` ` }` ` nn /= p` ` kk /= p` ` }` ` })` ` ` ` p.iterator(rootN+1, uint64(n/2), func(p uint64) {` ` if uint64(n)%p < uint64(k)%p {` ` factors[i] = p` ` i++` ` }` ` })` ` ` ` p.iterator(uint64(n-k+1), uint64(n), func(p uint64) {` ` factors[i] = p` ` i++` ` })` ` ` ` return product(factors[0:i])` `}` ```val start = System.nanoTime val binomial = new Binomial(n) for (k <- 33 until n by 10000) val expected = binomial.value(k) val stop = System.nanoTime val sec = DoubleFormat.fix2((stop - start)/1e9) + " sec " n = 234567 -> 5,45 sec ``` ```watch.Start(); var binomial = new Binomial(234567); for (uint k = 33; k < n; k += 10000) var expected = binomial.Value(k); watch.Stop(); n = 234567 -> 00:00:03.45 ``` `func testBinomial() {` ` const n = 234567` ` start := time.Nanoseconds()` ` p := makePrimes(n)` ` var a *big.Int` ` for k := uint(33); k <= n; k += 10000 {` ` b := p.binomial(n, k)` ` if k == 120033 {` ` a = b` ` }` ` }` ` stop := time.Nanoseconds()` ` fmt.Printf("n = %d -> %d Binomials %.2f sec\n",` ` n, (n-32)/10000, float64(stop-start)/1e9)` ` ` ` dtrunc := int64(float64(a.BitLen())*.30103) - 10` ` var first, rest big.Int` ` rest.Exp(first.SetInt64(10), rest.SetInt64(dtrunc), nil)` ` first.Quo(a, &rest)` ` fstr := first.String()` ` fmt.Printf("C(%d, 120033) begins %s... and has %d digits.\n",` ` n, fstr, int64(len(fstr))+dtrunc)` `}` DOUBLE-FACTORIAL SCALA DOUBLE-FACTORIAL C# DOUBLE-FACTORIAL GO `class DoubleFactorial(n: Int) {` ` ` ` def name() = "DoubleFactorial"` ` val N = if ((n & 1) == 0) n / 2 else n + 1` ` val swing: Swing = new Swing(N)` ` ` ` // OEIS: A006882 Double factorials n!!: a(n) = n*a(n-2)` ` def value: BigInt =` ` {` ` var dblFact = if (n < Small.oddDoubleFactorial.length)` ` BigInt(Small.oddDoubleFactorial(n)) else oddFactorial(N, n)` ` ` ` if ((n & 1) == 0)` ` {` ` val exp2 = n - MathFun.bitCount(n / 2)` ` dblFact = dblFact << exp2` ` }` ` ` ` return dblFact` ` }` ` ` ` def oddFactorial(n: Int, m: Int): BigInt =` ` {` ` if (n < Small.oddFactorial.length)` ` {` ` return BigInt(Small.oddFactorial(n))` ` }` ` ` ` var oddFact = oddFactorial(n / 2, m)` ` if (n < m) { oddFact *= oddFact }` ` ` ` return oddFact * swing.oddSwing(n)` ` }` `}` ` ` ` public class DoubleFactorial` ` {` ` private Swing swing;` ` ` ` public DoubleFactorial() { }` ` ` ` // OEIS: A006882 Double factorials n!!: a(n) = n*a(n-2).` ` public BigInt Value(uint n)` ` {` ` if (n < 0)` ` {` ` throw new ArgumentOutOfRangeException(` ` "DoubleFactorial: n has to be >= 0, but was " + n);` ` }` ` ` ` BigInt dblFact;` ` uint N = ((n & 1) == 0) ? n / 2 : n + 1;` ` ` ` if (n < Small.OddDoubleFactorial.Length)` ` {` ` dblFact = (BigInt)Small.OddDoubleFactorial[n];` ` }` ` else` ` {` ` if (N >= Small.OddSwing.Length)` ` {` ` swing = new Swing(N);` ` }` ` dblFact = OddFactorial(N, n);` ` }` ` ` ` if ((n & 1) == 0)` ` {` ` int exp2 = (int)(n - MathFun.BitCount(n / 2));` ` dblFact = dblFact << exp2;` ` }` ` return dblFact;` ` }` ` ` ` private BigInt OddFactorial(uint n, uint m)` ` {` ` if (n < Small.OddFactorial.Length)` ` {` ` return Small.OddFactorial[n];` ` }` ` ` ` var oddSwing = n < Small.OddSwing.Length ` ` ? Small.OddSwing[n] : swing.OddSwing(n); ` ` ` ` var oddFact = OddFactorial(n / 2, m);` ` if (n < m) { oddFact *= oddFact; }` ` ` ` return oddFact * oddSwing;` ` }` ` } ` `// returns nil if sieve not big enough` `func (p *primes) doubleFactorial(n uint) (r *big.Int) {` ` ` ` nEven := n&1 == 0` ` ` ` if n < uint(len(smallOddDoubleFactorial)) {` ` r = big.NewInt(smallOddDoubleFactorial[n])` ` } else {` ` var nn uint` ` if nEven {` ` nn = n / 2` ` } else {` ` nn = n + 1` ` }` ` ` ` if uint64(nn) > p.sieveLen && nn > uint(len(smallOddSwing)){` ` return nil` ` }` ` ` ` r = p.oddDoubleFactorial(nn, n)` ` }` ` ` ` if nEven {` ` r.Lsh(r, n - bitCount(uint64(n/2)))` ` }` ` ` ` return` `}` ` ` `func (p *primes) oddDoubleFactorial(n, m uint) *big.Int {` ` ` ` if n < uint(len(smallOddFactorial)) {` ` return big.NewInt(smallOddFactorial[n])` ` }` ` ` ` of := p.oddDoubleFactorial(n/2, m)` ` if n < m {` ` of.Mul(of, of)` ` }` ` ` ` return of.Mul(of, p.oddSwing(n))` `}` ```val start = System.nanoTime val df = new DoubleFactorial(n) val v = df.value val stop = System.nanoTime val sec = DoubleFormat.fix2((stop - start)/1e9) + " sec " n = 234566 -> 14,8 sec n = 234567 -> 20,6 sec ``` ```watch.Start(); var df = new DoubleFactorial(); var b = df.Value(n); watch.Stop(); Console.Write("primedouble " + watch.Elapsed + " "); n = 234566 -> 00:00:09.28 n = 234566 -> 00:00:13.19 ``` ` ` `func testDoubleFactorial() {` ` const nEven = 234566` ` const nOdd = 234567` ` start := time.Nanoseconds()` ` p := makePrimes(nOdd+1)` ` sieveDone := time.Nanoseconds()` ` aEven := p.doubleFactorial(nEven)` ` evenDone := time.Nanoseconds()` ` aOdd := p.doubleFactorial(nOdd)` ` oddDone := time.Nanoseconds()` ` fmt.Printf("sieve to %d -> %.2f sec\n",` ` nOdd+1, float64(sieveDone-start)/1e9)` ` fmt.Printf(" %d!! -> %.2f sec\n",` ` nEven, float64(evenDone-sieveDone)/1e9)` ` fmt.Printf(" %d!! -> %.2f sec\n",` ` nOdd, float64(oddDone-evenDone)/1e9)` ` ` ` dtrunc := int64(float64(aEven.BitLen())*.30103) - 10` ` var first, rest big.Int` ` rest.Exp(first.SetInt64(10), rest.SetInt64(dtrunc), nil)` ` first.Quo(aEven, &rest)` ` fstr := first.String()` ` fmt.Printf("%d!! begins %s... and has %d digits.\n",` ` nEven, fstr, int64(len(fstr))+dtrunc)` ` ` ` dtrunc = int64(float64(aOdd.BitLen())*.30103) - 10` ` rest.Exp(first.SetInt64(10), rest.SetInt64(dtrunc), nil)` ` first.Quo(aOdd, &rest)` ` fstr = first.String()` ` fmt.Printf("%d!! begins %s... and has %d digits.\n",` ` nOdd, fstr, int64(len(fstr))+dtrunc)` `}` ERATOSTHENES SCALA ERATOSTHENES C# ERATOSTHENES GO `// Prime number sieve, Eratosthenes (276-194 b.t. )` `// This implementation considers only multiples of primes` `// greater than 3, so the smallest value has to be mapped to 5.` `// Note: There is no multiplication operation in this function` `// and no call to a sqrt function.` `class Primes(val n: Int) {` ` ` ` private val primesOnBits =` ` Array[Long](1762821248, 848611808, 3299549660L, 2510511646L)` ` private val bitsPerInt: Int = 32` ` private val mask: Int = bitsPerInt - 1` ` private val log2Int: Int = 5` ` private var isComposite: Array[Long] = null` ` ` ` if (n < 386) {isComposite = primesOnBits; }` ` isComposite = new Array[Long]((n / (3 * bitsPerInt)) + 1: Int)` ` ` ` var d1 = 8; var d2 = 8; var p1 = 3; var p2 = 7; var s = 7` ` var s2 = 3; var l = 0; var c = 1; var max = n / 3; var inc = 0` ` var toggle: Boolean = false` ` ` ` while (s < max) // -- scan the sieve` ` {` ` // -- if a prime is found ...` ` if ((isComposite(l >> log2Int) & (1 << (l & mask))) == 0)` ` {` ` inc = p1 + p2 // -- ... cancel its multiples` ` ` ` // for (c <- s until max by inc) // slow!` ` c = s; while (c < max) // -- set c as composite` ` {` ` isComposite(c >> log2Int) |= 1 << (c & mask)` ` c += inc` ` }` ` ` ` c = s + s2; while (c < max) ` ` {` ` isComposite(c >> log2Int) |= 1 << (c & mask)` ` c += inc` ` }` ` }` ` ` ` l = l + 1` ` toggle = !toggle` ` if (toggle) {s += d2; d1 += 16; p1 += 2; p2 += 2; s2 = p2}` ` else {s += d1; d2 += 8; p1 += 2; p2 += 6; s2 = p1}` ` }` ` ` ` def iterator(min: Int, max: Int, visitor: (Int => Unit)): Unit =` ` {` ` // isComposite[0] ... isComposite[n] includes` ` // 5 <= prime numbers <= 96*(n+1) + 1` ` ` ` if (max > n) throw new IllegalArgumentException(` ` "Max not in sieve.")` ` val min2 = if (min < 2) 2 else min` ` ` ` if (min2 == 2) visitor(2)` ` if (min2 <= 3) visitor(3)` ` ` ` val absPos = (min2 + (min2 + 1) % 2) / 3 - 1` ` var index = absPos / bitsPerInt` ` var bitPos = absPos % bitsPerInt` ` var prime = 5+3*(bitsPerInt*index+bitPos)-(bitPos & 1)` ` var toggle: Boolean = (bitPos & 1) == 1` ` ` ` while (prime <= max)` ` {` ` var bitField = isComposite(index) >> bitPos` ` index += 1` ` while (bitPos < bitsPerInt)` ` {` ` if ((bitField & 1) == 0)` ` {` ` visitor(prime)` ` }` ` ` ` toggle = !toggle` ` prime = if (toggle) prime + 2 else prime + 4` ` if (prime > max) return` ` bitField >>= 1` ` bitPos += 1` ` }` ` bitPos = 0` ` }` ` }` `}` ` /// Prime number sieve, Eratosthenes (276-194 b.t. )` ` /// This implementation considers only multiples of primes` ` /// greater than 3, so the smallest value has to be mapped to 5.` ` /// Note: There is no multiplication operation in this function` ` /// and no call to a sqrt function.` ` public class Primes` ` {` ` const int bitsPerInt = 32;` ` const int mask = bitsPerInt - 1;` ` const int log2Int = 5;` ` private uint sieveLen = 0;` ` ` ` private static uint[] PrimesOnBits = {` ` 1762821248u, 848611808u, 3299549660u, 2510511646u };` ` ` ` private uint[] isComposite;` ` public delegate void Visitor(uint x);` ` ` ` public Primes(uint n)` ` {` ` sieveLen = n;` ` if (n < 386) { isComposite = PrimesOnBits; return; }` ` ` ` isComposite = new uint[(n / (3 * bitsPerInt)) + 1];` ` int d1 = 8, d2 = 8, p1 = 3, p2 = 7, s = 7, s2 = 3;` ` int l = 0, c = 1, max = (int)n / 3, inc;` ` bool toggle = false;` ` ` ` while (s < max) // -- scan the sieve` ` {` ` // -- if a prime is found ...` ` if ((isComposite[l >> log2Int] & (1u << (l++ & mask))) == 0)` ` {` ` inc = p1 + p2; // -- ... cancel its multiples` ` ` ` for (c = s; c < max; c += inc)` ` { // -- ... set c as composite` ` isComposite[c >> log2Int] |= 1u << (c & mask);` ` }` ` ` ` for (c = s + s2; c < max; c += inc)` ` {` ` isComposite[c >> log2Int] |= 1u << (c & mask);` ` }` ` }` ` ` ` if (toggle = !toggle) ` ` { s += d2; d1 += 16; p1 += 2; p2 += 2; s2 = p2; }` ` else { s += d1; d2 += 8; p1 += 2; p2 += 6; s2 = p1; }` ` }` ` }` ` ` ` public void Iterator(uint min, uint max, Visitor visitor)` ` {` ` // isComposite[0] ... isComposite[n] includes` ` // 5 <= primes numbers <= 96*(n+1) + 1` ` ` ` if (min < 2) min = 2;` ` if (max > sieveLen) throw ` ` new ArgumentOutOfRangeException("Max larger than sieve.");` ` ` ` if (min <= 2) visitor(2);` ` if (min <= 3) visitor(3);` ` ` ` int absPos = (int)((min + (min + 1) % 2) / 3 - 1);` ` int index = absPos / bitsPerInt;` ` int bitPos = absPos % bitsPerInt;` ` bool toggle = (bitPos & 1) == 1;` ` uint prime = (uint)(5+3*(bitsPerInt*index+bitPos)-(bitPos & 1));` ` ` ` while (prime <= max)` ` {` ` uint bitField = isComposite[index++] >> bitPos;` ` for (; bitPos < bitsPerInt; bitPos++)` ` {` ` if ((bitField & 1) == 0)` ` {` ` visitor(prime);` ` }` ` prime += (toggle = !toggle) ? 2u : 4u;` ` if (prime > max) return;` ` bitField >>= 1;` ` }` ` bitPos = 0;` ` }` ` }` `}` `// Prime number sieve, Eratosthenes (276-194 b.t. )` `// This implementation considers only multiples of primes` `// greater than 3, so the smallest value has to be mapped to 5.` `// word size dependent constants` `const ( ` ` bitsPerInt = 32` ` mask = bitsPerInt - 1` ` log2Int = 5` `)` `// holds completed sieve` `type primes struct {` ` sieveLen uint` ` isComposite []uint32` `}` `// constructor, completes sieve.` `func makePrimes(n uint) (ps *primes) {` ` ps = new(primes)` ` ps.sieveLen = n` ` if n < 386 {` ` ps.isComposite = []uint32{1762821248, 848611808, ` ` 3299549660, 2510511646}` ` return` ` }` ` ps.isComposite = make([]uint32, (n/(3*bitsPerInt))+1)` ` var (` ` d1, d2, p1, p2, s, s2 uint = 8, 8, 3, 7, 7, 3` ` l, c, max, inc uint = 0, 1, n / 3, 0` ` toggle bool` ` )` ` for s < max { // -- scan the sieve` ` // -- if a prime is found ...` ` if (ps.isComposite[l>>log2Int] & (1 << (l & mask))) == 0 {` ` inc = p1 + p2 // -- ... cancel its multiples` ` // -- ... set c as composite` ` for c = s; c < max; c += inc {` ` ps.isComposite[c>>log2Int] |= 1 << (c & mask)` ` }` ` for c = s + s2; c < max; c += inc {` ` ps.isComposite[c>>log2Int] |= 1 << (c & mask)` ` }` ` }` ` l++` ` toggle = !toggle` ` if toggle {` ` s += d2; d1 += 16; p1 += 2; p2 += 2; s2 = p2 ` ` } else {` ` s += d1; d2 += 8; p1 += 2; p2 += 6; s2 = p1 ` ` }` ` }` ` return` `}` `func (ps *primes) iterator(min, max uint, visitor func(uint)) bool {` ` // isComposite[0] ... isComposite[n] includes` ` // 5 <= primes numbers <= 96*(n+1) + 1` ` if min < 2 {` ` min = 2` ` }` ` if max > ps.sieveLen {` ` return false // Max larger than sieve` ` }` ` if min <= 2 {` ` visitor(2)` ` }` ` if min <= 3 {` ` visitor(3)` ` }` ` absPos := uint((min+(min+1)%2)/3 - 1)` ` index := absPos / bitsPerInt` ` bitPos := absPos % bitsPerInt` ` toggle := (bitPos & 1) == 1` ` prime := uint(5 + 3*(bitsPerInt*index+bitPos) - (bitPos & 1))` ` for prime <= max {` ` bitField := ps.isComposite[index] >> bitPos` ` index++` ` for ; bitPos < bitsPerInt; bitPos++ {` ` if (bitField & 1) == 0 {` ` visitor(prime)` ` }` ` toggle = !toggle` ` if toggle {` ` prime += 2` ` } else {` ` prime += 4` ` }` ` if prime > max {` ` return true` ` }` ` bitField >>= 1` ` }` ` bitPos = 0` ` }` ` return true` `}` ```var count = 0 def visitor(prime: Int): Unit = { count = count + 1 } def sieve(): Unit = { val n = 100000000 // Correct answer: 5,761,455 val start = System.nanoTime val p = new Primes(n) val stop = System.nanoTime println("Time taken to sieve = %s sec".format((stop-start)/1e9)) val start2 = System.nanoTime p.iterator(1, n, visitor) val stop2 = System.nanoTime println("Time taken to visit = %s sec".format((stop2-start2)/1e9)) println("Number of primes: " + count) } Time taken to sieve = 1.58 sec Time taken to visit = 0.44 sec Number of primes: 5761455 ``` ```static int count = 0; static void visitor(uint prime) { count = count + 1; } static void PrimeSieveTest() { var watch = new System.Diagnostics.Stopwatch(); var n = 100000000u; // Correct answer: 5,761,455 watch.Reset(); watch.Start(); var p = new Primes(n); watch.Stop(); Console.WriteLine("Time taken to sieve " + watch.Elapsed); watch.Reset(); watch.Start(); p.Iterator(1, n, visitor); watch.Stop(); Console.WriteLine("Time taken to visit " + watch.Elapsed); Console.WriteLine("Number of primes: " + count); } Time taken to sieve 00:00:01.10 Time taken to visit 00:00:00.49 Number of primes: 5761455 ``` `package main` `import (` ` "fmt"` ` "time"` `)` `const n = 100000000 // Correct answer: 5,761,455` `func main() {` ` var count uint` ` visitor := func(_ uint) {` ` count = count + 1` ` }` ` start := time.Nanoseconds()` ` p := makePrimes(n)` ` stop := time.Nanoseconds()` ` fmt.Printf("Time taken to sieve %.2f sec\n", ` ` float64(stop-start)/1e9)` ` start = time.Nanoseconds()` ` p.iterator(1, n, visitor)` ` stop = time.Nanoseconds()` ` fmt.Printf("Time taken to visit %.2f sec\n", ` ` float64(stop-start)/1e9)` ` fmt.Println("Number of primes: ", count)` `}` ` ` ` ` MathFun SCALA MathFun C# MathFun GO `object MathFun {` ` ` ` // Calibrate the treshhold` ` private val THRESHOLD_PRODUCT_SERIAL: Int = 28` ` ` ` def product(seq: Array[Long], start: Int, len: Int): BigInt =` ` {` ` if (len <= THRESHOLD_PRODUCT_SERIAL)` ` {` ` var sprod = BigInt(seq(start))` ` ` ` // for (i <- start + 1 until start + len) // slow!` ` val stop = start+len; var i = start+1; while(i < stop)` ` {` ` sprod *= BigInt(seq(i))` ` i += 1` ` }` ` ` ` return sprod` ` }` ` ` ` val halfLen = len >> 1` ` val lprod = product(seq, start, halfLen)` ` var rprod = product(seq, start + halfLen, len - halfLen)` ` return lprod * rprod` ` }` ` ` ` def bitCount(v: Int): Int =` ` {` ` var w = v - ((0xaaaaaaaa & v) >> 1)` ` w = (w & 0x33333333) + ((w >> 2) & 0x33333333)` ` w = w + (w >> 4) & 0x0f0f0f0f` ` w += w >> 8` ` w += w >> 16` ` ` ` return w & 0xff` ` }` ` ` ` def floorSqrt(n: Int): Int =` ` {` ` var a: Int = n` ` var b: Int = n` ` ` ` do {` ` a = b` ` b = (n / a + a) / 2` ` } while (b < a)` ` ` ` return a` ` }` `} ` ` public class MathFun` ` {` ` // Calibrate the treshholds` ` private const int THRESHOLD_PRODUCT_SERIAL = 24;` ` private const int THRESHOLD_PRODUCT_PARALLEL = 8000;` ` ` ` static public BigInt Product(long[] seq, int start, int len)` ` {` ` BigInt rprod, lprod;` ` int halfLen = len >> 1;` ` ` ` if (len > THRESHOLD_PRODUCT_PARALLEL)` ` {` ` rprod = BigInt.Zero; lprod = BigInt.Zero;` ` Parallel.Invoke(` ` () => { rprod = Product(seq, 0, halfLen); },` ` () => { lprod = Product(seq, halfLen, len - halfLen); }` ` );` ` return lprod * rprod;` ` }` ` ` ` if (len <= THRESHOLD_PRODUCT_SERIAL)` ` {` ` var sprod = new BigInt(seq[start]);` ` ` ` for (int i = start + 1; i < start + len; i++)` ` {` ` sprod *= seq[i];` ` }` ` return sprod;` ` }` ` ` ` lprod = Product(seq, start, halfLen);` ` rprod = Product(seq, start + halfLen, len - halfLen);` ` return lprod * rprod;` ` }` ` ` ` /// Bit count, loopfree!` ` public static uint BitCount(uint w)` ` {` ` w = w - ((0xaaaaaaaa & w) >> 1);` ` w = (w & 0x33333333) + ((w >> 2) & 0x33333333);` ` w = w + (w >> 4) & 0x0f0f0f0f;` ` w += w >> 8;` ` w += w >> 16;` ` ` ` return w & 0xff;` ` }` ` ` ` public static uint FloorSqrt(uint n)` ` {` ` uint a, b;` ` a = b = n;` ` do` ` {` ` a = b;` ` b = (n / a + a) / 2;` ` } while (b < a);` ` return a;` ` }` ` ` ` private MathFun() { }` ` }` `} ` `const productSerialThreshold = 24` ` ` `func product(seq []uint64) *big.Int {` ` if len(seq) <= productSerialThreshold {` ` var b big.Int` ` sprod := big.NewInt(int64(seq[0]))` ` for _, s := range seq[1:] {` ` b.SetInt64(int64(s))` ` sprod.Mul(sprod, &b)` ` }` ` return sprod` ` }` ` ` ` halfLen := len(seq) / 2` ` lprod := product(seq[0:halfLen])` ` rprod := product(seq[halfLen:])` ` return lprod.Mul(lprod, rprod)` `}` ` ` `func bitCount(w uint64) uint { // loopfree!` ` const (` ` ff = 1 << 64 - 1` ` mask1 = ff / 3` ` mask3 = ff / 5` ` maskf = ff / 17` ` maskp = maskf >> 3 & maskf` ` )` ` w -= w >> 1 & mask1` ` w = w & mask3 + w >> 2 & mask3` ` w = (w + w >> 4) & maskf` ` return uint(w * maskp >> 56)` `}` ` ` `func floorSqrt(n uint) uint {` ` for b := n; ; {` ` a := b` ` b = (n / a + a) / 2` ` if b >= a {` ` return a` ` }` ` }` ` return 0 // unreachable. required by current compiler.` `}` ~ ~ ~ `object Small {` ` val oddSwing = Array[Long](` ```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,9917826435L,583401555,20419054425L, 2268783825L,83945001525L,4418157975L,172308161025L, 34461632205L,1412926920405L,67282234305L,2893136075115L, 263012370465L,11835556670925L,514589420475L,24185702762325L, 8061900920775L,395033145117975L,15801325804719L, 805867616040669L,61989816618513L,3285460280781189L, 121683714103007L,6692604275665385L,956086325095055L, 54496920530418135L,1879204156221315L,110873045217057585L, 7391536347803839L,450883717216034179L,14544636039226909L, 916312070471295267L,916312070471295267L)``` ` ` ` val oddFactorial = Array[Long](` ```1,1,1,3,3,15,45,315,315,2835,14175,155925,467775, 6081075,42567525,638512875,638512875,10854718875L, 97692469875L,1856156927625L,9280784638125L,194896477400625L, 2143861251406875L,49308808782358125L,147926426347074375L, 3698160658676859375L)``` ` ` ` val oddDoubleFactorial = Array[Long](` ```1,1,1,3,1,15,3,105,3,945,15,10395,45,135135,315,2027025,315, 34459425,2835,654729075,14175,13749310575L,155925, 316234143225L,467775,7905853580625L,6081075,213458046676875L, 42567525,6190283353629375L,638512875,191898783962510625L, 638512875,6332659870762850625L,10854718875L)``` `}` ` /// Note: This is an attempt to optimize by accumulating` ` /// small products. Uses a trick I learned from Marco Bodrato.` `  ` ` class Factors ` ` {` ` private long[] factors;` ` private long maxProd, prod;` ` private int count;` ` ` ` public Factors(uint n)` ` {` ` int mem = (int)(0.63 * n / System.Math.Log(n));` ` factors = new long[mem];` ` }` ` ` ` public void Init()` ` {` ` maxProd = 1;` ` prod = 1;` ` count = 0;` ` }` ` ` ` public void SetMax(long max)` ` {` ` maxProd = long.MaxValue / max;` ` ` ` if (prod >= maxProd)` ` {` ` factors[count++] = prod;` ` prod = 1;` ` }` ` }` ` ` ` public void Add(long prime)` ` {` ` if (prod < maxProd)` ` { // the Bodrato trick:` ` prod *= prime; // accumulate small products first.` ` }` ` else` ` {` ` factors[count++] = prod;` ` prod = prime;` ` }` ` }` ` ` ` public BigInt Product()` ` {` ` factors[count++] = prod;` ` return MathFun.Product(factors, 0, count);` ` }` ` } ` `var smallOddSwing []int64 = []int64{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, 9917826435, 583401555, 20419054425,` ` 2268783825, 83945001525, 4418157975, 172308161025,` ` 34461632205, 1412926920405, 67282234305, 2893136075115,` ` 263012370465, 11835556670925, 514589420475, 24185702762325,` ` 8061900920775, 395033145117975, 15801325804719,` ` 805867616040669, 61989816618513, 3285460280781189,` ` 121683714103007, 6692604275665385, 956086325095055,` ` 54496920530418135, 1879204156221315, 110873045217057585,` ` 7391536347803839, 450883717216034179, 14544636039226909,` ` 916312070471295267, 916312070471295267}` ` ` `var smallOddFactorial []int64 = []int64{1, 1, 1, 3, 3,` ` 15, 45, 315, 315, 2835, 14175, 155925, 467775,` ` 6081075, 42567525, 638512875, 638512875, 10854718875,` ` 97692469875, 1856156927625, 9280784638125, 194896477400625,` ` 2143861251406875, 49308808782358125, 147926426347074375,` ` 3698160658676859375}` ` ` `var smallOddDoubleFactorial []int64 = []int64{1, 1, 1, 3, 1,` ` 15, 3, 105, 3, 945, 15, 10395, 45, 135135, 315, 2027025, 315,` ` 34459425, 2835, 654729075, 14175, 13749310575, 155925,` ` 316234143225, 467775, 7905853580625, 6081075, 213458046676875,` ` 42567525, 6190283353629375, 638512875, 191898783962510625,` ` 638512875, 6332659870762850625, 10854718875}` ` `