How fast one can calculate the compositorial

computational complexityelementary-number-theory

Definition. The compositorial of a composite number $n$ is the product of all composite numbers up to and including $n$.

Question. How fast one can calculate the compositorial $\displaystyle\prod_{i=1}^n C_i$ , where $n=20000000$ ?

PARI/GP computation gives the following results:

$$\begin{array}{|c|c|} n & time \\
\hline
100000&1,007\ \rm ms\\
200000&4,692\ \rm ms\\
300000&11,100\ \rm ms\\
400000&19,858\ \rm ms\\
500000&32,074\ \rm ms\\
\end{array}$$

Is it possible to estimate time needed for $n=20000000$ ?

EDIT

Using built-in PARI/GP function polinterpolate (see here) , I got estimated time: $93454273,148\ \rm ms$ , nearly $3$ years.

Best Answer

With a decent algorithm and compiler, computing that number is a matter of minutes (I'm not sure if optimisation can take it to be a matter of seconds, but some tweaking can certainly make it faster than it is now). The main point to make such computations fast is to reduce the number of multiplications involving at least one large factor. Thus, instead of sequentially multiplying all the composites to an accumulator it is much faster to arrange the multiplications in a tree, so that one has one multiplication with both factors about half the size of the final result, two multiplications where both factors have about a quarter the final size, and so on.

Haskell programme (note: My GHC is old, some changes may be necessary to make it compile with newer GHCs)

module Compositorial where

import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base

-- multiply all composite numbers <= n
-- Use a sieve for a fast compositeness test, and multiply in a tree-like fashion
compositorial :: Int -> Integer
compositorial n = treeProd (sieve n) 4 n

-- Sieve of Eratosthenes, basic
sieve :: Int -> UArray Int Bool
sieve end = runSTUArray (do
    ar <- newArray (0,end) False
    let mark step idx
            | idx > end = return ()
            | otherwise = unsafeWrite ar idx True >> mark step (idx+step)
        sift p
            | p*p > end = return ar
            | otherwise = do
                c <- unsafeRead ar p
                if c then return () else mark (2*p) (p*p)
                sift (p+2)
    mark 2 4
    sift 3)

-- The meat of it, tree product
-- For larger ranges, split roughly in half and recur,
-- for short ranges multiply directly
treeProd :: UArray Int Bool -> Int -> Int -> Integer
treeProd ar low high = go low high
  where
    go lo hi
        | lo + 4 < hi = let m = (hi + lo) `quot` 2 in (go lo m) * (go (m+1) hi)
        | otherwise   = product [fromIntegral n | n <- [lo .. hi], ar `unsafeAt` n]

This, complied with -O2 -fllvm, computed the compositorial of $21350730$ (since $\pi(21350730) = 1350729$, $21350730$ is the $20000000^{\text{th}}$ composite number), a number with $137945607$ digits, and printed it to a file, in

Total   time  114.71s  (129.96s elapsed)

just over two minutes.

Related Question