/// 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#)
/// Creative Commons Attribution-ShareAlike 3.0
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}
 

previous Back to the Homepage of Factorial Algorithms.