Primorial numbers
You are encouraged to solve this task according to the task description, using any language you may know.
Primorial numbers are those formed by multiplying successive prime numbers.
The primorial number series is:
- primorial(0) = 1 (by definition)
- primorial(1) = 2 (2)
- primorial(2) = 6 (2×3)
- primorial(3) = 30 (2×3×5)
- primorial(4) = 210 (2×3×5×7)
- primorial(5) = 2310 (2×3×5×7×11)
- primorial(6) = 30030 (2×3×5×7×11×13)
- ∙ ∙ ∙
To express this mathematically, primorialn is the product of the first n (successive) primes:
- ─── where is the kth prime number.
In some sense, generating primorial numbers is similar to factorials.
As with factorials, primorial numbers get large quickly.
- Task
- Show the first ten primorial numbers (0 ──► 9, inclusive).
- Show the length of primorial numbers whose index is: 10 100 1,000 10,000 and 100,000.
- Show the length of the one millionth primorial number (optional).
- Use exact integers, not approximations.
By length (above), it is meant the number of decimal digits in the numbers.
- Related tasks
- See also
11l
F get_primes(primes_count)
V limit = 17 * primes_count
V is_prime = [0B] * 2 [+] [1B] * (limit - 1)
L(n) 0 .< Int(limit ^ 0.5 + 1.5)
I is_prime[n]
L(i) (n * n .< limit + 1).step(n)
is_prime[i] = 0B
[Int] primes
L(prime) is_prime
I prime
primes.append(L.index)
I primes.len == primes_count
L.break
R primes
V primes = get_primes(100000)
F primorial(n)
BigInt r = 1
L(i) 0 .< n
r *= :primes[i]
R r
print(‘First ten primorials: ’(0.<10).map(n -> primorial(n)))
L(e) 6
V n = 10 ^ e
print(‘primorial(#.) has #. digits’.format(n, String(primorial(n)).len))- Output:
First ten primorials: [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870] primorial(1) has 1 digits primorial(10) has 10 digits primorial(100) has 220 digits primorial(1000) has 3393 digits primorial(10000) has 45337 digits primorial(100000) has 563921 digits
Arturo
primes: [2] ++ select.first: 99999 range.step: 2 3 ∞ => prime?
primorial: function [n][
if 0 = n -> return 1
product take primes n
]
print "First ten primorials:"
loop 0..9 => [print primorial &]
print ""
loop 1..5 'm -> print [
"primorial" 10^m "has" size ~"|primorial 10^m|" "digits"
]
- Output:
First ten primorials: 1 2 6 30 210 2310 30030 510510 9699690 223092870 primorial 10 has 10 digits primorial 100 has 220 digits primorial 1000 has 3393 digits primorial 10000 has 45337 digits primorial 100000 has 563921 digits
C
Uses a custom bit-sieve to generate primes, and GMP to keep track of the value of the primorial. Output takes ~3m to generate on a typical laptop.
#include <inttypes.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <gmp.h>
/* Eratosthenes bit-sieve */
int es_check(uint32_t *sieve, uint64_t n)
{
if ((n != 2 && !(n & 1)) || (n < 2))
return 0;
else
return !(sieve[n >> 6] & (1 << (n >> 1 & 31)));
}
uint32_t *es_sieve(const uint64_t nth, uint64_t *es_size)
{
*es_size = nth * log(nth) + nth * (log(log(nth)) - 0.9385f) + 1;
uint32_t *sieve = calloc((*es_size >> 6) + 1, sizeof(uint32_t));
for (uint64_t i = 3; i < sqrt(*es_size) + 1; i += 2)
if (!(sieve[i >> 6] & (1 << (i >> 1 & 31))))
for (uint64_t j = i * i; j < *es_size; j += (i << 1))
sieve[j >> 6] |= (1 << (j >> 1 & 31));
return sieve;
}
size_t mpz_number_of_digits(const mpz_t op)
{
char *opstr = mpz_get_str(NULL, 10, op);
const size_t oplen = strlen(opstr);
free(opstr);
return oplen;
}
#define PRIMORIAL_LIMIT 1000000
int main(void)
{
/* Construct a sieve of the first 1,000,000 primes */
uint64_t sieve_size;
uint32_t *sieve = es_sieve(PRIMORIAL_LIMIT, &sieve_size);
mpz_t primorial;
mpz_init_set_ui(primorial, 1);
uint64_t prime_count = 0;
int print = 1;
double unused;
for (uint64_t i = 2; i < sieve_size && prime_count <= PRIMORIAL_LIMIT; ++i) {
if (print) {
if (prime_count < 10)
gmp_printf("Primorial(%" PRIu64 ") = %Zd\n", prime_count, primorial);
/* Is the current number a power of 10? */
else if (!modf(log10(prime_count), &unused))
printf("Primorial(%" PRIu64 ") has %zu digits\n", prime_count, mpz_number_of_digits(primorial));
print = 0;
}
if (es_check(sieve, i)) {
mpz_mul_ui(primorial, primorial, i);
prime_count++;
print = 1;
}
}
free(sieve);
mpz_clear(primorial);
return 0;
}
- Output:
Primorial(0) = 1 Primorial(1) = 2 Primorial(2) = 6 Primorial(3) = 30 Primorial(4) = 210 Primorial(5) = 2310 Primorial(6) = 30030 Primorial(7) = 510510 Primorial(8) = 9699690 Primorial(9) = 223092870 Primorial(10) has 10 digits Primorial(100) has 220 digits Primorial(1000) has 3393 digits Primorial(10000) has 45337 digits Primorial(100000) has 563921 digits Primorial(1000000) has 6722809 digits
C#
The first 10 primorial numbers are calculated exactly, using integers. The number of digits are calculated exactly, using the summation of Log10()'s of the set of prime numbers. Nowhere in the task description does it say to calculate the larger-than-integer numbers, then count the digits, it just says to indicate the exact number of digits. This allows the program to execute rather quickly, as in just over 1/8 of a second at Tio.run. This is quite a bit quicker than the product tree algorithms presented elsewhere on this page that take 2 to 4 seconds on faster machines than are available at Tio.run.
using System;
class Program {
static int l;
static int[] gp(int n) {
var c = new bool[n]; var r = new int[(int)(1.28 * n)];
l = 0; r[l++] = 2; r[l++] = 3; int j, d, lim = (int)Math.Sqrt(n);
for (int i = 9; i < n; i += 6) c[i] = true;
for (j = 5, d = 4; j < lim; j += (d = 6 - d)) if (!c[j]) { r[l++] = j;
for (int k = j * j, ki = j << 1; k < n; k += ki) c[k] = true; }
for ( ; j < n; j += (d = 6 - d)) if (!c[j]) r[l++] = j; return r; }
static void Main(string[] args) {
var sw = System.Diagnostics.Stopwatch.StartNew();
var res = gp(15485864); sw.Stop();
double gpt = sw.Elapsed.TotalMilliseconds, tt;
var s = new string[19]; int si = 0;
s[si++] = String.Format("primes gen time: {0} ms", gpt); sw.Restart();
s[si++] = " Nth Primorial";
double y = 0; int x = 1, i = 0, lmt = 10;
s[si++] = String.Format("{0,7} {1}", 0, x);
while (true) {
if (i < 9) s[si++] = String.Format("{0,7} {1}", i + 1, x *= res[i]);
if (i == 8) s[si++] = " Nth Digits Time";
y += Math.Log10(res[i]);
if (++i == lmt) {
s[si++] = String.Format("{0,7} {1,-7} {2,7} ms", lmt, 1 + (int)y,
sw.Elapsed.TotalMilliseconds);
if ((lmt *= 10) > (int)1e6) break; } }
sw.Stop();
Console.WriteLine("{0}\n Tabulation: {1} ms", string.Join("\n", s),
tt = sw.Elapsed.TotalMilliseconds);
Console.Write(" Total:{0} ms", gpt + tt); } }
- Output:
primes gen time: 85.4715 ms
Nth Primorial
0 1
1 2
2 6
3 30
4 210
5 2310
6 30030
7 510510
8 9699690
9 223092870
Nth Digits Time
10 10 0.0695 ms
100 220 0.0898 ms
1000 3393 0.1335 ms
10000 45337 0.5501 ms
100000 563921 4.4646 ms
1000000 6722809 44.4457 ms
Tabulation: 44.4968 ms
Total:129.9683 ms
C++
#include <gmpxx.h>
#include <primesieve.hpp>
#include <cstdint>
#include <iomanip>
#include <iostream>
size_t digits(const mpz_class& n) { return n.get_str().length(); }
mpz_class primorial(unsigned int n) {
mpz_class p;
mpz_primorial_ui(p.get_mpz_t(), n);
return p;
}
int main() {
uint64_t index = 0;
primesieve::iterator pi;
std::cout << "First 10 primorial numbers:\n";
for (mpz_class pn = 1; index < 10; ++index) {
unsigned int prime = pi.next_prime();
std::cout << index << ": " << pn << '\n';
pn *= prime;
}
std::cout << "\nLength of primorial number whose index is:\n";
for (uint64_t power = 10; power <= 1000000; power *= 10) {
uint64_t prime = primesieve::nth_prime(power);
std::cout << std::setw(7) << power << ": "
<< digits(primorial(prime)) << '\n';
}
return 0;
}
- Output:
First 10 primorial numbers:
0: 1
1: 2
2: 6
3: 30
4: 210
5: 2310
6: 30030
7: 510510
8: 9699690
9: 223092870
Length of primorial number whose index is:
10: 10
100: 220
1000: 3393
10000: 45337
100000: 563921
1000000: 6722809
Clojure
Single Process
Using HashMap prime number generation from https://rosettacode.org/wiki/Sieve_of_Eratosthenes#Unbounded_Versions
(ns example
(:gen-class))
; Generate Prime Numbers (Implementation from RosettaCode--link above)
(defn primes-hashmap
"Infinite sequence of primes using an incremental Sieve or Eratosthenes with a Hashmap"
[]
(letfn [(nxtoddprm [c q bsprms cmpsts]
(if (>= c q) ;; only ever equal
; Update cmpsts with primes up to sqrt c
(let [p2 (* (first bsprms) 2),
nbps (next bsprms),
nbp (first nbps)]
(recur (+ c 2) (* nbp nbp) nbps (assoc cmpsts (+ q p2) p2)))
(if (contains? cmpsts c)
; Not prime
(recur (+ c 2) q bsprms
(let [adv (cmpsts c), ncmps (dissoc cmpsts c)]
(assoc ncmps
(loop [try (+ c adv)] ;; ensure map entry is unique
(if (contains? ncmps try)
(recur (+ try adv))
try))
adv)))
; prime
(cons c (lazy-seq (nxtoddprm (+ c 2) q bsprms cmpsts))))))]
(do (def baseoddprms (cons 3 (lazy-seq (nxtoddprm 5 9 baseoddprms {}))))
(cons 2 (lazy-seq (nxtoddprm 3 9 baseoddprms {}))))))
;; Generate Primorial Numbers
(defn primorial [n]
" Function produces the nth primorial number"
(if (= n 0)
1 ; by definition
(reduce *' (take n (primes-hashmap))))) ; multiply first n primes (retrieving primes from lazy-seq which generates primes as needed)
;; Show Results
(let [start (System/nanoTime)
elapsed-secs (fn [] (/ (- (System/nanoTime) start) 1e9))] ; System start time
(doseq [i (concat (range 10) [1e2 1e3 1e4 1e5 1e6])
:let [p (primorial i)]] ; Generate ith primorial number
(if (< i 10)
(println (format "primorial ( %7d ) = %10d" i (biginteger p))) ; Output for first 10
(println (format "primorial ( %7d ) has %8d digits\tafter %.3f secs" ; Output with time since starting for remainder
(long i) (count (str p)) (elapsed-secs))))))
- Output:
primorial ( 0 ) = 1 primorial ( 1 ) = 2 primorial ( 2 ) = 6 primorial ( 3 ) = 30 primorial ( 4 ) = 210 primorial ( 5 ) = 2310 primorial ( 6 ) = 30030 primorial ( 7 ) = 510510 primorial ( 8 ) = 9699690 primorial ( 9 ) = 223092870 primorial ( 100 ) has 220 digits after 0.012 secs primorial ( 1000 ) has 3393 digits after 0.048 secs primorial ( 10000 ) has 45337 digits after 0.284 secs primorial ( 100000 ) has 563921 digits after 7.731 secs primorial ( 1000000 ) has 6722809 digits after 706.593 secs Using: i7 920 @ 2.67 GHz CPU with Windows 10 /64 bit OS
Parallel Process
Using HashMap prime number generation from https://rosettacode.org/wiki/Sieve_of_Eratosthenes#Unbounded_Versions
(ns example
(:gen-class))
(defn primes-hashmap
"Infinite sequence of primes using an incremental Sieve or Eratosthenes with a Hashmap"
[]
(letfn [(nxtoddprm [c q bsprms cmpsts]
(if (>= c q) ;; only ever equal
; Update cmpsts with primes up to sqrt c
(let [p2 (* (first bsprms) 2),
nbps (next bsprms),
nbp (first nbps)]
(recur (+ c 2) (* nbp nbp) nbps (assoc cmpsts (+ q p2) p2)))
(if (contains? cmpsts c)
; Not prime
(recur (+ c 2) q bsprms
(let [adv (cmpsts c), ncmps (dissoc cmpsts c)]
(assoc ncmps
(loop [try (+ c adv)] ;; ensure map entry is unique
(if (contains? ncmps try)
(recur (+ try adv))
try))
adv)))
; prime
(cons c (lazy-seq (nxtoddprm (+ c 2) q bsprms cmpsts))))))]
(do (def baseoddprms (cons 3 (lazy-seq (nxtoddprm 5 9 baseoddprms {}))))
(cons 2 (lazy-seq (nxtoddprm 3 9 baseoddprms {}))))))
;; Number of workers (threads) based upon number of available processors
(def workers
(+ 2 (.. Runtime getRuntime availableProcessors)))
;; Generate of primorial numbers (using multiple processors)
(defn primorial [n]
(if (= n 0)
1
;(reduce mul+ (pmap #(reduce *' %) (partition-all (max workers (long (/ n workers))) (take n (primes-hashmap)))))));(*' allows for big integer arithmetic as needed)
(->> ; Threads (i.e. pipes) sequence of expressions
(take n (primes-hashmap) ) ; generate primes
(partition-all (max workers (long (/ n workers)))) ; partition primes amongst workers
(pmap #(reduce *' %)) ; Multiply primes in each worker in parallel
(reduce *')))) ; multiply results of all workers together
;; Generate and Time Output
(let [start (System/nanoTime)
elapsed-secs (fn [] (/ (- (System/nanoTime) start) 1e9))] ; System start time
(doseq [i (concat (range 10) [1e2 1e3 1e4 1e5 1e6])
:let [p (primorial i)]] ; Generate ith primorial number
(if (< i 10)
(println (format "primorial ( %7d ) = %10d" i (biginteger p))) ; Output for first 10
(println (format "primorial ( %7d ) has %8d digits\tafter %.3f secs" ; Output with time since starting for remainder
(long i) (count (str p)) (elapsed-secs))))))
- Output:
primorial ( 0 ) = 1 primorial ( 1 ) = 2 primorial ( 2 ) = 6 primorial ( 3 ) = 30 primorial ( 4 ) = 210 primorial ( 5 ) = 2310 primorial ( 6 ) = 30030 primorial ( 7 ) = 510510 primorial ( 8 ) = 9699690 primorial ( 9 ) = 223092870 primorial ( 100 ) has 220 digits after 0.016 secs primorial ( 1000 ) has 3393 digits after 0.050 secs primorial ( 10000 ) has 45337 digits after 0.364 secs primorial ( 100000 ) has 563921 digits after 2.619 secs primorial ( 1000000 ) has 6722809 digits after 69.812 secs Using: i7 920 @ 2.67 GHz CPU with Windows 10 /64 bit OS
CLU
% This program uses the 'bigint' cluster from
% the 'misc.lib' included with PCLU.
isqrt = proc (s: int) returns (int)
x0: int := s/2
if x0=0 then return(s) end
x1: int := (x0 + s/x0)/2
while x1 < x0 do
x0 := x1
x1 := (x0 + s/x0)/2
end
return(x0)
end isqrt
sieve = proc (n: int) returns (array[bool])
prime: array[bool] := array[bool]$fill(0,n+1,true)
prime[0] := false
prime[1] := false
for p: int in int$from_to(2, isqrt(n)) do
if prime[p] then
for c: int in int$from_to_by(p*p,n,p) do
prime[c] := false
end
end
end
return(prime)
end sieve
% Calculate the N'th primorial given a boolean array denoting primes
primorial = proc (n: int, prime: array[bool])
returns (bigint) signals (out_of_primes)
% 0'th primorial = 1
p: bigint := bigint$i2bi(1)
for i: int in array[bool]$indexes(prime) do
if ~prime[i] then continue end
if n=0 then break end
p := p * bigint$i2bi(i)
n := n-1
end
if n>0 then signal out_of_primes end
return(p)
end primorial
% Find the length in digits of a bigint without converting it to a string.
% The naive way takes over an hour to count the digits for p(100000),
% this one ~5 minutes.
n_digits = proc (n: bigint) returns (int)
own zero: bigint := bigint$i2bi(0)
own ten: bigint := bigint$i2bi(10)
digs: int := 1
dstep: int := 1
tenfac: bigint := ten
step: bigint := ten
while n >= tenfac do
digs := digs + dstep
step := step * ten
dstep := dstep + 1
next: bigint := tenfac*step
if n >= next then
tenfac := next
else
step, dstep := ten, 1
tenfac := tenfac*step
end
end
return(digs)
end n_digits
start_up = proc ()
po: stream := stream$primary_output()
% Sieve a million primes
prime: array[bool] := sieve(15485863)
% Show the first 10 primorials
for i: int in int$from_to(0,9) do
stream$puts(po, "primorial(" || int$unparse(i) || ") = ")
stream$putright(po, bigint$unparse(primorial(i, prime)), 15)
stream$putl(po, "")
end
% Show the length of some bigger primorial numbers
for tpow: int in int$from_to(1,5) do
p_ix: int := 10**tpow
stream$puts(po, "length of primorial(")
stream$putright(po, int$unparse(p_ix), 7)
stream$puts(po, ") = ")
stream$putright(po, int$unparse(n_digits(primorial(p_ix, prime))), 7)
stream$putl(po, "")
end
end start_up- Output:
primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 length of primorial( 10) = 10 length of primorial( 100) = 220 length of primorial( 1000) = 3393 length of primorial( 10000) = 45337 length of primorial( 100000) = 563921
Common Lisp
(defun primorial-number-length (n w)
(values (primorial-number n) (primorial-length w)))
(defun primorial-number (n)
(loop for a below n collect (primorial a)))
(defun primorial-length (w)
(loop for a in w collect (length (write-to-string (primorial a)))))
(defun primorial (n &optional (m 1) (k -1) (z 1) &aux (f (primep m)))
(if (= k n) z (primorial n (1+ m) (+ k (if f 1 0)) (if f (* m z) z))))
(defun primep (n)
(loop for a from 2 to (isqrt n) never (zerop (mod n a))))
- Output:
> (primorial-number-length 10 '(100 1000 10000 100000)) (1 2 6 30 210 2310 30030 510510 9699690 223092870) (220 3393 45337 563921)
D
import std.stdio;
import std.format;
import std.bigint;
import std.math;
import std.algorithm;
int sieveLimit = 1300_000;
bool[] notPrime;
void main()
{
// initialize
sieve(sieveLimit);
// output 1
foreach (i; 0..10)
writefln("primorial(%d): %d", i, primorial(i));
// output 2
foreach (i; 1..6)
writefln("primorial(10^%d) has length %d", i, count(format("%d", primorial(pow(10, i)))));
}
BigInt primorial(int n)
{
if (n == 0) return BigInt(1);
BigInt result = BigInt(1);
for (int i = 0; i < sieveLimit && n > 0; i++)
{
if (notPrime[i]) continue;
result *= BigInt(i);
n--;
}
return result;
}
void sieve(int limit)
{
notPrime = new bool[limit];
notPrime[0] = notPrime[1] = true;
auto max = sqrt(cast (float) limit);
for (int n = 2; n <= max; n++)
{
if (!notPrime[n])
{
for (int k = n * n; k < limit; k += n)
{
notPrime[k] = true;
}
}
}
}
- Output:
primorial(0): 1 primorial(1): 2 primorial(2): 6 primorial(3): 30 primorial(4): 210 primorial(5): 2310 primorial(6): 30030 primorial(7): 510510 primorial(8): 9699690 primorial(9): 223092870 primorial(10^1) has length 10 primorial(10^2) has length 220 primorial(10^3) has length 3393 primorial(10^4) has length 45337 primorial(10^5) has length 563921
Delphi
See Pascal.
EasyLang
fastfunc isprim num .
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
fastfunc nextprim prim .
repeat
prim += 1
until isprim prim = 1
.
return prim
.
func[] big n .
while n > 0
r[] &= n mod 10000000
n = n div 10000000
.
return r[]
.
func[] bnmul a[] b[] .
# this multiplication is limited to 300 digits accuracy
len r[] len a[] + len b[]
max = higher 1 (len a[] - 45)
for ia = max to len a[]
h = 0
for ib = 1 to len b[]
h += r[ia + ib - 1] + b[ib] * a[ia]
r[ia + ib - 1] = h mod 10000000
h = h div 10000000
.
r[ia + ib - 1] += h
.
while r[$] = 0
len r[] -1
.
return r[]
.
func$ str bn[] .
for i = len bn[] downto 1
s$ &= bn[i]
.
return s$
.
func bndlen b[] .
return (len b[] - 1) * 7 + floor log10 b[$] + 1
.
primor[] = [ 1 ]
prim = 2
i = 0
while i < 10
print str primor[]
primor[] = bnmul primor[] big prim
prim = nextprim prim
i += 1
.
print ""
print bndlen primor[]
for lim in [ 100 1000 10000 100000 ]
while i < lim
primor[] = bnmul primor[] big prim
prim = nextprim prim
i += 1
.
print bndlen primor[]
.
- Output:
1 2 6 30 210 2310 30030 510510 9699690 223092870 10 220 3393 45337 563921
Elixir
Prime generator works too inefficiently to generate 1 million in a short time, but it will work eventually. This solution is efficient up to 100,000 primes.
defmodule SieveofEratosthenes do
def init(lim) do
find_primes(2,lim,(2..lim))
end
def find_primes(count,lim,nums) when (count * count) > lim do
nums
end
def find_primes(count,lim,nums) when (count * count) <= lim do
find_primes(count+1,lim,Enum.reject(nums,&(rem(&1,count) == 0 and &1 > count)))
end
end
defmodule Primorial do
def first(n,primes) do
s = 0..9 |> Stream.map(fn n -> Enum.at(primes,n) end)
(0..n-1)
|> Enum.map(fn a -> s
|> Enum.take(a)
|> Enum.reduce(1, fn b,c -> b*c end)
|> format(a) end)
end
def numbers(lims,primes) do
numbers(lims,primes,[])
end
def numbers([],_primes,vals) do
vals
|> Enum.reverse
|> Enum.map(fn {m,n} -> str_fr(n,m) end)
end
def numbers([lim|lims],primes,vals) do
numbers(lims,primes,[{lim,number_length(primes,lim)}] ++ vals)
end
defp number_length(primes,n) do
primes
|> Enum.take(n)
|> Enum.reduce(fn a,b -> a * b end)
|> Integer.to_string
|> String.length
end
defp format(pri,i), do: IO.puts("Primorial #{i}: #{pri}")
defp str_fr(pri,i), do: IO.puts("Primorial #{i} has length: #{pri}")
end
Primorial.first(10,SieveofEratosthenes.init(50))
Primorial.numbers([10,100,1_000,10_000,100_000],SieveofEratosthenes.init(1_300_000))
- Output:
Primorial 0: 1
Primorial 1: 2
Primorial 2: 6
Primorial 3: 30
Primorial 4: 210
Primorial 5: 2310
Primorial 6: 30030
Primorial 7: 510510
Primorial 8: 9699690
Primorial 9: 223092870
Primorial 10 has length: 10
Primorial 100 has length: 220
Primorial 1000 has length: 3393
Primorial 10000 has length: 45337
Primorial 100000 has length: 563921
F#
This task uses Extensible Prime Generator (F#)
// Primorial Numbers. Nigel Galloway: August 3rd., 2021
primes32()|>Seq.scan((*)) 1|>Seq.take 10|>Seq.iter(printf "%d "); printfn "\n"
[10;100;1000;10000;100000]|>List.iter(fun n->printfn "%d" ((int)(System.Numerics.BigInteger.Log10 (Seq.item n (primesI()|>Seq.scan((*)) 1I)))+1))
- Output:
1 2 6 30 210 2310 30030 510510 9699690 223092870 10 220 3393 45337 563921
Factor
USING: formatting kernel literals math math.functions
math.primes sequences ;
IN: rosetta-code.primorial-numbers
CONSTANT: primes $[ 1,000,000 nprimes ]
: digit-count ( n -- count ) log10 floor >integer 1 + ;
: primorial ( n -- m ) primes swap head product ;
: .primorial ( n -- ) dup primorial "Primorial(%d) = %d\n"
printf ;
: .digit-count ( n -- ) dup primorial digit-count
"Primorial(%d) has %d digits\n" printf ;
: part1 ( -- ) 10 iota [ .primorial ] each ;
: part2 ( -- ) { 10 100 1000 10000 100000 1000000 }
[ .digit-count ] each ;
: main ( -- ) part1 part2 ;
MAIN: main
- Output:
Primorial(0) = 1 Primorial(1) = 2 Primorial(2) = 6 Primorial(3) = 30 Primorial(4) = 210 Primorial(5) = 2310 Primorial(6) = 30030 Primorial(7) = 510510 Primorial(8) = 9699690 Primorial(9) = 223092870 Primorial(10) has 10 digits Primorial(100) has 220 digits Primorial(1000) has 3393 digits Primorial(10000) has 45337 digits Primorial(100000) has 563921 digits Primorial(1000000) has 6722809 digits
Fortran
The Plan
Since exact integer arithmetic is required, there can be no use of various approximations and the primorial numbers very soon exceed the integer limits, even faster than do factorials. As the one who provided the example of calculating factorials in the wikipædia article for multi-precision arithmetic, I am hoist by my own petard.
Prepare a list of prime numbers
The first task is to prepare a collection of the first million prime numbers. I already have disc files for up to 4,294,984,663 but this would not be portable so instead a subroutine to prepare the values, which is possibly faster than reading from a disc file anyway. Since this is a largeish set, the division method is probably not as good as the sieve method. Both methods could be employed to generate consecutive prime numbers, interlaced with their use to calculate the next primorial number, but divide-and-conquer is more attractive, so, first, prepare the primes. The sieve is arranged to represent only odd numbers so that there is no tedious pass with every other element, and further, by arranging the length of the sieve to be a primorial number, the pattern of the knock-outs for the first few primes is the same for each usage. A span of 30030 (spanned by an array of 15015 elements) seems reasonable, as the next size up would be 510510 which is a bit big. It would be nice if the initial-value specification protocol could generate that pattern via the repetition of the components in an expression that would be evaluated by the compiler - something like (3-pattern) & (5-pattern) & (7-pattern) & (11-pattern) & (13-pattern), with each pattern being generated by a suitable repetition count - I am certainly not going to write out a sequence of 15015 constants, and there are limits on a statement length anyway so devising a programme to generate the required source statements is not promising. But alas, I can't see how to do this either via a DATA statement or a PARAMETER specification. So, array START is initialised by executable statements that sieve out multiples of 3, 5, 7, 11, and 13 - including their own appearances. Then, a proper sieve process is started, that augments the PRIME collection as it goes at the beginning since the first sieve pass needs prime numbers that are not initially in array PRIME. On the other hand, by abandoning the special initial value array one could double the size of the sieve using no more storage, which might be better. When the sieve process is working with large numbers, there will be fewer and fewer "hits" for each sieve number. With a sieve span of 30,030, a prime of the order of 300,000 will have only one chance in ten of falling within a particular surge (and if so delivers only one hit within it), and the effort for the misses will have been wasted. For this task, the largest prime required for sieving is less than 4,000 so each prime is still making a few hits.
Because the sieve does not represent even numbers, the determination of a starting point is a little tricky. For a given prime P, the requirement is to find the first odd multiple of P beyond N0 (the start of a span, an even number), and if this value is beyond the end of the span (marked by NN) then there is no need to step through the span with P. Conflating this test with the test for the start point of P*P being beyond the end of the span is a mistake that, with a span of 30, led to 361 being declared a prime, because 360/17 = 21+ so 22 would be the start multiple for 17 except that 22 is even so 23 is used, and 23*17 = 391 which is beyond the end of the span so, wrongly, the sieving was thereby deemed complete and 19 was not tried, a mistake because 19*19 = 361.
The sieve array is declared LOGICAL*1 as it is painful to use the default of a 32-bit storage word to represent a boolean variable. Alas, accessing a single byte of storage may take more time, especially when placing a value. Still worse would be the unpacking of say a 32-bit integer to get at individual bits in isolation. With large arrays, a lot of storage could be saved, but execution speed would be poor unless the hardware offered support for single-bit access.
And after all that fuss, the production of the first million primes is completed in a blink.
Multi-precision multiplication
The multiplication of the multi-precision number is done on the cheap by passing the prime as an ordinary integer, not as a big number. The millionth prime is 15,485,863 which fits easily into a 32-bit integer that can hold values up to 2,147,483,647 but this implies a limit on the base of the arithmetic of the big number scheme. If it works in base 1,000 then a digit multiplied by that large prime would yield up to 15,485,863,000 - which exceeds the 32-bit integer limit, while if it works in base 10, the largest number would be up to 154,858,630, well within range. But at the cost of many more digits in the number, and correspondingly longer computation time. Some experiments with calculating Primorial(100000):
Base: 10 100 1,000 10,000 100,000 Secs: 554 278 185 117 52 - but wrong! 64-bit 300 241
The overflow happens in D*N + C, and converting D from INTEGER*4 to INTEGER*8 prevented the overflow, but introduces its own slowness. Activating "maximum optimisations" reduced 241 to 199, and reverting to base 100 with no INTEGER*8 converted 278 to 133. Then, abandoning the array bound checking reduced it further, to 121. The end result of the trick is to limit the base to 100. Given that, INTEGER*2 could be used for DIGIT, but a trial run showed almost the same time taken for Primorial(100000).
Hopefully, the compiler recognises the connection between the consecutive statements D = B.DIGIT(I) then D = D*N + C (rather than D = B.DIGIT(I)*N + C) which was done to facilitate trying INTEGER*8 D without needing to blabber on about conversion routines. But more importantly, one can dream that the compiler will recognise the classic sequence
B.DIGIT(I) = MOD(D,BIGBASE)
C = D/BIGBASE
and avoid performing two divisions. When an integer division is performed, the quotient appears in one register and the remainder in another and when writing in assembler with access to such registers there is no need for the two divisions as is forced by the use of high-level languages. An alternative method, once consecutive primorials are not required, would be to compute the big number product of say ten consecutive primes then multiply the main big number by that, using multi-precision arithmetic for both. Abandoning the trick that restricts the base to 100 (for the current upper limit of the millionth prime) would allow the use of larger bases. Another possibility would be to calculate in a base more directly matching that of the computer (since the variable-length decimal arithmetic of the IBM1620 et al is no longer available), for instance base 65536 with sixteen-bit words, etc. Inspection of the result to calculate the number of decimal digits required would not be troublesome. In such a case, one might try something like
INTEGER*4 D !A 32-bit product.
INTEGER*2 II(2),C,R !Some 16-bit variables.
EQUIVALENCE (D,II) !Align.
EQUIVALENCE (II(1),C),(II(2),R) !Carry in the high order half of D, Result in the low.
Supposing that unsigned 16-bit arithmetic was available then the product in D need not be split into a carry and a digit via the MOD and divide as above. But this is unlikely. Only in assembly would good code be possible.
Source
The source style takes advantage of F90 for the MODULE protocol that saves the annoyance of otherwise having to prepare a collection of COMMON statements for those who would not write a single mainline for the job. With F90 the lower bound of array BIT need no longer be one, otherwise, with older Fortran the expressions would have to be recast and one added or subtracted as appropriate. This is a simple (but error-prone) clerical task: computers excel at clerical tasks, so let the computer do it.
Also helpful is the TYPE statement to define a big number with the use of somewhat explanatory nomenclature. Associated names (such as BIGORDER, BIGBASE, etc.) could also have been defined as a part of a TYPE definition, but the extra blabber didn't seem worthwhile when one set only would be used. So, for them the old-style usage of names with a structure in the name's form. More serious is the assistance offered by the abilities of the PARAMETER statement to define constants with interrelationships, here especially with the use of BIGORDER.
Various F90 array statements provide a convenience that could otherwise be achieved via explicit DO-loops, but rather more difficult to do without are the I0 format code (which should be just "I") that reveals an integer value with a size to fit the number, and the I<ORDER>.<ORDER> usage of two features: the parameterisation of what would otherwise be constants, and the F90 augmentation of Iw to Iw.d whereby leading zero digits are supplied. This enables the value of the big number to be revealed as a proper string of digits. Before F90, this would have to be done via writing to a text variable and messing about. Thus, the topmost digit of a big number, B.DIGIT(B.LAST), is sent out with I0 format so that there are neither leading spaces nor zeroes, then the subsequent digits are rolled with I<ORDER>.<ORDER> so that any necessary leading zeroes for each of them are provided - else there would be gaps within the digit string as presented should the base be greater than ten and a B.DIGIT value be less than ten. Then there's the usual trick that once a FORMAT sequence is exhausted by the elements of an I/O list, a new line is started and the scan of the codes (of FORMAT 101) resumes at the rightmost open bracket. This is not exercised for the smaller numbers that are shown in full, but was during testing.
On the other hand, the modifying prefix nP for E (and F) format codes has long been available. This shifts the position of the decimal point left or right by n places, and with the E format code modifies the exponent accordingly. Thus, 1.2345E+5 can become .12345E+6 via the prefix -1P as in FORMAT 113. Evidently, this is the same value. When used with a F code there is no exponent part to modify accordingly, but here, the exponent is calculated separately (and presented via the I0 format code), and, one is added in the WRITE statement's output expressions, thus I4 + 1 and I8 + 1. The point of all this is that the resulting exponent part gives the number of decimal digits in the number, as per the specification.
MODULE BIGNUMBERS !Limited services: decimal integers, no negative numbers.
INTEGER BIGORDER !A limit attempt at generality.
PARAMETER (BIGORDER = 2) !This is the order of the base of the big number arithmetic.
INTEGER BIGBASE,BIGLIMIT !Sized thusly.
PARAMETER (BIGBASE = 10**BIGORDER, BIGLIMIT = 8888888/BIGORDER) !Enough?
TYPE BIGNUM !So, a big number is simple.
INTEGER LAST !This many digits (of size BIGBASE) are in use.
INTEGER DIGIT(BIGLIMIT) !The digits, in ascending power order.
END TYPE BIGNUM !So much for that.
CONTAINS !Now for some assistants.
SUBROUTINE BIGMULT(B,N) !B:=B*N; Multiply by an integer possibly bigger than the base.
TYPE(BIGNUM) B !The worker.
INTEGER N !A computer number, not a multi-digit number.
INTEGER D !Must be able to hold (BIGBASE - 1)*N + C
INTEGER C !The carry to the next digit.
INTEGER I !A stepper.
C = 0 !No previous digit to carry from.
DO I = 1,B.LAST !Step through the digits, upwards powers.
D = B.DIGIT(I) !Grab a digit.
D = D*N + C !Apply the multiply.
B.DIGIT(I) = MOD(D,BIGBASE) !Place the resulting digit.
C = D/BIGBASE !Agony! TWO divisions per step!!
END DO !On to the next digit up.
DO WHILE(C .GT. 0) !Now spread the last carry to further digits.
B.LAST = B.LAST + 1 !Up one more.
IF (B.LAST .GT. BIGLIMIT) STOP "Overflow by multiply!" !Perhaps not.
B.DIGIT(B.LAST) = MOD(C,BIGBASE) !The digit.
C = C/BIGBASE !The carry may be large, if N is large.
END DO !So slog on until it is gone.
END SUBROUTINE BIGMULT !Primary school stuff.
END MODULE BIGNUMBERS !No fancy tricks.
MODULE ERATOSTHENES !Prepare an array of prime numbers.
Considers odd numbers only as the pattern is very simple. Some trickery as a consequence.
INTEGER NP,LASTP !Counters.
PARAMETER (LASTP = 1000000) !The specified need.
INTEGER PRIME(0:LASTP),PREZAP !Initialisation is rather messy.
PARAMETER (PREZAP = 6) !Up to PRIME(6) = 13.
DATA NP/PREZAP/, PRIME(0:PREZAP)/1,2,3,5,7,11,13/ !Not counting the "zeroth" prime, 1.
CONTAINS !Some tricky stuff/
SUBROUTINE PREPARE PRIMES !Fetch a limited copy of the Platonic ideal.
INTEGER SURGE,LB !A sieve has a certain rather special size.
PARAMETER (SURGE = 30030, LB = SURGE/2 - 1) != 2*3*5*7*11*13.
LOGICAL*1 BIT(0:LB),START(0:LB)!Two such arrays, thanks.
INTEGER N0,NN !Bounds for the current sieve span.
INTEGER I,P,IP !Assistants.
C The scheme for a cycle of 2*3*5 = 30, remembering that even numbers are not involved so BIT(0:14).
C | surge 1 | surge 2 | surge 3 |
C N = | 1 1 1 1 1 2 2 2 2 2|3 3 3 3 3 4 4 4 4 4 5 5 5 5 5|6 6 6 6 6 7 7 7 7 7 8 8 8 8 8|9 9 9 9 9...
C |1 3 5 7 9 1 3 5 7 9 1 3 5 7 9|1 3 5 7 9 1 3 5 7 9 1 3 5 7 9|1 3 5 7 9 1 3 5 7 9 1 3 5 7 9|1 3 5 7 9...
C BIT(index) | 1 1 1 1 1| 1 1 1 1 1| 1 1 1 1 1|
C |0 1 2 3 4 5 6 7 8 9 0 1 2 3 4|0 1 2 3 4 5 6 7 8 9 0 1 2 3 4|0 1 2 3 4 5 6 7 8 9 0 1 2 3 4|0 1 2 3 4...
c 3 step | * * * * * | * * * * * | * * * * * | * *
c 5 step | * * * | * * * | * * * | *
c 7 step | x x | x * | * * |*
Concoct the initial state, once only, that repeats every SURGE.
START = .TRUE. !Prepare the field.
DO I = 2,PREZAP !Only odd numbers are represented, so no PRIME(1) = 2..
P = PRIME(I) !Select a step.
START(P/2:LB:P) = .FALSE. !Knock out multiples of P.
END DO !This pattern is palindromic.
NN = 0 !Syncopation. Where the previous surge ended.
Commence a pass through the BIT sieve.
10 N0 = NN !BIT(0) corresponds to N0 + 1, BIT(i) to N0 + 1 + 2i.
NN = NN + SURGE !BIT(LB) to NN - 1.
BIT = START !Pre-zapped for lesser primes.
IP = PREZAP !Syncopation. The last pre-zapped prime.
11 IP = IP + 1 !The next prime to sieve with.
IF (IP.GT.NP) CALL SCANFOR(.TRUE.) !Whoops, not yet to hand!
P = PRIME(IP) !Now grab it.
IF (P*P.GE.NN) GO TO 12 !If P*P exceeds the end, so will larger P.
I = N0/P + 1 !First multiple of P past N0.
IF (I.LT.P) I = P !Less than P is superfluous: the position was zapped by earlier action.
IF (MOD(I,2).EQ.0) I = I + 1 !If even, advance to the next odd multiple. Such as P.
I = I*P !The first number to zap. It will always be odd.
IF (I.LT.NN) THEN !Within the span?
I = (I - N0 - 1)/2 !Yes. Its offset into the current span.
BIT(I:LB:P) = .FALSE. !Zap every P'th position.
END IF !So much for that P.
GO TO 11 !On to the next.
Completed the passes. Scan for survivors.
12 CALL SCANFOR(.FALSE.) !All, not just the first new prime.
IF (NP.LT.LASTP) GO TO 10 !Another batch?
RETURN !Done.
CONTAINS !Fold two usages into one routine.
SUBROUTINE SCANFOR(ONE) !Finds survivors.
LOGICAL ONE !Perhaps only the first is desired.
INTEGER I,P !Assistants.
DO I = 0,LB !Scan the current state.
IF (BIT(I)) THEN !Is this one unsullied?
P = N0 + 1 + 2*I !Yes! This is its value.
IF (P.LE.PRIME(NP)) CYCLE !But we may have it already.
IF (NP.GE.LASTP) RETURN !Whoops, perhaps too many!
NP = NP + 1 !But if not, another new prime!
PRIME(NP) = P !So, stash it. Extract now just this one.
IF (ONE) RETURN !Later candidates may yet be unzapped.
END IF !So much for that value.
END DO !On to the next.
END SUBROUTINE SCANFOR !An odd IF allows for two usages in one routine.
END SUBROUTINE PREPARE PRIMES !Faster than reading from a disc file?
END MODULE ERATOSTHENES !Certainly, less storage is required this way.
PROGRAM PRIMORIAL !Simple enough, with some assistants.
USE ERATOSTHENES !Though probably not as he expected.
USE BIGNUMBERS !Just so.
TYPE(BIGNUM) B !I'll have one.
INTEGER P,MARK !Step stuff.
INTEGER E,D !Assistants for the floating-point analogue...
INTEGER TASTE,IT !Additional stuff for its rounding.
PARAMETER (TASTE = 8/BIGORDER) !Sufficient digits to show.
INTEGER LEAD(TASTE) !With a struggle.
REAL T0,T1 !Some CPU time attempts.
REAL*4 F4 !I'll also have a go via logs.
REAL*8 F8 !In two precisions.
INTEGER I4,I8 !Not much hope for single precision, though.
WRITE (6,1) LASTP,BIGBASE !Announce.
1 FORMAT ("Calculates primorial numbers up to prime ",I0,
1 ", working in base ",I0)
CALL PREPARE PRIMES !First, catch your rabbit.
Commence prime mashing.
100 B.LAST = 1 !Begin at the beginning.
B.DIGIT(1) = 1 !With one.
DO P = 0,9 !Step up to the ninth prime, thus the first ten values as specified.
CALL BIGMULT(B,PRIME(P)) !Multiply by a possibly large integer.
WRITE (6,101) P,PRIME(P),P,B.DIGIT(B.LAST:1:-1) !Digits in Arabic/Hindu order.
101 FORMAT ("Prime(",I0,") = ",I0,", Primorial(",I0,") = ",
1 I0,9I<BIGORDER>.<BIGORDER>,/,(10I<BIGORDER>.<BIGORDER>))
END DO !On to the next prime.
Convert to logarithmic striders.
CALL CPU_TIME(T0) !Start the clock.
MARK = 10 !To be remarked upon in passing.
DO P = 10,LASTP !Step through additional primes.
CALL BIGMULT(B,PRIME(P)) !Bigger, ever bigger the big number grows.
IF (P.EQ.MARK) THEN !A report point?
MARK = MARK*10 !Yes. Prepare to note the next.
CALL CPU_TIME(T1) !Where are we at?
E = (B.LAST - 1)*BIGORDER !Convert from 10**BIGORDER to base 10.
D = B.DIGIT(B.LAST) !Grab the high-order digit.
DO WHILE(D.GT.0) !It is not zero..
E = E + 1 !So it is at least one base ten digit.
D = D/10 !Snip.
END DO !And perhaps there will be more.
Contemplate the rounding of the floating-point analogue.
I4 = MIN(TASTE,B.LAST) !I'm looking to taste the top digits.
LEAD(1:I4) = B.DIGIT(B.LAST:B.LAST - I4 + 1:-1) !Reverse, to have normal order.
IF (B.LAST.GT.TASTE) THEN !Are there even more digits?
IT = I4 !Yes. This is now the low-order digit tasted.
D = 0 !We should consider rounding up.
IF (B.DIGIT(B.LAST - I4).GE.BIGBASE/2) D = 1 !If the next digit is big enough.
DO WHILE (D.GT.0) !Spread the carry.
D = 0 !This one is used up.
LEAD(IT) = LEAD(IT) + 1 !Thusly.
IF (LEAD(IT).GT.BIGBASE) THEN !But, maybe, overflow!
IF (IT.GT.1) THEN !Is there a higher-order to carry to?
LEAD(IT) = LEAD(IT) - BIGBASE !Yes!
IT = IT - 1 !Step back to it,
D = 1 !Reassert a carry.
END IF !But only if there was a recipient available.
END IF !If not, the carry will still be zero.
END DO !And the loop won't continue.
END IF !So, no test for IT > 0 in a compound "while".
Cast forth the results.
WRITE (6,102) P,PRIME(P), !Name the step and its prime.
1 P,LEAD(1:I4),E, !The step and the leading few DIGIT of its primorial.
2 T1 - T0 !CPU advance.
102 FORMAT ("Prime(",I0,") = ",I0,", Primorial(",I0,") ~ 0." !Approximately.
1 I0,<I4 - 1>I<BIGORDER>.<BIGORDER>,"E+",I0, !No lead zero digits, then with lead zero digits.
2 T80,F12.3," seconds.") !Append some CPU time information.
T0 = T1 !Ready for the next popup.
END IF !So much for a report.
END DO !On to the next prime.
Chew some logarithms.
110 WRITE (6,111) !Some explanation.
111 FORMAT (/,"Via summing logarithms: Single Double") !Ah, layout.
MARK = 10 !Start somewhere interesting.
112 F4 = SUM(LOG10( FLOAT(PRIME(1:MARK)))) !Whee!
F8 = SUM(LOG10(DFLOAT(PRIME(1:MARK)))) !Generic function names too.
I4 = F4 !Grab the integer part.
I8 = F8 !The idea being to isolate the fractional part.
WRITE (6,113) MARK,"10",10**(F4 - I4),I4 + 1,10**(F8 - I8),I8 + 1 !Reconstitute the number in extended E-format.
113 FORMAT (I8,"#...in base ",A2,-1PF9.5,"E+",I0,T40,-1PF13.7,"E+",I0) !As if via E-format.
F4 = SUM(LOG( FLOAT(PRIME(1:MARK))))/LOG(10.0) !Do it again in Naperian logs.
F8 = SUM(LOG(DFLOAT(PRIME(1:MARK))))/LOG(10D0) !Perhaps more accurately?
I4 = F4
I8 = F8
WRITE (6,113) MARK,"e ",10**(F4 - I4),I4 + 1,10**(F8 - I8),I8 + 1 !We'll see.
MARK = MARK*10 !The next reporting point.
IF (MARK.LE.LASTP) GO TO 112 !Are we there yet?
END !So much for that.
Results
For the smaller numbers, the number of digits in the results is apparent. For the larger values, it is given by the value of the exponent part.
Calculates primorial numbers up to prime 1000000, working in base 100
Prime(0) = 1, Primorial(0) = 1
Prime(1) = 2, Primorial(1) = 2
Prime(2) = 3, Primorial(2) = 6
Prime(3) = 5, Primorial(3) = 30
Prime(4) = 7, Primorial(4) = 210
Prime(5) = 11, Primorial(5) = 2310
Prime(6) = 13, Primorial(6) = 30030
Prime(7) = 17, Primorial(7) = 510510
Prime(8) = 19, Primorial(8) = 9699690
Prime(9) = 23, Primorial(9) = 223092870
Prime(10) = 29, Primorial(10) ~ 0.64696932E+10 0.000 seconds.
Prime(100) = 541, Primorial(100) ~ 0.47119308E+220 0.000 seconds.
Prime(1000) = 7919, Primorial(1000) ~ 0.6786296E+3393 0.000 seconds.
Prime(10000) = 104729, Primorial(10000) ~ 0.9063360E+45337 0.875 seconds.
Prime(100000) = 1299709, Primorial(100000) ~ 0.1909604E+563921 121.547 seconds.
Prime(1000000) = 15485863, Primorial(1000000) ~ 0.1147175E+6722809 17756.797 seconds.
Via summing logarithms: Single Double
10#...in base 10 0.64697E+10 0.6469693E+10
10#...in base e 0.64697E+10 0.6469693E+10
100#...in base 10 0.47123E+220 0.4711931E+220
100#...in base e 0.47120E+220 0.4711931E+220
1000#...in base 10 0.67887E+3393 0.6786296E+3393
1000#...in base e 0.67849E+3393 0.6786296E+3393
10000#...in base 10 0.10650E+45338 0.9063360E+45337
10000#...in base e 0.82047E+45337 0.9063360E+45337
100000#...in base 10 0.15399E+563940 0.1909604E+563921
100000#...in base e 0.48697E+563884 0.1909604E+563921
1000000#...in base 10 0.31623E+669099 0.1147174E+6722809
1000000#...in base e 0.10000E+669683 0.1147175E+6722809
This is a rather severe demonstration of a faster-than-exponential growth rate in CPU time - the problem size increases exponentially as 10, 100, 1000, 10000, etc. (so that another 90, 900, 9000, etc. multiplies must be done) but the CPU time increases by more than a factor of ten. This is on an AMD FX6300 "six" core cpu at 3·76GHz that is steadily running the Great Internet Mersenne Prime computation on all six. Despite being at "idle" priority, it is presumably their contention for memory access that roughly triples the running time of any other individual computation.
Progress messages would have been helpful, but prior experience prompts caution: an attempt at having a collection of subprograms each accumulate their own CPU time usage for later assessment caused a worse than tenfold increase in CPU time for test runs and bizarre values as well. Another attempt at a "wait until a given time of day" that in the absence of access to any system routine for that became a loop on a call to the time-of-day routine promptly crashed the system - presumably due to too many (software-instituted) interrupts per second. Windows...
By contrast, there is only a slight pause for the even the largest summation of logarithms. Single precision demonstrates again its limited precision and 1000# suggests that base e logarithms are computed slightly more accurately. The double-precision results are the same as the first few decimal digits of the exact calculation, once some effort had been expended to have the values shown rounded correctly.
FreeBASIC
' version 22-09-2015
' compile with: fbc -s console
Const As UInteger Base_ = 1000000000
ReDim Shared As UInteger primes()
Sub sieve(need As UInteger)
' estimate is to high, but ensures that we have enough primes
Dim As UInteger max = need * (Log(need) + Log(Log(need)))
Dim As UInteger t = 1 ,x , x2
Dim As Byte p(max)
ReDim primes (need + need \ 3) ' we trim the array later
primes(0) = 1 ' by definition
primes(1) = 2 ' first prime, the only even prime
' only consider the odd number
For x = 3 To Sqr(max) Step 2
If p(x) = 0 Then
For x2 = x * x To max Step x * 2
p(x2) = 1
Next
End If
Next
' move found primes to array
For x = 3 To max Step 2
If p(x) = 0 Then
t += 1
primes(t) = x
EndIf
Next
'ReDim Preserve primes(t)
ReDim Preserve primes(need)
End Sub
' ------=< MAIN >=------
Dim As UInteger n, i, pow, primorial
Dim As String str_out, buffer = Space(10)
Dim As UInteger max = 100000 ' maximum number of primes we need
sieve(max)
primorial = 1
Print
For n = 0 To 9
primorial = primorial * primes(n)
Print Using " primorial(#) ="; n;
RSet buffer, Str(primorial)
str_out = buffer
Print str_out
Next
' could use GMP, but why not make are own big integer routine
Dim As UInteger bigint(max), first = max, last = max
Dim As UInteger l, p, carry, low = 9, high = 10
Dim As ULongInt result
Dim As UInteger Ptr big_i
' start at the back, number grows to the left like normal number
bigint(last) = primorial
Print
For pow = 0 To Len(Str(max)) -2
If pow > 0 Then
low = high
high = high * 10
End If
For n = low + 1 To high
carry = 0
big_i = @bigint(last)
For i = last To first Step -1
result = CULngInt(primes(n)) * *big_i + carry
carry = result \ Base_
*big_i = result - carry * Base_
big_i = big_i -1
Next i
If carry <> 0 Then
first = first -1
*big_i = carry
End If
Next n
l = Len(Str(bigint(first))) + (last - first) * 9
Print " primorial("; high; ") has "; l ;" digits"
Next pow
' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
- Output:
primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 primorial(10) has 10 digits primorial(100) has 220 digits primorial(1000) has 3393 digits primorial(10000) has 45337 digits primorial(100000) has 563921 digits
Frink
Simple version
primorial[n] := product[first[primes[], n]]
for n = 0 to 9
println["primorial[$n] = " + primorial[n]]
for n = [10, 100, 1000, 10000, 100000, million]
println["Length of primorial $n is " + length[toString[primorial[n]]] + " decimal digits."]- Output:
primorial[0] = 1 primorial[1] = 2 primorial[2] = 6 primorial[3] = 30 primorial[4] = 210 primorial[5] = 2310 primorial[6] = 30030 primorial[7] = 510510 primorial[8] = 9699690 primorial[9] = 223092870 Length of primorial 10 is 10 decimal digits. Length of primorial 100 is 220 decimal digits. Length of primorial 1000 is 3393 decimal digits. Length of primorial 10000 is 45337 decimal digits. Length of primorial 100000 is 563921 decimal digits. Length of primorial 1000000 is 6722809 decimal digits.
Binary splitting
The program above can be made much faster by using a "binary splitting" algorithm that recursively breaks the number range into halves, and then multiplies these smaller numbers together, which will be faster when run under a Java Virtual Machine version 1.8 and later which have a better-than-O(n2) multiply operation. (Did you know that Frink's author submitted the changes to Java that made BigInteger and BigDecimal much faster with lots of digits?)
/** Calculate the primes and set up for the recursive call. */
primorialSplitting[n] :=
{
if n == 0
return 1
primes = array[first[primes[], n]]
return primorialSplitting[n, 0, n-1, primes]
}
/** The actual recursive algorithm. */
primorialSplitting[n, begin, end, primes] :=
{
range = (end-begin)
if range >= 2
{
middle = (begin + end) div 2
return primorialSplitting[n, begin, middle, primes] * primorialSplitting[n, middle+1, end, primes]
}
if range == 1
return primes@begin * primes@end
if range == 0
return primes@begin
}
for n = 0 to 9
println["primorial[$n] = " + primorialSplitting[n]]
for n = [10, 100, 1000, 10000, 100000, million]
println["Length of primorial $n is " + length[toString[primorialSplitting[n]]] + " decimal digits."]Go
Basic multiplication
Since this task isn't specifically about generating primes, a fast external package is used to generate the primes. The Go standard math/big package is used to multiply these as exact integers.
package main
import (
"fmt"
"math/big"
"time"
"github.com/jbarham/primegen.go"
)
func main() {
start := time.Now()
pg := primegen.New()
var i uint64
p := big.NewInt(1)
tmp := new(big.Int)
for i <= 9 {
fmt.Printf("primorial(%v) = %v\n", i, p)
i++
p = p.Mul(p, tmp.SetUint64(pg.Next()))
}
for _, j := range []uint64{1e1, 1e2, 1e3, 1e4, 1e5, 1e6} {
for i < j {
i++
p = p.Mul(p, tmp.SetUint64(pg.Next()))
}
fmt.Printf("primorial(%v) has %v digits", i, len(p.String()))
fmt.Printf("\t(after %v)\n", time.Since(start))
}
}
- Output:
primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 primorial(10) has 10 digits (after 367.273µs) primorial(100) has 220 digits (after 8.460177ms) primorial(1000) has 3393 digits (after 8.760826ms) primorial(10000) has 45337 digits (after 26.838309ms) primorial(100000) has 563921 digits (after 2.049290216s) primorial(1000000) has 6722809 digits (after 4m19.931513819s)
Since this includes timing information this is also relevant:
% go version ; sysctl -n hw.model go version go1.4.2 freebsd/amd64 Intel(R) Xeon(R) CPU E3-1245 v3 @ 3.40GHz
Product tree
As noted elsewhere, it is far quicker to use a product tree to multiply the primes together. In the following version, the vecprod() function follows the lines of the Phix example with a GMP wrapper (rather than Go's native big.Int type) being used for extra speed.
package main
import (
"fmt"
"github.com/jbarham/primegen"
big "github.com/ncw/gmp"
"time"
)
func vecprod(primes []uint64) *big.Int {
if len(primes) == 0 {
return big.NewInt(1)
}
s := make([]*big.Int, len(primes))
le := len(s)
for i := 0; i < le; i++ {
s[i] = new(big.Int).SetUint64(primes[i])
}
for le > 1 {
for i := 0; i < le/2; i++ {
s[i].Mul(s[i], s[le-i-1])
}
c := le / 2
if le&1 == 1 {
c++
}
s = s[0:c]
le = c
}
return s[0]
}
func main() {
start := time.Now()
pg := primegen.New()
var primes []uint64
for i := uint64(0); i < 1e6; i++ {
primes = append(primes, pg.Next())
}
for i := 0; i < 10; i++ {
fmt.Printf("primorial(%d) = %d\n", i, vecprod(primes[0:i]))
}
fmt.Println()
for _, i := range []uint64{1e1, 1e2, 1e3, 1e4, 1e5, 1e6} {
fmt.Printf("primorial(%d) has length %d\n", i, len(vecprod(primes[0:i]).String()))
}
fmt.Printf("\nTook %s\n", time.Since(start))
}
- Output:
Timed using a Core i7-8565U machine running Go 1.14.1 on Ubuntu 18.04:
primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 primorial(10) has length 10 primorial(100) has length 220 primorial(1000) has length 3393 primorial(10000) has length 45337 primorial(100000) has length 563921 primorial(1000000) has length 6722809 Took 2.11913871s
Haskell
import Control.Arrow ((&&&))
import Data.List (scanl1, foldl1')
getNthPrimorial :: Int -> Integer
getNthPrimorial n = foldl1' (*) (take n primes)
primes :: [Integer]
primes = 2 : filter isPrime [3,5..]
isPrime :: Integer -> Bool
isPrime = isPrime_ primes
where isPrime_ :: [Integer] -> Integer -> Bool
isPrime_ (p:ps) n
| p * p > n = True
| n `mod` p == 0 = False
| otherwise = isPrime_ ps n
primorials :: [Integer]
primorials = 1 : scanl1 (*) primes
main :: IO ()
main = do
-- Print the first 10 primorial numbers
let firstTen = take 10 primorials
putStrLn $ "The first 10 primorial numbers are: " ++ show firstTen
-- Show the length of the primorials with index 10^[1..6]
let powersOfTen = [1..6]
primorialTens = map (id &&& (length . show . getNthPrimorial . (10^))) powersOfTen
calculate = mapM_ (\(a,b) -> putStrLn $ "Primorial(10^"++show a++") has "++show b++" digits")
calculate primorialTens
- Output:
The first 10 primorial numbers are: [1,2,6,30,210,2310,30030,510510,9699690,223092870] Primorial(10^1) has 10 digits Primorial(10^2) has 220 digits Primorial(10^3) has 3393 digits Primorial(10^4) has 45337 digits Primorial(10^5) has 563921 digits Primorial(10^6) has 6722809 digits
J
Implementation:
primorial=:*/@:p:@i."0
Task examples:
primorial i. 10 NB. first 10 primorial numbers
1 2 6 30 210 2310 30030 510510 9699690 223092870
#":primorial 10x NB. lengths (of decimal representations)...
10
#":primorial 100x
220
#":primorial 1000x
3393
#":primorial 10000x
45337
#":primorial 100000x
563921
#":primorial 1000000x
6722809
Note that the x suffix on a decimal number indicates that calculations should be exact (what the lisp community has dubbed "bignums"). This is a bit slow (internally, every number winds up being represented as a list of what might be thought of as "digits" in a very large base), but it gets the job done.
Java
import java.math.BigInteger;
public class PrimorialNumbers {
final static int sieveLimit = 1300_000;
static boolean[] notPrime = sieve(sieveLimit);
public static void main(String[] args) {
for (int i = 0; i < 10; i++)
System.out.printf("primorial(%d): %d%n", i, primorial(i));
for (int i = 1; i < 6; i++) {
int len = primorial((int) Math.pow(10, i)).toString().length();
System.out.printf("primorial(10^%d) has length %d%n", i, len);
}
}
static BigInteger primorial(int n) {
if (n == 0)
return BigInteger.ONE;
BigInteger result = BigInteger.ONE;
for (int i = 0; i < sieveLimit && n > 0; i++) {
if (notPrime[i])
continue;
result = result.multiply(BigInteger.valueOf(i));
n--;
}
return result;
}
public static boolean[] sieve(int limit) {
boolean[] composite = new boolean[limit];
composite[0] = composite[1] = true;
int max = (int) Math.sqrt(limit);
for (int n = 2; n <= max; n++) {
if (!composite[n]) {
for (int k = n * n; k < limit; k += n) {
composite[k] = true;
}
}
}
return composite;
}
}
primorial(0): 1 primorial(1): 2 primorial(2): 6 primorial(3): 30 primorial(4): 210 primorial(5): 2310 primorial(6): 30030 primorial(7): 510510 primorial(8): 9699690 primorial(9): 223092870 primorial(10^1) has length 10 primorial(10^2) has length 220 primorial(10^3) has length 3393 primorial(10^4) has length 45337 primorial(10^5) has length 563921
JavaScript
Exact integers
Uses the native BigInt type availiable in Node.js. As this takes a while, some values between 10 000 and 100 000 are shown to help suggest it is still doing stuff...
{ // calculate and show some primorial numbers
'use strict'
let primorial = 1n
let prime = 1n
let pn = []
const maxNumber = 2000000
let isPrime = []
for( let i = 1; i <= maxNumber; i ++ ){ isPrime[ i ] = i % 2 != 0 }
isPrime[ 1 ] = false
isPrime[ 2 ] = true
const rootMaxNumber = Math.floor( Math.sqrt( maxNumber ) )
for( let s = 3; s <= rootMaxNumber; s += 2 ){
if( isPrime[ s ] ){
for( let p = s * s; p <= maxNumber; p += s ){ isPrime[ p ] = false }
}
}
const primeMax = 100000
pn[ 0 ] = 1
let nextToShow = 10
for( let i = 1; i <= primeMax; i ++ ){
// find the next prime
prime += 1n
while( ! isPrime[ prime ] ){ prime += 1n }
primorial *= prime
if( i < 10 ){
pn[ i ] = primorial
}
else if( i == nextToShow ){
if( nextToShow < 10000 ){
nextToShow *= 10
}
else{
nextToShow += 10000
}
if( i == 10 ){ console.log( "primorials 0-9: ", pn.toString() ) }
// show the number of digits in the primorial
let p = primorial
let length = 0
while( p > 0 ){
length += 1
p /= 10n
}
console.log( "length of primorial " + i + " is "+ length )
}
}
}
- Output:
primorials 0-9: 1,2,6,30,210,2310,30030,510510,9699690,223092870 length of primorial 10 is 10 length of primorial 100 is 220 length of primorial 1000 is 3393 length of primorial 10000 is 45337 length of primorial 20000 is 97389 length of primorial 30000 is 151937 length of primorial 40000 is 208100 length of primorial 50000 is 265460 length of primorial 60000 is 323772 length of primorial 70000 is 382876 length of primorial 80000 is 442655 length of primorial 90000 is 503026 length of primorial 100000 is 563921
Runtime is around 10 minutes on the Windows 11 laptop I'm using (too long for TIO.RUN).
Using a reduced number of digits
A modification of the Exact integers version, also using the native BigInt type. As we only need to show digit counts, this deviates from the task requirement to use exact integers by only keeping around 500 leading digits., which is MUCH faster.
{ // calculate and show some primorial numbers
'use strict'
let primorial = 1n
let prime = 1n
let pn = []
const maxNumber = 2000000
let isPrime = []
for( let i = 1; i <= maxNumber; i ++ ){ isPrime[ i ] = i % 2 != 0 }
isPrime[ 1 ] = false
isPrime[ 2 ] = true
const rootMaxNumber = Math.floor( Math.sqrt( maxNumber ) )
for( let s = 3; s <= rootMaxNumber; s += 2 ){
if( isPrime[ s ] ){
for( let p = s * s; p <= maxNumber; p += s ){ isPrime[ p ] = false }
}
}
const primeMax = 100000
pn[ 0 ] = 1
let nextToShow = 10
let reducedDigits = 0
const n500 = 10n**500n
const n1000 = 10n**1000n
for( let i = 1; i <= primeMax; i ++ ){
// find the next prime
prime += 1n
while( ! isPrime[ prime ] ){ prime += 1n }
primorial *= prime
if( i < 10 ){
pn[ i ] = primorial
}
else if( i == nextToShow ){
if( nextToShow < 10000 ){
nextToShow *= 10
}
else{
nextToShow += 10000
}
if( i == 10 ){ console.log( "primorials 0-9: ", pn.toString() ) }
// show the number of digits in the primorial
let p = primorial
let length = 0
while( p > 0 ){
length += 1
p /= 10n
}
console.log( "length of primorial " + i + " is "+ ( reducedDigits + length ) )
}
if( primorial > n1000 ){
// the number has more than 1000 digits - reduce it to 500-ish
primorial /= n500
reducedDigits += 500
}
}
}
- Output:
primorials 0-9: 1,2,6,30,210,2310,30030,510510,9699690,223092870 length of primorial 10 is 10 length of primorial 100 is 220 length of primorial 1000 is 3393 length of primorial 10000 is 45337 length of primorial 20000 is 97389 length of primorial 30000 is 151937 length of primorial 40000 is 208100 length of primorial 50000 is 265460 length of primorial 60000 is 323772 length of primorial 70000 is 382876 length of primorial 80000 is 442655 length of primorial 90000 is 503026 length of primorial 100000 is 563921
Runtime is around 1 second on TIO.RUN.
jq
Works with gojq, the Go implementation of jq
See Erdős-primes#jq for a suitable definition of `is_prime` as used here.
def primes:
2, range(3; infinite; 2) | select(is_prime);
# generate an infinite stream of primorials beginning with primorial(0)
def primorials:
0, foreach primes as $p (1; .*$p; .);
"The first ten primorial numbers are:",
limit(10; primorials),
"\nThe primorials with the given index have the lengths shown:",
([10, 100, 1000, 10000, 100000] as $sample
| limit($sample|length;
foreach primes as $p ([0,1]; # [index, primorial]
.[0]+=1 | .[1] *= $p;
select(.[0]|IN($sample[])) | [.[0], (.[1]|tostring|length)] ) ))- Output:
The first ten primorial numbers are: 0 2 6 30 210 2310 30030 510510 9699690 223092870 The primorials with the given index have the lengths shown: [10,10] [100,220] [1000,3393] [10000,45337] [100000,563921]
Julia
using Primes
primelist = primes(300000001) # primes to 30 million
primorial(n) = foldr(*, primelist[1:n], init=BigInt(1))
println("The first ten primorials are: $([primorial(n) for n in 1:10])")
for i in 1:6
n = 10^i
p = primorial(n)
plen = Int(floor(log10(p))) + 1
println("primorial($n) has length $plen digits in base 10.")
end
- Output:
The first ten primorials are: BigInt[2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870, 6469693230] primorial(10) has length 10 digits in base 10. primorial(100) has length 220 digits in base 10. primorial(1000) has length 3393 digits in base 10. primorial(10000) has length 45337 digits in base 10. primorial(100000) has length 563921 digits in base 10. primorial(1000000) has length 6722809 digits in base 10.
Kotlin
// version 1.0.6
import java.math.BigInteger
const val LIMIT = 1000000 // expect a run time of about 20 minutes on a typical laptop
fun isPrime(n: Int): Boolean {
if (n < 2) return false
if (n % 2 == 0) return n == 2
if (n % 3 == 0) return n == 3
var d : Int = 5
while (d * d <= n) {
if (n % d == 0) return false
d += 2
if (n % d == 0) return false
d += 4
}
return true
}
fun countDigits(bi: BigInteger): Int = bi.toString().length
fun main(args: Array<String>) {
println("Primorial(0) = 1")
println("Primorial(1) = 2")
var count = 1
var p = 3
var prod = BigInteger.valueOf(2)
var target = 10
while(true) {
if (isPrime(p)) {
count++
prod *= BigInteger.valueOf(p.toLong())
if (count < 10) {
println("Primorial($count) = $prod")
if (count == 9) println()
}
else if (count == target) {
println("Primorial($target) has ${countDigits(prod)} digits")
if (count == LIMIT) break
target *= 10
}
}
p += 2
}
}
- Output:
Primorial(0) = 1 Primorial(1) = 2 Primorial(2) = 6 Primorial(3) = 30 Primorial(4) = 210 Primorial(5) = 2310 Primorial(6) = 30030 Primorial(7) = 510510 Primorial(8) = 9699690 Primorial(9) = 223092870 Primorial(10) has 10 digits Primorial(100) has 220 digits Primorial(1000) has 3393 digits Primorial(10000) has 45337 digits Primorial(100000) has 563921 digits Primorial(1000000) has 6722809 digits
Lingo
For generating the list of primes an auto-extending Sieve of Eratosthenes lib is used (fast), and for the larger primorials a simple custom big-integer lib (very slow).
-- libs
sieve = script("math.primes").new()
bigint = script("bigint").new()
cnt = 1000 * 100
primes = sieve.getNPrimes(cnt)
pr = 1
put "Primorial 0: " & pr
repeat with i = 1 to 9
pr = pr*primes[i]
put "Primorial " & i & ": " & pr
end repeat
pow10 = 10
repeat with i = 10 to cnt
pr = bigint.mul(pr, primes[i])
if i mod pow10=0 then
put "Primorial " & i & " has length: " & pr.length
pow10 = pow10 * 10
end if
end repeat- Output:
-- "Primorial 0: 1" -- "Primorial 1: 2" -- "Primorial 2: 6" -- "Primorial 3: 30" -- "Primorial 4: 210" -- "Primorial 5: 2310" -- "Primorial 6: 30030" -- "Primorial 7: 510510" -- "Primorial 8: 9699690" -- "Primorial 9: 223092870" -- "Primorial 10 has length: 10" -- "Primorial 100 has length: 220" -- "Primorial 1000 has length: 3393" -- "Primorial 10000 has length: 45337" -- "Primorial 100000 has length: 563921"
Maple
with(NumberTheory):
primorial := proc(n::integer)
local total := 1:
local count:
for count from 1 to n do
total *= ithprime(count):
end:
return total;
end proc:
primorialDigits := proc(n::integer)
local logSum := 0:
local count:
for count from 1 to n do
logSum += log10(ithprime(count)):
end:
return ceil(logSum);
end proc:
print("The first 10 primorial numbers");
for count from 0 to 9 do
cat("primorial(", count, ") = ", primorial(count))
end;
for expon from 1 to 5 do
cat("primorial(", 10^expon, ") has ", primorialDigits(10^expon), " digits");
end;- Output:
"The first 10 primorial numbers""primorial(0) = 1""primorial(1) = 2""primorial(2) = 6""primorial(3) = 30""primorial(4) = 210""primorial(5) = 2310""primorial(6) = 30030""primorial(7) = 510510""primorial(8) = 9699690""primorial(9) = 223092870""primorial(10) has 10 digits""primorial(100) has 220 digits""primorial(1000) has 3393 digits""primorial(10000) has 45337 digits""primorial(100000) has 563921 digits"
Mathematica /Wolfram Language
The first 10 primorial numbers are:
FoldList[Times, 1, Prime @ Range @ 9]
- Output:
{1,2,6,30,210,2310,30030,510510,9699690,223092870}
From the first million primes:
primes = Prime @ Range[10^6];
define the primorial:
primorial[n_]:= Times @@ primes[[;;n]]
The lengths of selected primorial numbers are:
Grid@Table[{"primorial(10^" <> ToString[n] <> ") has ",
{timing,answer}=AbsoluteTiming[IntegerLength@primorial[10^n]];answer,
" digits in "<>ToString@timing<>" seconds"}, {n,6}]
- Output:
primorial(10^1) has 10 digits in 0.000055 seconds primorial(10^2) has 220 digits in 0.00007 seconds primorial(10^3) has 3393 digits in 0.000425 seconds primorial(10^4) has 45337 digits in 0.006795 seconds primorial(10^5) has 563921 digits in 0.076371 seconds primorial(10^6) has 6722809 digits in 1.29307 seconds
Nickle
library "prime_sieve.5c"
# For 1 million primes
# int val = 15485867;
int val = 1299743;
int start = millis();
int [*] primes = PrimeSieve::primes(val);
printf("%d primes (%d) in %dms\n", dim(primes), primes[dim(primes)-1], millis() - start);
int primorial(int n) {
if (n == 0) return 1;
if (n == 1) return 2;
int v = 2;
for (int i = 2; i <= n; i++) {
v *= primes[i-2];
}
return v;
}
for (int i = 0; i < 10; i++) {
printf("primorial(%d) = %d\n", i, primorial(i));
}
for (int i = 1; i < 6; i++) {
start = millis();
int p = 10**i;
int pn = primorial(p);
int digits = floor(Math::log10(pn)) + 1;
printf("primorial(%d) has %d digits, in %dms\n", p, digits, millis() - start);
}
- Output:
prompt$ nickle primorial.5c 100001 primes (1299743) in 1838ms primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 primorial(10) has 10 digits, in 4ms primorial(100) has 220 digits, in 2ms primorial(1000) has 3393 digits, in 3ms primorial(10000) has 45337 digits, in 128ms primorial(100000) has 563921 digits, in 12304ms
Nim
Single thread
We use a sieve of Erathosthenes to generate the primes, but with only the odd numbers to reduce memory usage. Even for one millions primes, the execution time is negligible (0.12s)
To compute the primorial numbers, we use an iterator. Performance is good with 1.1s for primorial(10000). But for one million, this takes much more time.
import times
let t0 = cpuTime()
####################################################################################################
# Build list of primes.
const
NPrimes = 1_000_000
N = 16 * NPrimes
var sieve: array[(N - 1) div 2 + 1, bool] # False (default) means prime.
for i, composite in sieve:
if not composite:
let n = 2 * i + 3
for k in countup(n * n, N, 2 * n):
sieve[(k - 3) div 2] = true
var primes = @[2]
for i, composite in sieve:
if not composite:
primes.add 2 * i + 3
if primes.len < NPrimes:
quit "Not enough primes. Please, increase value of N."
####################################################################################################
# Compute primorial.
import strformat
import bignum
const LastToPrint = NPrimes
iterator primorials(): Int =
## Yield successive primorial numbers.
var prim = newInt(1)
yield prim
for p in primes:
prim *= p
yield prim
var n = 0
for prim in primorials():
echo &"primorial({n}) = {prim}"
inc n
if n == 10: break
n = 0
var nextToPrint = 10
for prim in primorials():
if n == nextToPrint:
echo &"primorial({n}) has {($prim).len} digits"
if nextToPrint == LastToPrint: break
nextToPrint *= 10
inc n
echo ""
echo &"Total time: {cpuTime() - t0:.2f} s"
- Output:
primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 primorial(10) has 10 digits primorial(100) has 220 digits primorial(1000) has 3393 digits primorial(10000) has 45337 digits primorial(100000) has 563921 digits primorial(1000000) has 6722809 digits Total time: 127.65 s
CPU: Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz, 2607 MHz with 8GB RAM.
Multiple threads
With several threads, we cannot use an iterator anymore. For N threads (workers), we split the list of primes in N sublists and let each worker compute a product. We terminate with a final multiplication of the N results. With four workers, we get the result in about 10 seconds. With eight workers, we get it in about five seconds.
Note that to get execution time, we can no longer use “cpuTime” as it returns only the CPU time used by the main thread. So, we use “getMonoTime” instead.
import std/monotimes
let t0 = getMonoTime()
####################################################################################################
# Build list of primes.
const
NPrimes = 1_000_000
N = 16 * NPrimes
var sieve: array[(N - 1) div 2 + 1, bool] # False (default) means prime.
for i, composite in sieve:
if not composite:
let n = 2 * i + 3
for k in countup(n * n, N, 2 * n):
sieve[(k - 3) div 2] = true
var primes = @[2]
for i, composite in sieve:
if not composite:
primes.add 2 * i + 3
if primes.len < NPrimes:
quit "Not enough primes. Please, increase value of N."
####################################################################################################
# Compute primorial.
import strformat, threadpool
import bignum
const NWorkers = 8
proc computeProduct(a: openArray[int]): Int =
result = newInt(1)
for n in a: result *= n
proc primorial(n: int): Int =
if n == 0: return newInt(1)
# Prepare sublists.
var input: array[NWorkers, seq[int]]
for i in 0..<n:
input[i mod NWorkers].add primes[i]
# Spawn workers and get partial products.
var responses: array[NWorkers, FlowVar[Int]]
for i in 0..<NWorkers:
responses[i] = spawn computeProduct(input[i])
# Compute final product.
result = ^responses[0]
for i in 1..<NWorkers:
result *= ^responses[i]
for n in 0..9:
echo &"primorial({n}) = {primorial(n)}"
for n in [10, 100, 1_000, 10_000, 1_000_000]:
echo &"primorial({n}) has {len($primorial(n))} digits"
echo ""
echo &"Total time: {(getMonoTime() - t0)}"
- Output:
primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 primorial(10) has 10 digits primorial(100) has 220 digits primorial(1000) has 3393 digits primorial(10000) has 45337 digits primorial(1000000) has 6722809 digits Total time: (seconds: 5, nanosecond: 270238339)
PARI/GP
nthprimorial(n)=prod(i=1,n,prime(i));
vector(10,i,nthprimorial(i-1))
vector(5,n,#Str(nthprimorial(10^n)))
#Str(nthprimorial(10^6))- Output:
%1 = [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870] %2 = [10, 220, 3393, 45337, 563921] %3 = 433637
With vecprod
Pari/GP 2.11 added the vecprod command, which makes a significantly faster version possible. We use primes(n) to get a vector of the first n primes, which vecprod multiplies using a product tree.
nthprimorial(n)=vecprod(primes(n));
vector(6,n,#Str(nthprimorial(10^n)))- Output:
[10, 220, 3393, 45337, 563921, 6722809] ? ## *** last result computed in 1,663 ms.
Pascal
Using unit like GMP. MPArith of wolfgang-ehrhardt.de mp_primorial {-Primorial of n; a = n# = product of primes <= n.} so i sieve to get the right n
{$H+}
uses
sysutils,mp_types,mp_base,mp_prime,mp_numth;
var
x: mp_int;
t0,t1: TDateTime;
s: AnsiString;
var
i,cnt : NativeInt;
ctx :TPrimeContext;
begin
mp_init(x);
cnt := 1;
i := 2;
FindFirstPrime32(i,ctx);
i := 10;
t0 := time;
repeat
repeat
FindNextPrime32(ctx);
inc(cnt);
until cnt = i;
mp_primorial(ctx.prime,x);
s:= mp_adecimal(x);
writeln('MaxPrime ',ctx.prime:10,length(s):8,' digits');
i := 10*i;
until i > 1000*1000;
t1 := time;
Writeln((t1-t0)*86400.0:10:3,' s');
end.
- Output:
MaxPrime 29 10 digits
MaxPrime 541 220 digits
MaxPrime 7919 3393 digits
MaxPrime 104729 45337 digits
MaxPrime 1299709 563921 digits
MaxPrime 15485863 6722809 digits
111.290 s
using i4330 3.5 Ghz Win7 32-bit fpc 2.6.4 a little outdated
uups: 174.468 s Linux 32-bit fpc 3.0.1
real 2m54.470s user 1m59.133s sys 0m55.188s!
alternative
getPrimorialExact uses multiplication to Base 1e9.Heavily inspired by FreeBASIC It is extremly slow.PrimorialExact (1000000) should take about 2400 seconds on my computer.See GO which takes 2 secs instead of my 38 secs ( 3.4 Ghz 64Bit vs 3.5 Ghz 32 Bit ). Tested x64 -> runtime 16m48 x64 div is substituted by mul and shift.Maybe first Mul then base convertion will be faster. Obviously GMP will do it by far better.
program Primorial;
{$IFDEF FPC} {$MODE DELPHI} {$ENDIF}
uses
sysutils;
var
primes : array[0..1000000] of LongInt;
procedure InitSieve;
const
HiSieve = 15485864;
var
sieve: array of boolean;
i, j: NativeInt;
Begin
setlength(sieve,HiSieve);
fillchar(sieve[0],HiSieve,chr(ord(True)));
For i := 2 to Trunc(sqrt(HiSieve)) do
IF sieve[i] then Begin
j := i*i;repeat sieve[j]:= false;inc(j,i);until j>= HiSieve-1;end;
i:= 2;j:= 1;
repeat
IF sieve[i] then begin primes[j]:= i;inc(j) end;
inc(i);
until i > HiSieve;
primes[0] := 1;setlength(sieve,0);
end;
function getPrimorial(n:NativeInt):Uint64;
Begin
result := ORD(n>=0);
IF (n >= 0) AND (n < 16) then
repeat result := result*primes[n]; dec(n); until n < 1;
end;
function getPrimorialDecDigits(n:NativeInt):NativeInt;
var
res: extended;
Begin
result := -1;
IF (n > 0) AND (n <= 1000*1000) then
Begin
res := 0;
repeat res := res+ln(primes[n]); dec(n); until n < 1;
result := trunc(res/ln(10))+1;
end;
end;
function getPrimorialExact(n:NativeInt):NativeInt;
const
LongWordDec = 1000000000;
var
MulArr : array of LongWord;
pMul : ^LongWord;
Mul1,prod,carry : Uint64;
i,j,ul : NativeInt;
begin
i := getPrimorialDecDigits(n) DIV 9 +10;
Setlength(MulArr,i);
Ul := 0;
MulArr[Ul]:= 1;
i := 1;
repeat
Mul1 := 1;
//Make Mul1 as large as possible
while (i<= n) AND ((LongWordDec DIV MUL1) >= primes[i]) do
Begin Mul1 := Mul1*primes[i]; inc(i); end;
carry := 0;
pMul := @MulArr[0];
For j := 0 to UL do
Begin
prod := Mul1*pMul^+Carry;
Carry := prod Div LongWordDec;
pMul^ := Prod - Carry*LongWordDec;
inc(pMul);
end;
IF Carry <> 0 then Begin inc(Ul);pMul^:= Carry; End;
until i> n;
//count digits
i := Ul*9;
Carry := MulArr[Ul];
repeat
Carry := Carry DIV 10;
inc(i);
until Carry = 0;
result := i;
end;
var
i: NativeInt;
Begin
InitSieve;
write('Primorial (0->9) ');
For i := 0 to 9 do
write(getPrimorial(i),',');
writeln(#8#32#13#10);
i:= 10;
repeat
writeln('Primorial (',i,') = digits ',
getPrimorialDecDigits(i),' digits');
i := i*10;
until i> 1000000;
writeln;
i:= 10;
repeat
writeln('PrimorialExact (',i,') = digits ',
getPrimorialExact(i),' digits');
i := i*10;
until i> 100000;
end.
- Output:
Primorial (0->9) 1,2,6,30,210,2310,30030,510510,9699690,223092870 Primorial (10) = digits 10 digits Primorial (100) = digits 220 digits Primorial (1000) = digits 3393 digits Primorial (10000) = digits 45337 digits Primorial (100000) = digits 563921 digits Primorial (1000000) = digits 6722809 digits //real 0m0.143s without PrimorialExact PrimorialExact (10) = digits 10 digits PrimorialExact (100) = digits 220 digits PrimorialExact (1000) = digits 3393 digits PrimorialExact (10000) = digits 45337 digits PrimorialExact (100000) = digits 563921 digits real 0m38.311s for x64 I tried it PrimorialExact (1000000) = digits 6722809 digits real 16m48.240s
Perl
use ntheory qw(pn_primorial);
say "First ten primorials: ", join ", ", map { pn_primorial($_) } 0..9;
say "primorial(10^$_) has ".(length pn_primorial(10**$_))." digits" for 1..6;
The pn_primorial function is smart enough to return a Math::BigInt object if the result won't fit into a native integer, so it all works out.
- Output:
First ten primorials: 1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870 primorial(10^1) has 10 digits primorial(10^2) has 220 digits primorial(10^3) has 3393 digits primorial(10^4) has 45337 digits primorial(10^5) has 563921 digits primorial(10^6) has 6722809 digits
Still using the library for the core activities, we can do the same in two steps. primes($n) returns a reference to a list of primes up to n, nth_prime($n) returns the nth prime, and vecprod does an efficient product tree. So for 10^6 we can do:
use ntheory ":all";
say length( vecprod( @{primes( nth_prime(10**6) )} ) );
returning the same result in only slightly more time. Both examples take under 4 seconds.
Phix
Multiplying the vector elements in pairs is much faster for essentially the same reason that merge sort is faster than insertion sort.
You could probably optimise the number of mpz_init() this invokes a fair bit, if you really wanted to.
with javascript_semantics include mpfr.e function vecprod(sequence s) if s={} then s = {mpz_init(1)} else for i=1 to length(s) do s[i] = mpz_init(s[i]) end for while length(s)>1 do for i=1 to floor(length(s)/2) do mpz_mul(s[i],s[i],s[-i]) end for s = s[1..ceil(length(s)/2)] end while end if return s[1] end function constant t0 = time(), max10 = iff(platform()=JS?4:6), tests = tagset(9,0)&sq_power(10,tagset(max10)), primes = get_primes(1_000_000) for i=1 to length(tests) do integer ti = tests[i] mpz primorial = vecprod(primes[1..ti]) string ps = iff(ti<10?sprintf("= %s",{mpz_get_str(primorial,comma_fill:=true)}) :sprintf("has %,d digits",mpz_sizeinbase(primorial,10))) printf(1,"Primorial(%,d) %s\n", {ti, ps}) end for ?elapsed(time()-t0)
- Output:
Primorial(0) = 1 Primorial(1) = 2 Primorial(2) = 6 Primorial(3) = 30 Primorial(4) = 210 Primorial(5) = 2,310 Primorial(6) = 30,030 Primorial(7) = 510,510 Primorial(8) = 9,699,690 Primorial(9) = 223,092,870 Primorial(10) has 10 digits Primorial(100) has 220 digits Primorial(1,000) has 3,393 digits Primorial(10,000) has 45,337 digits Primorial(100,000) has 563,921 digits Primorial(1,000,000) has 6,722,809 digits "6.2s"
Under pwa/p2js Primorial(10,000) is obtained in 2s but 10^5 appears to be well beyond reach.
PicoLisp
This code uses prime? and take functions from Extensible Prime Generator(PicoLisp)
(de prime? (N Lst)
(let S (sqrt N)
(for D Lst
(T (> D S) T)
(T (=0 (% N D)) NIL) ) ) )
(de take (N)
(let I 1
(make
(link 2)
(do (dec N)
(until (prime? (inc 'I 2) (made)))
(link I) ) ) ) )
# This is a simple approach to calculate primorial may not be the fastest one
(de primorial (N)
(apply * (take N)) )
#print 1st 10 primorial numbers
(for M 10 (prinl "primorial: "(primorial M)))
# print the length of primorial numbers.
[prinl (length (primorial (** 10 1)]
[prinl (length (primorial (** 10 2)]
[prinl (length (primorial (** 10 3)]
[prinl (length (primorial (** 10 4)]
#The last one takes a very long time to compute.
[prinl (length (primorial (** 10 5)]- Output:
primorial: 2 primorial: 6 primorial: 30 primorial: 210 primorial: 2310 primorial: 30030 primorial: 510510 primorial: 9699690 primorial: 223092870 primorial: 6469693230 10 220 3393 45337 563921
Python
Uses the pure python library pyprimes .
from pyprimes import nprimes
from functools import reduce
primelist = list(nprimes(1000001)) # [2, 3, 5, ...]
def primorial(n):
return reduce(int.__mul__, primelist[:n], 1)
if __name__ == '__main__':
print('First ten primorals:', [primorial(n) for n in range(10)])
for e in range(7):
n = 10**e
print('primorial(%i) has %i digits' % (n, len(str(primorial(n)))))
- Output:
First ten primorials: [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870] primorial(1) has 1 digits primorial(10) has 10 digits primorial(100) has 220 digits primorial(1000) has 3393 digits primorial(10000) has 45337 digits primorial(100000) has 563921 digits primorial(1000000) has 6722809 digits
Quackery
eratosthenes and isprime are defined at Sieve of Eratosthenes#Quackery.
[ 0 swap
[ dip 1+
10 /
dup 0 = until ]
drop ] is digits ( n --> n )
[ stack ] is primorials ( --> s )
1299710 eratosthenes
' [ 1 ]
1299710 times
[ i^ isprime if
[ i^ over -1 peek * join ] ]
primorials put
primorials share 10 split drop echo
cr
[] 6 times
[ primorials share
10 i^ ** peek
digits join ]
echo- Output:
[ 1 2 6 30 210 2310 30030 510510 9699690 223092870 ] [ 1 10 220 3393 45337 563921 ]
Racket
We have to reimplement nth-prime and replace it with a memorized version, to make it faster. But we can't memorize primorial because it would use too much memory.
#lang racket
(require (except-in math/number-theory nth-prime))
(define-syntax-rule (define/cache (name arg) body ...)
(begin
(define cache (make-hash))
(define (name arg)
(hash-ref! cache arg (lambda () body ...)))))
(define (num-length n)
;warning: this defines (num-length 0) as 0
(if (zero? n)
0
(add1 (num-length (quotient n 10)))))
(define/cache (nth-prime n)
(if (zero? n)
2
(for/first ([p (in-naturals (add1 (nth-prime (sub1 n))))]
#:when (prime? p))
p)))
(define (primorial n)
(if (zero? n)
1
(* (primorial (sub1 n))
(nth-prime (sub1 n)))))
(displayln
(for/list ([i (in-range 10)])
(primorial i)))
(for ([i (in-range 1 6)])
(printf "Primorial(10^~a) has ~a digits.\n"
i
(num-length (primorial (expt 10 i)))))
- Output:
(1 2 6 30 210 2310 30030 510510 9699690 223092870) Primorial(10^1) has 10 digits. Primorial(10^2) has 220 digits. Primorial(10^3) has 3393 digits. Primorial(10^4) has 45337 digits. Primorial(10^5) has 563921 digits.
Raku
(formerly Perl 6)
Native Raku
With the module Math::Primesieve, this runs in about 1/3 the time vs the built in prime generator.
use Math::Primesieve;
my $sieve = Math::Primesieve.new;
my @primes = $sieve.primes(10_000_000);
sub primorial($n) { [*] @primes[^$n] }
say "First ten primorials: {(primorial $_ for ^10)}";
say "primorial(10^$_) has {primorial(10**$_).chars} digits" for 1..5;
- Output:
First ten primorials: 1 2 6 30 210 2310 30030 510510 9699690 223092870 primorial(10^1) has 10 digits primorial(10^2) has 220 digits primorial(10^3) has 3393 digits primorial(10^4) has 45337 digits primorial(10^5) has 563921 digits
Imported library
For a real speed boost, load the Perl 5 ntheory library.
use Lingua::EN::Numbers;
use ntheory:from<Perl5> <pn_primorial>;
say "First ten primorials: ", ^10 .map( { pn_primorial($_) } ).join: ', ';
for 1..8 {
my $now = now;
printf "primorial(10^%d) has %-11s digits - %s\n", $_,
comma(pn_primorial(10**$_).Str.chars),
"Elapsed seconds: {(now - $now).round: .001}";
}
First ten primorials: 1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870 primorial(10^1) has 10 digits - Elapsed seconds: 0.002 primorial(10^2) has 220 digits - Elapsed seconds: 0.015 primorial(10^3) has 3,393 digits - Elapsed seconds: 0.001 primorial(10^4) has 45,337 digits - Elapsed seconds: 0.005 primorial(10^5) has 563,921 digits - Elapsed seconds: 0.139 primorial(10^6) has 6,722,809 digits - Elapsed seconds: 3.015 primorial(10^7) has 77,919,922 digits - Elapsed seconds: 56.483 primorial(10^8) has 885,105,237 digits - Elapsed seconds: 981.71
REXX
Version 1
/*REXX program computes some primorial numbers for low numbers, and for various 10^n.*/
parse arg N H . /*get optional arguments: N, L, H */
if N=='' | N==',' then N= 10 /*Not specified? Then use the default.*/
if H=='' | H==',' then H= 100000 /* " " " " " " */
numeric digits 600000 /*be able to handle gihugic numbers. */
w= length( commas( digits() ) ) /*W: width of the largest commatized #*/
@.=.; @.0= 1; @.1= 2; @.2= 3; @.3= 5; @.4= 7; @.5= 11; @.6= 13 /*some low primes.*/
s.1= 4; s.2= 9; s.3= 25; s.4= 49; s.5= 121; s.6= 169 /*squared primes. */
#= 6 /*number of primes*/
do j=0 for N /*calculate the first N primorial #s.*/
say right(j, length(N))th(j) " primorial is: " right(commas(primorial(j) ), N+2)
end /*j*/
say
iw= length( commas(H) ) + 2 /*IW: width of largest commatized index*/
p= 1 /*initialize the first multiplier for P*/
do k=1 for H /*process a large range of numbers. */
p= p * prime(k) /*calculate the next primorial number. */
parse var k L 2 '' -1 R /*get the left and rightmost dec digits*/
if R\==0 then iterate /*if right─most decimal digit\==0, skip*/
if L\==1 then iterate /* " left─most " " \==1, " */
if strip(k, , 0)\==1 then iterate /*Not a power of 10? Then skip this K.*/
say right( commas(k), iw)th(k) ' primorial number length in decimal digits is:' ,
right( commas( length(p) ), w)
end /*k*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _; do ?=length(_)-3 to 1 by -3; _=insert(',', _, ?); end; return _
th: parse arg th; return word('th st nd rd', 1+ (th//10)*(th//100%10\==1)*(th//10<4))
/*──────────────────────────────────────────────────────────────────────────────────────*/
primorial: procedure expose @. s. #; parse arg y; != 1 /*obtain the arg Y. */
do p=0 to y; != ! * prime(p) /*calculate product. */
end /*p*/; return ! /*return with the #. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
prime: procedure expose @. s. #; parse arg n; if @.n\==. then return @.n
numeric digits 9 /*limit digs to min.*/
do j=@.#+2 by 2 /*start looking at #*/
if j//2==0 then iterate; if j//3==0 then iterate /*divisible by 2│3 ?*/
parse var j '' -1 _; if _==5 then iterate /*right─most dig≡5? */
if j//7==0 then iterate; if j//11==0 then iterate /*divisible by 7│11?*/
do k=6 while s.k<=j; if j//@.k==0 then iterate j /*divide by primes. */
end /*k*/
#= # + 1; @.#= j; s.#= j * j; return j /*next prime; return*/
end /*j*/
- output when using the default inputs:
0th primorial is: 1
1st primorial is: 2
2nd primorial is: 6
3rd primorial is: 30
4th primorial is: 210
5th primorial is: 2,310
6th primorial is: 30,030
7th primorial is: 510,510
8th primorial is: 9,699,690
9th primorial is: 223,092,870
10th primorial number length in decimal digits is: 10
100th primorial number length in decimal digits is: 220
1,000th primorial number length in decimal digits is: 3,393
10,000th primorial number length in decimal digits is: 45,337
100,000th primorial number length in decimal digits is: 563,921
950 seconds (Regina)
Version 2
Libraries: How to use
Library: Functions
Library: Numbers
Library: Constants
Library: Settings
Library: Abend
Library: Sequences
This task is challenging, since it involves calculations with very big integers (numeric digits) or very big magnitude floats (exponent). Primorial(10^4) already has 45k digits, Primorial(10^5) 564k digits and Primorial(10^6) even almost 7m digits... Let's see. REXX has here in practice no limitations. Part 2 of the task asks for the numbers of digits. For this part, my program first collects all primes up to the 1000000th, using Primes() (a wheeled sieve) with argument Nthprime() (lower bound for the Nth prime). See Numbers and Sequences. Then, I apply 3 methods for finding the number of digits: 1. Multiply all primes in floating point with a precision of say 10 digits and fetch the exponent. 2. Sum all values log10(prime) also in 10 digits and use the outcome. 3. Multiply all primes as integers in the required precision (45k, 563k, ... see above) and fetch the length.
include Settings
say 'Primorial numbers'
parse version version; say version; say
call GetPrimes -9
call ShowFirst10
call GetPrimes -1e6
call UseFloat 1e6
call UseLog 1e6
call UseInt 1e4
exit
GetPrimes:
arg x
call Time('r')
say 'Get primes up to the' Abs(x)'th:'
numeric digits 10
prim. = 0
call Primes(x)
say Format(Time('e'),,3) 'seconds'; say
return
ShowFirst10:
call Time('r')
say 'The first 10 primorials:'
numeric digits 10
say 'Primorial(0) = 1'
a = 1
do i = 1 to 9
a = a*prim.prime.i
say 'Primorial('i') =' a
end
say Format(Time('e'),,3) 'seconds'; say
return
UseFloat:
arg x
call Time('r')
say 'Number of digits by real arithmetic and taking exponent:'
numeric digits 10
a = 1; e = 0
do i = 1 to x
a = a*prim.prime.i
f = Xpon(i)
if f <> e then do
say 'Primorial(10^'f') has' Xpon(a)+1 'digits'
e = f
end
end
say Format(Time('e'),,3) 'seconds'; say
return
UseLog:
arg x
call Time('r')
say 'Number of digits by summing log10:'
numeric digits 10
a = 0; e = 0
do i = 1 to x
a = a+Log(prim.prime.i)
f = Xpon(i)
if f <> e then do
say 'Primorial(10^'f') has' Trunc(a)+1 'digits'
e = f
end
end
say Format(Time('e'),,3) 'seconds'; say
return
UseInt:
arg x
call Time('r')
say 'Number of digits by integer arithmetic and taking length:'
select
when x = 10 then
numeric digits 12
when x = 100 then
numeric digits 230
when x = 1000 then
numeric digits 3400
when x = 10000 then
numeric digits 46000
when x = 100000 then
numeric digits 570000
otherwise
nop
end
a = 1; e = 0
do i = 1 to x
a = a*prim.prime.i
f = Xpon(i)
if f <> e then do
say 'Primorial(10^'f') has' Length(a) 'digits'
e = f
end
end
say Format(Time('e'),,3) 'seconds'; say
return
include Numbers
include Functions
include Sequences
include Constants
include Abend
Notice UseInt gets only to 1e4, while the other two up to 1e6.
- Output:
REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024 Primorial numbers Get primes up to the 9th: 0.000 seconds The first 10 primorials: Primorial(0) = 1 Primorial(1) = 2 Primorial(2) = 6 Primorial(3) = 30 Primorial(4) = 210 Primorial(5) = 2310 Primorial(6) = 30030 Primorial(7) = 510510 Primorial(8) = 9699690 Primorial(9) = 223092870 0.000 seconds Get primes up to the 1000000th: 6.182 seconds Number of digits by real arithmetic and taking exponent: Primorial(10^1) has 10 digits Primorial(10^2) has 220 digits Primorial(10^3) has 3393 digits Primorial(10^4) has 45337 digits Primorial(10^5) has 563921 digits Primorial(10^6) has 6722809 digits 7.310 seconds Number of digits by summing log10: Primorial(10^1) has 10 digits Primorial(10^2) has 220 digits Primorial(10^3) has 3393 digits Primorial(10^4) has 45337 digits Primorial(10^5) has 563921 digits Primorial(10^6) has 6722809 digits 65.473 seconds Number of digits by integer arithmetic and taking length: Primorial(10^1) has 10 digits Primorial(10^2) has 220 digits Primorial(10^3) has 3393 digits Primorial(10^4) has 45337 digits 5.778 seconds
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 Primorial numbers Get primes up to the 9th: 0.000 seconds The first 10 primorials: Same as Regina 0.000 seconds Get primes up to the 1000000th: 12.777 seconds Number of digits by real arithmetic and taking exponent: Same as Regina 2.631 seconds Number of digits by summing log10: Same as Regina 33.613 seconds Number of digits by integer arithmetic and taking length: Same as Regina 1013.082 seconds
These results are interesting. ooRexx is considerable faster with floating point calculations, but much MUCH slower in working with big integers (here up to 45k digits).
I tried Regina with 'call UseInt 1e5', requiring the handling of integers with up to 560k digits.
REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024 Primorial numbers Number of digits by integer arithmetic and taking length: Primorial(10^1) has 10 digits Primorial(10^2) has 220 digits Primorial(10^3) has 3393 digits Primorial(10^4) has 45337 digits Primorial(10^5) has 563921 digits 960.408 seconds
Hardly.... but it came through! Procedure UseInt might probable be improved by dynamically increasing numeric digits ongoing the calculations, so that the lower primorials are calculated with only the needed precision in stead of 560k digits.
Ring
# Project: Primorial numbers
load "bignumber.ring"
decimals(0)
num = 0
prim = 0
limit = 10000000
see "working..." + nl
see "wait for done..." + nl
while num < 100001
prim = prim + 1
prime = []
primorial(prim)
end
see "done..." + nl
func primorial(pr)
n = 1
n2 = 0
flag = 1
while flag = 1 and n < limit
nr = isPrime(n)
if n=1
nr=1
ok
if nr=1
n2 = n2 + 1
add(prime,n)
ok
if n2=pr
flag=0
num = num + 1
ok
n = n + 1
end
pro = 1
str = ""
for n=1 to len(prime)
pro = FuncMultiply("" + pro,"" + prime[n])
str = str + prime[n] + "*"
next
str = left(str,len(str)-1)
if pr < 11
see "primorial(" + string(pr-1) + ") : " + pro + nl
ok
if pr = 11
see "primorial(" + string(pr-1) + ") " + "has " + len(pro) + " digits"+ nl
but pr = 101
see "primorial(" + string(pr-1) + ") " + "has " + len(pro) + " digits"+ nl
but pr = 1001
see "primorial(" + string(pr-1) + ") " + "has " + len(pro) + " digits"+ nl
but pr = 10001
see "primorial(" + string(pr-1) + ") " + "has " + len(pro) + " digits"+ nl
but pr = 100001
see "primorial(" + string(pr-1) + ") " + "has " + len(pro) + " digits"+ nl
okworking... wait for done... primorial(0) = 1 primorial(1) = 2 primorial(2) = 6 primorial(3) = 30 primorial(4) = 210 primorial(5) = 2310 primorial(6) = 30030 primorial(7) = 510510 primorial(8) = 9699690 primorial(9) = 223092870 primorial(10) has 10 digits primorial(100) has 220 digits primorial(1000) has 3393 digits primorial(10000) has 45337 digits primorial(100000) has 563921 digits done...
Ruby
require 'prime'
def primorial_number(n)
pgen = Prime.each
(1..n).inject(1){|p,_| p*pgen.next}
end
puts "First ten primorials: #{(0..9).map{|n| primorial_number(n)}}"
(1..5).each do |n|
puts "primorial(10**#{n}) has #{primorial_number(10**n).to_s.size} digits"
end
- Output:
First ten primorials: [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870] primorial(10**1) has 10 digits primorial(10**2) has 220 digits primorial(10**3) has 3393 digits primorial(10**4) has 45337 digits primorial(10**5) has 563921 digits
Rust
extern crate primal;
extern crate rayon;
extern crate rug;
use rayon::prelude::*;
use rug::Integer;
fn partial(p1 : usize, p2 : usize) -> String {
let mut aux = Integer::from(1);
let (_, hi) = primal::estimate_nth_prime(p2 as u64);
let sieve = primal::Sieve::new(hi as usize);
let prime1 = sieve.nth_prime(p1);
let prime2 = sieve.nth_prime(p2);
for i in sieve.primes_from(prime1).take_while(|i| *i <= prime2) {
aux = Integer::from(aux * i as u32);
}
aux.to_string_radix(10)
}
fn main() {
let mut j1 = Integer::new();
for k in [2,3,5,7,11,13,17,19,23,29].iter() {
j1.assign_primorial(*k);
println!("Primorial : {}", j1);
}
println!("Digits of primorial 10 : {}", partial(1, 10).chars().fold(0, |n, _| n + 1));
println!("Digits of primorial 100 : {}", partial(1, 100).chars().fold(0, |n, _| n + 1));
println!("Digits of primorial 1_000 : {}", partial(1, 1_000).chars().fold(0, |n, _| n + 1));
println!("Digits of primorial 10_000 : {}", partial(1, 10_000).chars().fold(0, |n, _| n + 1));
println!("Digits of primorial 100_000 : {}", partial(1, 100_000).chars().fold(0, |n, _| n + 1));
let mut auxi = Integer::from(1);
let ranges = vec![[1, 300_000], [300_001, 550_000], [550_001, 800_000], [800_001, 1_000_000]];
let v = ranges.par_iter().map(|value| partial(value[0], value[1])).collect::<Vec<_>>();
for i in v.iter() {
auxi =Integer::from(&auxi * i.parse::<Integer>().unwrap());
}
let result = auxi.to_string_radix(10).chars().fold(0, |n, _| n+1);
println!("Digits of primorial 1_000_000 : {}",result);
}
using Intel(R) Core(TM) i7-5500U CPU @ 2.40GHz
- Output:
Primorial : 2 Primorial : 6 Primorial : 30 Primorial : 210 Primorial : 2310 Primorial : 30030 Primorial : 510510 Primorial : 9699690 Primorial : 223092870 Primorial : 6469693230 Digits of primorial 10 : 10 Digits of primorial 100 : 220 Digits of primorial 1_000 : 3393 Digits of primorial 10_000 : 45337 Digits of primorial 100_000 : 563921 Digits of primorial 1_000_000 : 6722809 real 0m19.448s user 1m2.306s sys 0m0.054s
Scala
This example uses Spire's SafeLong for arbitrarily large integers and parallel vectors to accelerate bulk operations on large collections. The primorial is calculated by simply building a list of primes and taking the product. The prime generator here is relatively inefficient, but as it isn't the focus of this task the priority is brevity while retaining acceptable performance.
import spire.math.SafeLong
import spire.implicits._
import scala.collection.parallel.immutable.ParVector
object Primorial {
def main(args: Array[String]): Unit = {
println(
s"""|First 10 Primorials:
|${LazyList.range(0, 10).map(n => f"$n: ${primorial(n).toBigInt}%,d").mkString("\n")}
|
|Lengths of Primorials:
|${LazyList.range(1, 7).map(math.pow(10, _).toInt).map(i => f"$i%,d: ${primorial(i).toString.length}%,d").mkString("\n")}
|""".stripMargin)
}
def primorial(num: Int): SafeLong = if(num == 0) 1 else primesSL.take(num).to(ParVector).reduce(_*_)
lazy val primesSL: Vector[SafeLong] = 2 +: ParVector.range(3, 20000000, 2).filter(n => !Iterator.range(3, math.sqrt(n).toInt + 1, 2).exists(n%_ == 0)).toVector.sorted.map(SafeLong(_))
}
- Output:
First 10 Primorials: 0: 1 1: 2 2: 6 3: 30 4: 210 5: 2,310 6: 30,030 7: 510,510 8: 9,699,690 9: 223,092,870 Lengths of Primorials: 10: 10 100: 220 1,000: 3,393 10,000: 45,337 100,000: 563,921 1,000,000: 6,722,809
Sidef
say (
'First ten primorials: ',
{|i| pn_primorial(i) }.map(^10).join(', ')
)
{ |i|
say ("primorial(10^#{i}) has " + pn_primorial(10**i).len + ' digits')
} << 1..6
- Output:
First ten primorials: 1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870 primorial(10^1) has 10 digits primorial(10^2) has 220 digits primorial(10^3) has 3393 digits primorial(10^4) has 45337 digits primorial(10^5) has 563921 digits primorial(10^6) has 6722809 digits
Wren
Wren CLI
It takes very little time to sieve for the first 100,000 primes (< 0.2 seconds) but, despite using the 'binary splitting' algorithm and a trick of my own, multiplying them all together is very slow compared to the GMP-based solutions taking more than 4 minutes on my machine.
import "./big" for BigInt
import "./math" for Int
import "./fmt" for Fmt
var vecprod = Fn.new { |primes|
var le = primes.count
if (le == 0) return BigInt.one
var s = List.filled(le, null)
for (i in 0...le) s[i] = BigInt.new(primes[i])
while (le > 1) {
var c = (le/2).floor
for(i in 0...c) s[i] = s[i] * s[le-i-1]
if (le & 1 == 1) c = c + 1
le = c
}
return s[0]
}
var primes = Int.primeSieve(1.3e6) // enough to generate first 100,000 primes
var prod = 1
System.print("The first ten primorial numbers are:")
for (i in 0..9) {
System.print("%(i): %(prod)")
prod = prod * primes[i]
}
System.print("\nThe following primorials have the lengths shown:")
// first multiply the first 100,000 primes together in pairs to reduce BigInt conversions needed
var primes2 = List.filled(50000, 0)
for (i in 0...50000) primes2[i] = primes[2*i] * primes[2*i+1]
for (i in [10, 100, 1000, 10000, 100000]) {
Fmt.print("$6d: $d", i, vecprod.call(primes2[0...i/2]).toString.count)
}
- Output:
The first ten primorial numbers are:
0: 1
1: 2
2: 6
3: 30
4: 210
5: 2310
6: 30030
7: 510510
8: 9699690
9: 223092870
The following primorials have the lengths shown:
10: 10
100: 220
1000: 3393
10000: 45337
100000: 563921
Embedded
Far quicker, of course than the CLI version, completing in around 2.0 seconds.
import "./math" for Int
import "./gmp" for Mpz
import "./fmt" for Fmt
var limit = 16 * 1e6 // more than enough to find first million primes
var primes = Int.primeSieve(limit-1)
primes.insert(0, 1)
System.print("The first ten primorial numbers are:")
var z = Mpz.new()
for (i in 0..9) System.print("%(i): %(z.primorial(primes[i]))")
System.print("\nThe following primorials have the lengths shown:")
for (i in [1e1, 1e2, 1e3, 1e4, 1e5, 1e6]) {
Fmt.print("$7d: $d", i, z.primorial(primes[i]).digitsInBase(10))
}
- Output:
The first ten primorial numbers are:
0: 1
1: 2
2: 6
3: 30
4: 210
5: 2310
6: 30030
7: 510510
8: 9699690
9: 223092870
The following primorials have the lengths shown:
10: 10
100: 220
1000: 3393
10000: 45337
100000: 563921
1000000: 6722809
zkl
Using Extensible prime generator#zkl and the GMP big int library.
sieve:=Import("sieve.zkl",False,False,False).postponed_sieve;
primes:=Utils.Generator(sieve).walk(0d10); // first 10 primes
foreach n in (10)
{ primes[0,n].reduce('*,1):println("primorial(%d)=%d".fmt(n,_)); }
var [const] BN=Import("zklBigNum");
primes:=Utils.Generator(sieve).walk(0d1_000_000);
foreach n in ([1..6]){ n=(10).pow(n);
primes[0,n].pump(BN(1).mul)
:println("primorial(%,d)=%,d digits".fmt(n,_.numDigits));
}Big int multiplication is done in place to minimize garbage. Also, subsets of read only lists (which the list of primes is) are not copies (ie they are a small header that points into the original list).
Strangely, it takes much longer to do one million big int multiplications than to generate one million primes, so it would be a good idea to "generate a prime, multiply, print at intervals" but I'm too lazy.
- Output:
primorial(0)=1 primorial(1)=2 primorial(2)=6 primorial(3)=30 primorial(4)=210 primorial(5)=2310 primorial(6)=30030 primorial(7)=510510 primorial(8)=9699690 primorial(9)=223092870 primorial(10)=10 digits primorial(100)=220 digits primorial(1,000)=3,393 digits primorial(10,000)=45,338 digits primorial(100,000)=563,921 digits primorial(1,000,000)=6,722,809 digits
- Programming Tasks
- Prime Numbers
- 11l
- Arturo
- C
- GMP
- C sharp
- C++
- Primesieve
- Clojure
- CLU
- Common Lisp
- D
- Delphi
- EasyLang
- Elixir
- F Sharp
- Factor
- Fortran
- FreeBASIC
- Frink
- Go
- GMP(Go wrapper)
- Haskell
- J
- Java
- JavaScript
- Jq
- Julia
- Kotlin
- Lingo
- Maple
- Mathematica
- Wolfram Language
- Nickle
- Nim
- Bignum
- PARI/GP
- Pascal
- Perl
- Ntheory
- Phix
- Phix/mpfr
- PicoLisp
- Python
- Quackery
- Racket
- Raku
- REXX
- Ring
- Ruby
- Rust
- Scala
- Sidef
- Wren
- Wren-big
- Wren-math
- Wren-fmt
- Wren-gmp
- Zkl