Skip to content

Latest commit

 

History

History
182 lines (141 loc) · 12.7 KB

factorial.md

File metadata and controls

182 lines (141 loc) · 12.7 KB

Factorial

Function signature

The factorial function is defined as:

bigint_t factorial(size_t n);

This calculates, obviously, $n!$.

The in-parameter n is limited to type size_t because:

  • Calculating a factorial anything near $10^9!$ is very long in practice, let alone values closer to $10^{19}!$ which is still far from what size_t can store on 64-bit architectures ($\approx 1.845 \times 10^{19}$).
  • Since size_t is linked to how much memory a system holds, the space required to store the result of factorial bigger than size_t integers exceeds the capabilities. In practice, that limit is hit long before reaching the maximum value a size_t can express.

Naive implementation

The naive implementation for factorial calculation is along the lines of:

bigint::bigint_t bigint::factorial(size_t n)
{
  bigint::bigint_t result = 1;
  for (size_t i = 2; i <= n; ++i)
    result *= i;
  return result;
}

This approach suffers from several downsides:

  • It requires n - 1 multiplications (from 2 to n) and importantly, some of these products are done more than once; with the algorithmic complexity of the multiplication, this is something to avoid.
    For instance, we would want $3$, $6$, $12$, $24$ (i.e. $3$'s multiplied by a power of $2$) to be processed as $64 \times (3^2)^2$: This requires 2 multiplications (The multiplication by $64$ is done through bit shifting) instead of 3.
  • The operands passed in the successive multiplications are very different in size.
    Having one of the operand consiting of 1 or 2 limbs prevents from using the more efficient multiplication algorithms.
  • As was mentioned above, we can make use of bit shifting operations when multiplying even numbers.
    In fact, it is possible to account for all the $2$'s encountered between $1$ and $n$ with a single shift.

The library's approach

Principle

Instead of the naive implementation, the implemented algorithm works "backward", with only odd numbers ($3$, $5$, $7$, $9$ ...). All the even numbers will be managed at once, by shifting bits left.
For the sake of clarity, we will illustrate each step with $n = 20$. The table below represents the factors; when used, they get crossed out.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Algorithm

  1. Initialize the variables $\text{result} \leftarrow 1$, $\text{product} \leftarrow 1$ and $\text{c} \leftarrow 1$.

  2. Determine $h$ so that $2^h \leq n &lt; 2^{h+1}$.
    Initialize $\text{shift} \leftarrow h \times (h+1) / 2$

  3. Determine $p$ so that $2^p \times 3 \leq n$. It is possible that $2^p \times 5 \leq n$.
    Example: with $n = 20$: $2^3 \times 3 = 24 \gt 20$ and $2^2 \times 3 = 12 \leq 20$. This means $p = 2$.

  4. Find the highest odd integer $i$ that verifies $2^p \times i \leq n$. With $n = 20$ and $p = 2$, that number is $5$ ($2^2 \times 5 \leq 20$).
    Do:
    $\text{shift} \leftarrow \text{shift} + p \times (i-1)/2$
    $\displaystyle \text{product} \leftarrow \text{product} \times \prod_{\begin{array}\text{ } k = c+2 \\ k \text{ is odd} \end{array}}^i(k)\text{ }$
    $c \leftarrow i$.

  5. Do $\text{result} \leftarrow \text{result} \times \text{product}$

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  1. If $p &gt; 0$, do $p \leftarrow p - 1$ then go back to step 4.
  • At $p = 1$, the numbers to add to product are $7$ and $9$ and we consume all the factors that are multiples of $2$ without being multiples of $4$.
    This includes $6$ and $10$ which were already covered by $\text{product}$ for no computing cost at all.
    First: $\text{product} \leftarrow \text{product} \times (7 \times 9)$
    Then $\text{result} \leftarrow \text{result} \times \text{product}$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  • At $p = 0$, the numbers are $11$, $13$, $15$, $17$ and $19$ and we consume all the factors that are not multiples of $2$.
    Like above, $3$, $5$, $7$, $9$ were already included in $\text{product}$, thus are included again. First: $\text{product} \leftarrow \text{product} \times (11 \times 13 \times 15 \times 17 \times 19)$
    Then $\text{result} \leftarrow \text{result} \times \text{product}$
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  1. All that remains is shift the bits of $\text{result}$ by $\text{shift}$ to account for all the even numbers.
    Doing so at the last step saves some execution time as all the previous multiplications did not have to browse through the additional limbs that shift creates.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Time savings

Compared to the naive approach and for big enough factorials:

  • The result is attained with half the number of multiplications required by the naive approach.
    The computational complexity is not improved by that but when every multiplication takes a measurable time, cutting their number by half is an important improvement.
  • It makes better use of the optimized multiplication algorithms, as both $\text{product}$ and $\text{result}$ grow.
    This part of the improvement results in better computational complexity: the algorithm repeats $n$ times a step with an overall better $\text{O}(\text{M}(n))$.

A complete example

With $n = 60$:

$\text{product} \leftarrow 1$
$\text{result} \leftarrow 1$

  1. First loop:
    $p \leftarrow 4$ ($2^p = 16 $)
    $\text{product} \leftarrow \text{product} \times 3 = 3$
    $\text{result} \leftarrow \text{result} \times \text{product} = 3$
    1 value, $48$, gets used up.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
  1. Second loop:
    $p \leftarrow 3$ ($2^p = 8 $)
    $\text{product} \leftarrow \text{product} \times (5 \times 7) = 105$
    $\text{result} \leftarrow \text{result} \times \text{product} = 315$
    3 values, $24$, $40$ and $56$, get used up, with 3 multiplications done.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
  1. Third loop:
    $p \leftarrow 2$ ($2^p = 4 $)
    $\text{product} \leftarrow \text{product} \times (9 \times 11 \times 13 \times 15) = 2{,}027{,}025$
    $\text{result} \leftarrow \text{result} \times \text{product} = 638{,}512{,}875$
    7 values, $12$, $20$, $28$, $36$, $44$, $52$ and $60$, get used up, with 5 multiplications done.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
  1. Fourth loop:
    $p \leftarrow 1$ ($2^p = 2$)
    $\text{product} \leftarrow \text{product} \times (17 \times 19 \times 21 \times 23 \times 25 \times 27 \times 29) = 6{,}190{,}283{,}353{,}629{,}375 $
    $\text{result} \leftarrow \text{result} \times \text{product} = 3{,}952{,}575{,}621{,}190{,}533{,}915{,}703{,}125$
    14 values, $6$, $10$, $14$, $18$, $22$, $26$, $30$, $34$, $38$, $42$, $46$, $50$, $54$ and $58$, get used up, with 8 multiplications done.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
  1. Fifth loop:
    $p \leftarrow 0$ ($2^p = 1$)
    $\text{product} \leftarrow \text{product} \times (31 \times 33 \times 35 \times 37 \times \text{...} \times 57 \times 59) = 29{,}215{,}606{,}371{,}473{,}169{,}285{,}018{,}060{,}091{,}249{,}259{,}296{,}875$
    $\text{result} \leftarrow \text{result} \times \text{product} = 115{,}476{,}893{,}502{,}183{,}682{,}653{,}166{,}335{,}352{,}659{,}171{,}719{,}555{,}028{,}600{,}718{,}376{,}458{,}740{,}234{,}375$
    29 values, $3$, $5$, $7$, $9$, $11$, $13$, $15$, $17$, $19$, ..., $51$, $53$, $55$, $57$ and $59$, get used up, with 16 multiplications done.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
  1. Finalizing:
    The powers of 2 under $n=60$ account for 15 bits to shift ($2^0$, $2^1$, $2^2$, $2^3$, $2^4$, $2^5$).
    During the calculation, we consumed:
    • When $p$ was $4$: $1$ value.
    • When $p$ was $3$: $3$ values.
    • When $p$ was $2$: $7$ values.
    • When $p$ was $1$: $14$ values.

 This makes a total of $15+4 \times 1+3 \times 3 +2 \times 7 + 14 = 56$ bits. Finally, $\text{result} \leftarrow \text{result} \times 2^{56}$ by shifting bits.

$$ \begin{flalign} & 60! = 8{,}320{,}987{,}112{,}741{,}390{,}144{,}276{,}341{,}183{,}223{,}364{,}380{,}754{,}172{,}606{,}361{,}245{,}952{,}449{,}277{,}696{,}409{,}600{,}000{,}000{,}000{,}000 & \end{flalign} $$

$60$ being relatively small, the number of multiplications used for the calculation ($34$) is still significantly higher than the $n/2$ ratio that the algorithm approaches. We can see each successive loop gets more efficient than the previous one, though.