88
\chapter*{Prefaces to the Draft Edition}
89
I started this text in April 2003 to complement my LibTomMath library. That is, explain how to implement the functions
90
contained in LibTomMath. The goal is to have a textbook that any Computer Science student can use when implementing their
91
own multiple precision arithmetic. The plan I wanted to follow was flesh out all the
92
ideas and concepts I had floating around in my head and then work on it afterwards refining a little bit at a time. Chance
93
would have it that I ended up with my summer off from Algonquin College and I was given four months solid to work on the
96
Choosing to not waste any time I dove right into the project even before my spring semester was finished. I wrote a bit
97
off and on at first. The moment my exams were finished I jumped into long 12 to 16 hour days. The result after only
98
a couple of months was a ten chapter, three hundred page draft that I quickly had distributed to anyone who wanted
99
to read it. I had Jean-Luc Cooke print copies for me and I brought them to Crypto'03 in Santa Barbara. So far I have
100
managed to grab a certain level of attention having people from around the world ask me for copies of the text was certain
103
Now we are past December 2003. By this time I had pictured that I would have at least finished my second draft of the text.
104
Currently I am far off from this goal. I've done partial re-writes of chapters one, two and three but they are not even
105
finished yet. I haven't given up on the project, only had some setbacks. First O'Reilly declined to publish the text then
106
Addison-Wesley and Greg is tried another which I don't know the name of. However, at this point I want to focus my energy
107
onto finishing the book not securing a contract.
109
So why am I writing this text? It seems like a lot of work right? Most certainly it is a lot of work writing a textbook.
110
Even the simplest introductory material has to be lined with references and figures. A lot of the text has to be re-written
111
from point form to prose form to ensure an easier read. Why am I doing all this work for free then? Simple. My philosophy
112
is quite simply ``Open Source. Open Academia. Open Minds'' which means that to achieve a goal of open minds, that is,
113
people willing to accept new ideas and explore the unknown you have to make available material they can access freely
116
I've been writing free software since I was about sixteen but only recently have I hit upon software that people have come
117
to depend upon. I started LibTomCrypt in December 2001 and now several major companies use it as integral portions of their
118
software. Several educational institutions use it as a matter of course and many freelance developers use it as
119
part of their projects. To further my contributions I started the LibTomMath project in December 2002 aimed at providing
120
multiple precision arithmetic routines that students could learn from. That is write routines that are not only easy
121
to understand and follow but provide quite impressive performance considering they are all in standard portable ISO C.
123
The second leg of my philosophy is ``Open Academia'' which is where this textbook comes in. In the end, when all is
124
said and done the text will be useable by educational institutions as a reference on multiple precision arithmetic.
126
At this time I feel I should share a little information about myself. The most common question I was asked at
127
Crypto'03, perhaps just out of professional courtesy, was which school I either taught at or attended. The unfortunate
128
truth is that I neither teach at or attend a school of academic reputation. I'm currently at Algonquin College which
129
is what I'd like to call ``somewhat academic but mostly vocational'' college. In otherwords, job training.
131
I'm a 21 year old computer science student mostly self-taught in the areas I am aware of (which includes a half-dozen
132
computer science fields, a few fields of mathematics and some English). I look forward to teaching someday but I am
133
still far off from that goal.
135
Now it would be improper for me to not introduce the rest of the texts co-authors. While they are only contributing
136
corrections and editorial feedback their support has been tremendously helpful in presenting the concepts laid out
137
in the text so far. Greg has always been there for me. He has tracked my LibTom projects since their inception and even
138
sent cheques to help pay tuition from time to time. His background has provided a wonderful source to bounce ideas off
139
of and improve the quality of my writing. Mads is another fellow who has just ``been there''. I don't even recall what
140
his interest in the LibTom projects is but I'm definitely glad he has been around. His ability to catch logical errors
141
in my written English have saved me on several occasions to say the least.
143
What to expect next? Well this is still a rough draft. I've only had the chance to update a few chapters. However, I've
144
been getting the feeling that people are starting to use my text and I owe them some updated material. My current tenative
145
plan is to edit one chapter every two weeks starting January 4th. It seems insane but my lower course load at college
146
should provide ample time. By Crypto'04 I plan to have a 2nd draft of the text polished and ready to hand out to as many
147
people who will take it.
89
When I tell people about my LibTom projects and that I release them as public domain they are often puzzled.
90
They ask why I did it and especially why I continue to work on them for free. The best I can explain it is ``Because I can.''
91
Which seems odd and perhaps too terse for adult conversation. I often qualify it with ``I am able, I am willing.'' which
92
perhaps explains it better. I am the first to admit there is not anything that special with what I have done. Perhaps
93
others can see that too and then we would have a society to be proud of. My LibTom projects are what I am doing to give
94
back to society in the form of tools and knowledge that can help others in their endeavours.
96
I started writing this book because it was the most logical task to further my goal of open academia. The LibTomMath source
97
code itself was written to be easy to follow and learn from. There are times, however, where pure C source code does not
98
explain the algorithms properly. Hence this book. The book literally starts with the foundation of the library and works
99
itself outwards to the more complicated algorithms. The use of both pseudo--code and verbatim source code provides a duality
100
of ``theory'' and ``practice'' that the computer science students of the world shall appreciate. I never deviate too far
101
from relatively straightforward algebra and I hope that this book can be a valuable learning asset.
103
This book and indeed much of the LibTom projects would not exist in their current form if it was not for a plethora
104
of kind people donating their time, resources and kind words to help support my work. Writing a text of significant
105
length (along with the source code) is a tiresome and lengthy process. Currently the LibTom project is four years old,
106
comprises of literally thousands of users and over 100,000 lines of source code, TeX and other material. People like Mads and Greg
107
were there at the beginning to encourage me to work well. It is amazing how timely validation from others can boost morale to
108
continue the project. Definitely my parents were there for me by providing room and board during the many months of work in 2003.
110
To my many friends whom I have met through the years I thank you for the good times and the words of encouragement. I hope I
111
honour your kind gestures with this project.
113
Open Source. Open Academia. Open Minds.
149
115
\begin{flushright} Tom St Denis \end{flushright}
1398
1371
EXAM,bn_mp_set.c
1400
Line @21,mp_zero@ calls mp\_zero() to clear the mp\_int and reset the sign. Line @22,MP_MASK@ copies the digit
1401
into the least significant location. Note the usage of a new constant \textbf{MP\_MASK}. This constant is used to quickly
1402
reduce an integer modulo $\beta$. Since $\beta$ is of the form $2^k$ for any suitable $k$ it suffices to perform a binary AND with
1403
$MP\_MASK = 2^k - 1$ to perform the reduction. Finally line @23,a->used@ will set the \textbf{used} member with respect to the
1404
digit actually set. This function will always make the integer positive.
1373
First we zero (line @21,mp_zero@) the mp\_int to make sure that the other members are initialized for a
1374
small positive constant. mp\_zero() ensures that the \textbf{sign} is positive and the \textbf{used} count
1375
is zero. Next we set the digit and reduce it modulo $\beta$ (line @22,MP_MASK@). After this step we have to
1376
check if the resulting digit is zero or not. If it is not then we set the \textbf{used} count to one, otherwise
1379
We can quickly reduce modulo $\beta$ since it is of the form $2^k$ and a quick binary AND operation with
1380
$2^k - 1$ will perform the same operation.
1406
1382
One important limitation of this function is that it will only set one digit. The size of a digit is not fixed, meaning source that uses
1407
1383
this function should take that into account. Only trivially small constants can be set using this function.
1665
1643
EXAM,bn_s_mp_add.c
1667
Lines @27,if@ to @35,}@ perform the initial sorting of the inputs and determine the $min$ and $max$ variables. Note that $x$ is a pointer to a
1668
mp\_int assigned to the largest input, in effect it is a local alias. Lines @37,init@ to @42,}@ ensure that the destination is grown to
1669
accomodate the result of the addition.
1645
We first sort (lines @27,if@ to @35,}@) the inputs based on magnitude and determine the $min$ and $max$ variables.
1646
Note that $x$ is a pointer to an mp\_int assigned to the largest input, in effect it is a local alias. Next we
1647
grow the destination (@37,init@ to @42,}@) ensure that it can accomodate the result of the addition.
1671
1649
Similar to the implementation of mp\_copy this function uses the braced code and local aliases coding style. The three aliases that are on
1672
1650
lines @56,tmpa@, @59,tmpb@ and @62,tmpc@ represent the two inputs and destination variables respectively. These aliases are used to ensure the
1673
1651
compiler does not have to dereference $a$, $b$ or $c$ (respectively) to access the digits of the respective mp\_int.
1675
The initial carry $u$ is cleared on line @65,u = 0@, note that $u$ is of type mp\_digit which ensures type compatibility within the
1676
implementation. The initial addition loop begins on line @66,for@ and ends on line @75,}@. Similarly the conditional addition loop
1677
begins on line @81,for@ and ends on line @90,}@. The addition is finished with the final carry being stored in $tmpc$ on line @94,tmpc++@.
1678
Note the ``++'' operator on the same line. After line @94,tmpc++@ $tmpc$ will point to the $c.used$'th digit of the mp\_int $c$. This is useful
1679
for the next loop on lines @97,for@ to @99,}@ which set any old upper digits to zero.
1653
The initial carry $u$ will be cleared (line @65,u = 0@), note that $u$ is of type mp\_digit which ensures type
1654
compatibility within the implementation. The initial addition (line @66,for@ to @75,}@) adds digits from
1655
both inputs until the smallest input runs out of digits. Similarly the conditional addition loop
1656
(line @81,for@ to @90,}@) adds the remaining digits from the larger of the two inputs. The addition is finished
1657
with the final carry being stored in $tmpc$ (line @94,tmpc++@). Note the ``++'' operator within the same expression.
1658
After line @94,tmpc++@, $tmpc$ will point to the $c.used$'th digit of the mp\_int $c$. This is useful
1659
for the next loop (line @97,for@ to @99,}@) which set any old upper digits to zero.
1681
1661
\subsection{Low Level Subtraction}
1682
1662
The low level unsigned subtraction algorithm is very similar to the low level unsigned addition algorithm. The principle difference is that the
1760
1740
EXAM,bn_s_mp_sub.c
1762
Line @24,min@ and @25,max@ perform the initial hardcoded sorting of the inputs. In reality the $min$ and $max$ variables are only aliases and are only
1763
used to make the source code easier to read. Again the pointer alias optimization is used within this algorithm. Lines @42,tmpa@, @43,tmpb@ and @44,tmpc@ initialize the aliases for
1764
$a$, $b$ and $c$ respectively.
1766
The first subtraction loop occurs on lines @47,u = 0@ through @61,}@. The theory behind the subtraction loop is exactly the same as that for
1767
the addition loop. As remarked earlier there is an implementation reason for using the ``awkward'' method of extracting the carry
1768
(\textit{see line @57, >>@}). The traditional method for extracting the carry would be to shift by $lg(\beta)$ positions and logically AND
1769
the least significant bit. The AND operation is required because all of the bits above the $\lg(\beta)$'th bit will be set to one after a carry
1770
occurs from subtraction. This carry extraction requires two relatively cheap operations to extract the carry. The other method is to simply
1771
shift the most significant bit to the least significant bit thus extracting the carry with a single cheap operation. This optimization only works on
1772
twos compliment machines which is a safe assumption to make.
1774
If $a$ has a larger magnitude than $b$ an additional loop (\textit{see lines @64,for@ through @73,}@}) is required to propagate the carry through
1775
$a$ and copy the result to $c$.
1742
Like low level addition we ``sort'' the inputs. Except in this case the sorting is hardcoded
1743
(lines @24,min@ and @25,max@). In reality the $min$ and $max$ variables are only aliases and are only
1744
used to make the source code easier to read. Again the pointer alias optimization is used
1745
within this algorithm. The aliases $tmpa$, $tmpb$ and $tmpc$ are initialized
1746
(lines @42,tmpa@, @43,tmpb@ and @44,tmpc@) for $a$, $b$ and $c$ respectively.
1748
The first subtraction loop (lines @47,u = 0@ through @61,}@) subtract digits from both inputs until the smaller of
1749
the two inputs has been exhausted. As remarked earlier there is an implementation reason for using the ``awkward''
1750
method of extracting the carry (line @57, >>@). The traditional method for extracting the carry would be to shift
1751
by $lg(\beta)$ positions and logically AND the least significant bit. The AND operation is required because all of
1752
the bits above the $\lg(\beta)$'th bit will be set to one after a carry occurs from subtraction. This carry
1753
extraction requires two relatively cheap operations to extract the carry. The other method is to simply shift the
1754
most significant bit to the least significant bit thus extracting the carry with a single cheap operation. This
1755
optimization only works on twos compliment machines which is a safe assumption to make.
1757
If $a$ has a larger magnitude than $b$ an additional loop (lines @64,for@ through @73,}@) is required to propagate
1758
the carry through $a$ and copy the result to $c$.
1777
1760
\subsection{High Level Addition}
1778
1761
Now that both lower level addition and subtraction algorithms have been established an effective high level signed addition algorithm can be
2465
2464
EXAM,bn_s_mp_mul_digs.c
2467
Lines @31,if@ to @35,}@ determine if the Comba method can be used first. The conditions for using the Comba routine are that min$(a.used, b.used) < \delta$ and
2468
the number of digits of output is less than \textbf{MP\_WARRAY}. This new constant is used to control
2469
the stack usage in the Comba routines. By default it is set to $\delta$ but can be reduced when memory is at a premium.
2471
Of particular importance is the calculation of the $ix+iy$'th column on lines @64,mp_word@, @65,mp_word@ and @66,mp_word@. Note how all of the
2472
variables are cast to the type \textbf{mp\_word}, which is also the type of variable $\hat r$. That is to ensure that double precision operations
2473
are used instead of single precision. The multiplication on line @65,) * (@ makes use of a specific GCC optimizer behaviour. On the outset it looks like
2474
the compiler will have to use a double precision multiplication to produce the result required. Such an operation would be horribly slow on most
2475
processors and drag this to a crawl. However, GCC is smart enough to realize that double wide output single precision multipliers can be used. For
2476
example, the instruction ``MUL'' on the x86 processor can multiply two 32-bit values and produce a 64-bit result.
2466
First we determine (line @30,if@) if the Comba method can be used first since it's faster. The conditions for
2467
sing the Comba routine are that min$(a.used, b.used) < \delta$ and the number of digits of output is less than
2468
\textbf{MP\_WARRAY}. This new constant is used to control the stack usage in the Comba routines. By default it is
2469
set to $\delta$ but can be reduced when memory is at a premium.
2471
If we cannot use the Comba method we proceed to setup the baseline routine. We allocate the the destination mp\_int
2472
$t$ (line @36,init@) to the exact size of the output to avoid further re--allocations. At this point we now
2473
begin the $O(n^2)$ loop.
2475
This implementation of multiplication has the caveat that it can be trimmed to only produce a variable number of
2476
digits as output. In each iteration of the outer loop the $pb$ variable is set (line @48,MIN@) to the maximum
2477
number of inner loop iterations.
2479
Inside the inner loop we calculate $\hat r$ as the mp\_word product of the two mp\_digits and the addition of the
2480
carry from the previous iteration. A particularly important observation is that most modern optimizing
2481
C compilers (GCC for instance) can recognize that a $N \times N \rightarrow 2N$ multiplication is all that
2482
is required for the product. In x86 terms for example, this means using the MUL instruction.
2484
Each digit of the product is stored in turn (line @68,tmpt@) and the carry propagated (line @71,>>@) to the
2478
2487
\subsection{Faster Multiplication by the ``Comba'' Method}
2481
One of the huge drawbacks of the ``baseline'' algorithms is that at the $O(n^2)$ level the carry must be computed and propagated upwards. This
2482
makes the nested loop very sequential and hard to unroll and implement in parallel. The ``Comba'' \cite{COMBA} method is named after little known
2483
(\textit{in cryptographic venues}) Paul G. Comba who described a method of implementing fast multipliers that do not require nested
2484
carry fixup operations. As an interesting aside it seems that Paul Barrett describes a similar technique in
2485
his 1986 paper \cite{BARRETT} written five years before.
2487
At the heart of the Comba technique is once again the long-hand algorithm. Except in this case a slight twist is placed on how
2488
the columns of the result are produced. In the standard long-hand algorithm rows of products are produced then added together to form the
2489
final result. In the baseline algorithm the columns are added together after each iteration to get the result instantaneously.
2491
In the Comba algorithm the columns of the result are produced entirely independently of each other. That is at the $O(n^2)$ level a
2492
simple multiplication and addition step is performed. The carries of the columns are propagated after the nested loop to reduce the amount
2493
of work requiored. Succintly the first step of the algorithm is to compute the product vector $\vec x$ as follows.
2490
One of the huge drawbacks of the ``baseline'' algorithms is that at the $O(n^2)$ level the carry must be
2491
computed and propagated upwards. This makes the nested loop very sequential and hard to unroll and implement
2492
in parallel. The ``Comba'' \cite{COMBA} method is named after little known (\textit{in cryptographic venues}) Paul G.
2493
Comba who described a method of implementing fast multipliers that do not require nested carry fixup operations. As an
2494
interesting aside it seems that Paul Barrett describes a similar technique in his 1986 paper \cite{BARRETT} written
2497
At the heart of the Comba technique is once again the long-hand algorithm. Except in this case a slight
2498
twist is placed on how the columns of the result are produced. In the standard long-hand algorithm rows of products
2499
are produced then added together to form the final result. In the baseline algorithm the columns are added together
2500
after each iteration to get the result instantaneously.
2502
In the Comba algorithm the columns of the result are produced entirely independently of each other. That is at
2503
the $O(n^2)$ level a simple multiplication and addition step is performed. The carries of the columns are propagated
2504
after the nested loop to reduce the amount of work requiored. Succintly the first step of the algorithm is to compute
2505
the product vector $\vec x$ as follows.
2495
2507
\begin{equation}
2496
2508
\vec x_n = \sum_{i+j = n} a_ib_j, \forall n \in \lbrace 0, 1, 2, \ldots, i + j \rbrace
2584
2596
\textbf{Input}. mp\_int $a$, mp\_int $b$ and an integer $digs$ \\
2585
2597
\textbf{Output}. $c \leftarrow \vert a \vert \cdot \vert b \vert \mbox{ (mod }\beta^{digs}\mbox{)}$. \\
2587
Place an array of \textbf{MP\_WARRAY} double precision digits named $\hat W$ on the stack. \\
2599
Place an array of \textbf{MP\_WARRAY} single precision digits named $W$ on the stack. \\
2588
2600
1. If $c.alloc < digs$ then grow $c$ to $digs$ digits. (\textit{mp\_grow}) \\
2589
2601
2. If step 1 failed return(\textit{MP\_MEM}).\\
2591
Zero the temporary array $\hat W$. \\
2592
3. for $n$ from $0$ to $digs - 1$ do \\
2593
\hspace{3mm}3.1 $\hat W_n \leftarrow 0$ \\
2595
Compute the columns. \\
2596
4. for $ix$ from $0$ to $a.used - 1$ do \\
2597
\hspace{3mm}4.1 $pb \leftarrow \mbox{min}(b.used, digs - ix)$ \\
2598
\hspace{3mm}4.2 If $pb < 1$ then goto step 5. \\
2599
\hspace{3mm}4.3 for $iy$ from $0$ to $pb - 1$ do \\
2600
\hspace{6mm}4.3.1 $\hat W_{ix+iy} \leftarrow \hat W_{ix+iy} + a_{ix}b_{iy}$ \\
2602
Propagate the carries upwards. \\
2603
5. $oldused \leftarrow c.used$ \\
2604
6. $c.used \leftarrow digs$ \\
2605
7. If $digs > 1$ then do \\
2606
\hspace{3mm}7.1. for $ix$ from $1$ to $digs - 1$ do \\
2607
\hspace{6mm}7.1.1 $\hat W_{ix} \leftarrow \hat W_{ix} + \lfloor \hat W_{ix-1} / \beta \rfloor$ \\
2608
\hspace{6mm}7.1.2 $c_{ix - 1} \leftarrow \hat W_{ix - 1} \mbox{ (mod }\beta\mbox{)}$ \\
2610
\hspace{3mm}8.1 $ix \leftarrow 0$ \\
2611
9. $c_{ix} \leftarrow \hat W_{ix} \mbox{ (mod }\beta\mbox{)}$ \\
2613
Zero excess digits. \\
2614
10. If $digs < oldused$ then do \\
2615
\hspace{3mm}10.1 for $n$ from $digs$ to $oldused - 1$ do \\
2616
\hspace{6mm}10.1.1 $c_n \leftarrow 0$ \\
2617
11. Clamp excessive digits of $c$. (\textit{mp\_clamp}) \\
2618
12. Return(\textit{MP\_OKAY}). \\
2603
3. $pa \leftarrow \mbox{MIN}(digs, a.used + b.used)$ \\
2605
4. $\_ \hat W \leftarrow 0$ \\
2606
5. for $ix$ from 0 to $pa - 1$ do \\
2607
\hspace{3mm}5.1 $ty \leftarrow \mbox{MIN}(b.used - 1, ix)$ \\
2608
\hspace{3mm}5.2 $tx \leftarrow ix - ty$ \\
2609
\hspace{3mm}5.3 $iy \leftarrow \mbox{MIN}(a.used - tx, ty + 1)$ \\
2610
\hspace{3mm}5.4 for $iz$ from 0 to $iy - 1$ do \\
2611
\hspace{6mm}5.4.1 $\_ \hat W \leftarrow \_ \hat W + a_{tx+iy}b_{ty-iy}$ \\
2612
\hspace{3mm}5.5 $W_{ix} \leftarrow \_ \hat W (\mbox{mod }\beta)$\\
2613
\hspace{3mm}5.6 $\_ \hat W \leftarrow \lfloor \_ \hat W / \beta \rfloor$ \\
2614
6. $W_{pa} \leftarrow \_ \hat W (\mbox{mod }\beta)$ \\
2616
7. $oldused \leftarrow c.used$ \\
2617
8. $c.used \leftarrow digs$ \\
2618
9. for $ix$ from $0$ to $pa$ do \\
2619
\hspace{3mm}9.1 $c_{ix} \leftarrow W_{ix}$ \\
2620
10. for $ix$ from $pa + 1$ to $oldused - 1$ do \\
2621
\hspace{3mm}10.1 $c_{ix} \leftarrow 0$ \\
2624
12. Return MP\_OKAY. \\
2627
2633
\textbf{Algorithm fast\_s\_mp\_mul\_digs.}
2628
This algorithm performs the unsigned multiplication of $a$ and $b$ using the Comba method limited to $digs$ digits of precision. The algorithm
2629
essentially peforms the same calculation as algorithm s\_mp\_mul\_digs, just much faster.
2631
The array $\hat W$ is meant to be on the stack when the algorithm is used. The size of the array does not change which is ideal. Note also that
2632
unlike algorithm s\_mp\_mul\_digs no temporary mp\_int is required since the result is calculated directly in $\hat W$.
2634
The $O(n^2)$ loop on step four is where the Comba method's advantages begin to show through in comparison to the baseline algorithm. The lack of
2635
a carry variable or propagation in this loop allows the loop to be performed with only single precision multiplication and additions. Now that each
2636
iteration of the inner loop can be performed independent of the others the inner loop can be performed with a high level of parallelism.
2634
This algorithm performs the unsigned multiplication of $a$ and $b$ using the Comba method limited to $digs$ digits of precision.
2636
The outer loop of this algorithm is more complicated than that of the baseline multiplier. This is because on the inside of the
2637
loop we want to produce one column per pass. This allows the accumulator $\_ \hat W$ to be placed in CPU registers and
2638
reduce the memory bandwidth to two \textbf{mp\_digit} reads per iteration.
2640
The $ty$ variable is set to the minimum count of $ix$ or the number of digits in $b$. That way if $a$ has more digits than
2641
$b$ this will be limited to $b.used - 1$. The $tx$ variable is set to the to the distance past $b.used$ the variable
2642
$ix$ is. This is used for the immediately subsequent statement where we find $iy$.
2644
The variable $iy$ is the minimum digits we can read from either $a$ or $b$ before running out. Computing one column at a time
2645
means we have to scan one integer upwards and the other downwards. $a$ starts at $tx$ and $b$ starts at $ty$. In each
2646
pass we are producing the $ix$'th output column and we note that $tx + ty = ix$. As we move $tx$ upwards we have to
2647
move $ty$ downards so the equality remains valid. The $iy$ variable is the number of iterations until
2648
$tx \ge a.used$ or $ty < 0$ occurs.
2650
After every inner pass we store the lower half of the accumulator into $W_{ix}$ and then propagate the carry of the accumulator
2651
into the next round by dividing $\_ \hat W$ by $\beta$.
2638
2653
To measure the benefits of the Comba method over the baseline method consider the number of operations that are required. If the
2639
2654
cost in terms of time of a multiply and addition is $p$ and the cost of a carry propagation is $q$ then a baseline multiplication would require
2644
2659
EXAM,bn_fast_s_mp_mul_digs.c
2646
The memset on line @47,memset@ clears the initial $\hat W$ array to zero in a single step. Like the slower baseline multiplication
2647
implementation a series of aliases (\textit{lines @67, tmpx@, @70, tmpy@ and @75,_W@}) are used to simplify the inner $O(n^2)$ loop.
2648
In this case a new alias $\_\hat W$ has been added which refers to the double precision columns offset by $ix$ in each pass.
2650
The inner loop on lines @83,for@, @84,mp_word@ and @85,}@ is where the algorithm will spend the majority of the time, which is why it has been
2651
stripped to the bones of any extra baggage\footnote{Hence the pointer aliases.}. On x86 processors the multiplication and additions amount to at the
2652
very least five instructions (\textit{two loads, two additions, one multiply}) while on the ARMv4 processors they amount to only three
2653
(\textit{one load, one store, one multiply-add}). For both of the x86 and ARMv4 processors the GCC compiler performs a good job at unrolling the loop
2654
and scheduling the instructions so there are very few dependency stalls.
2656
In theory the difference between the baseline and comba algorithms is a mere $O(qn)$ time difference. However, in the $O(n^2)$ nested loop of the
2657
baseline method there are dependency stalls as the algorithm must wait for the multiplier to finish before propagating the carry to the next
2658
digit. As a result fewer of the often multiple execution units\footnote{The AMD Athlon has three execution units and the Intel P4 has four.} can
2659
be simultaneously used.
2661
As per the pseudo--code we first calculate $pa$ (line @47,MIN@) as the number of digits to output. Next we begin the outer loop
2662
to produce the individual columns of the product. We use the two aliases $tmpx$ and $tmpy$ (lines @61,tmpx@, @62,tmpy@) to point
2663
inside the two multiplicands quickly.
2665
The inner loop (lines @70,for@ to @72,}@) of this implementation is where the tradeoff come into play. Originally this comba
2666
implementation was ``row--major'' which means it adds to each of the columns in each pass. After the outer loop it would then fix
2667
the carries. This was very fast except it had an annoying drawback. You had to read a mp\_word and two mp\_digits and write
2668
one mp\_word per iteration. On processors such as the Athlon XP and P4 this did not matter much since the cache bandwidth
2669
is very high and it can keep the ALU fed with data. It did, however, matter on older and embedded cpus where cache is often
2670
slower and also often doesn't exist. This new algorithm only performs two reads per iteration under the assumption that the
2671
compiler has aliased $\_ \hat W$ to a CPU register.
2673
After the inner loop we store the current accumulator in $W$ and shift $\_ \hat W$ (lines @75,W[ix]@, @78,>>@) to forward it as
2674
a carry for the next pass. After the outer loop we use the final carry (line @82,W[ix]@) as the last digit of the product.
2661
2676
\subsection{Polynomial Basis Multiplication}
2662
2677
To break the $O(n^2)$ barrier in multiplication requires a completely different look at integer multiplication. In the following algorithms
2977
2992
EXAM,bn_mp_toom_mul.c
2994
The first obvious thing to note is that this algorithm is complicated. The complexity is worth it if you are multiplying very
2995
large numbers. For example, a 10,000 digit multiplication takes approximaly 99,282,205 fewer single precision multiplications with
2996
Toom--Cook than a Comba or baseline approach (this is a savings of more than 99$\%$). For most ``crypto'' sized numbers this
2997
algorithm is not practical as Karatsuba has a much lower cutoff point.
2999
First we split $a$ and $b$ into three roughly equal portions. This has been accomplished (lines @40,mod@ to @69,rshd@) with
3000
combinations of mp\_rshd() and mp\_mod\_2d() function calls. At this point $a = a2 \cdot \beta^2 + a1 \cdot \beta + a0$ and similiarly
3003
Next we compute the five points $w0, w1, w2, w3$ and $w4$. Recall that $w0$ and $w4$ can be computed directly from the portions so
3004
we get those out of the way first (lines @72,mul@ and @77,mul@). Next we compute $w1, w2$ and $w3$ using Horners method.
3006
After this point we solve for the actual values of $w1, w2$ and $w3$ by reducing the $5 \times 5$ system which is relatively
2980
3009
\subsection{Signed Multiplication}
2981
3010
Now that algorithms to handle multiplications of every useful dimensions have been developed, a rather simple finishing touch is required. So far all
2982
3011
of the multiplication algorithms have been unsigned multiplications which leaves only a signed multiplication algorithm to be established.
2984
\newpage\begin{figure}[!here]
3013
\begin{figure}[!here]
2987
3016
\begin{tabular}{l}
3122
3150
EXAM,bn_s_mp_sqr.c
3124
Inside the outer loop (\textit{see line @32,for@}) the square term is calculated on line @35,r =@. Line @42,>>@ extracts the carry from the square
3125
term. Aliases for $a_{ix}$ and $t_{ix+iy}$ are initialized on lines @45,tmpx@ and @48,tmpt@ respectively. The doubling is performed using two
3126
additions (\textit{see line @57,r + r@}) since it is usually faster than shifting,if not at least as fast.
3152
Inside the outer loop (line @32,for@) the square term is calculated on line @35,r =@. The carry (line @42,>>@) has been
3153
extracted from the mp\_word accumulator using a right shift. Aliases for $a_{ix}$ and $t_{ix+iy}$ are initialized
3154
(lines @45,tmpx@ and @48,tmpt@) to simplify the inner loop. The doubling is performed using two
3155
additions (line @57,r + r@) since it is usually faster than shifting, if not at least as fast.
3157
The important observation is that the inner loop does not begin at $iy = 0$ like for multiplication. As such the inner loops
3158
get progressively shorter as the algorithm proceeds. This is what leads to the savings compared to using a multiplication to
3128
3161
\subsection{Faster Squaring by the ``Comba'' Method}
3129
3162
A major drawback to the baseline method is the requirement for single precision shifting inside the $O(n^2)$ nested loop. Squaring has an additional
3147
3180
\textbf{Input}. mp\_int $a$ \\
3148
3181
\textbf{Output}. $b \leftarrow a^2$ \\
3150
Place two arrays of \textbf{MP\_WARRAY} mp\_words named $\hat W$ and $\hat {X}$ on the stack. \\
3183
Place an array of \textbf{MP\_WARRAY} mp\_digits named $W$ on the stack. \\
3151
3184
1. If $b.alloc < 2a.used + 1$ then grow $b$ to $2a.used + 1$ digits. (\textit{mp\_grow}). \\
3152
3185
2. If step 1 failed return(\textit{MP\_MEM}). \\
3153
3. for $ix$ from $0$ to $2a.used + 1$ do \\
3154
\hspace{3mm}3.1 $\hat W_{ix} \leftarrow 0$ \\
3155
\hspace{3mm}3.2 $\hat {X}_{ix} \leftarrow 0$ \\
3156
4. for $ix$ from $0$ to $a.used - 1$ do \\
3157
\hspace{3mm}Compute the square.\\
3158
\hspace{3mm}4.1 $\hat {X}_{ix+ix} \leftarrow \left ( a_{ix} \right )^2$ \\
3160
\hspace{3mm}Compute the double products.\\
3161
\hspace{3mm}4.2 for $iy$ from $ix + 1$ to $a.used - 1$ do \\
3162
\hspace{6mm}4.2.1 $\hat W_{ix+iy} \leftarrow \hat W_{ix+iy} + a_{ix}a_{iy}$ \\
3163
5. $oldused \leftarrow b.used$ \\
3164
6. $b.used \leftarrow 2a.used + 1$ \\
3166
Double the products and propagate the carries simultaneously. \\
3167
7. $\hat W_0 \leftarrow 2 \hat W_0 + \hat {X}_0$ \\
3168
8. for $ix$ from $1$ to $2a.used$ do \\
3169
\hspace{3mm}8.1 $\hat W_{ix} \leftarrow 2 \hat W_{ix} + \hat {X}_{ix}$ \\
3170
\hspace{3mm}8.2 $\hat W_{ix} \leftarrow \hat W_{ix} + \lfloor \hat W_{ix - 1} / \beta \rfloor$ \\
3171
\hspace{3mm}8.3 $b_{ix-1} \leftarrow W_{ix-1} \mbox{ (mod }\beta\mbox{)}$ \\
3172
9. $b_{2a.used} \leftarrow \hat W_{2a.used} \mbox{ (mod }\beta\mbox{)}$ \\
3173
10. if $2a.used + 1 < oldused$ then do \\
3174
\hspace{3mm}10.1 for $ix$ from $2a.used + 1$ to $oldused$ do \\
3175
\hspace{6mm}10.1.1 $b_{ix} \leftarrow 0$ \\
3176
11. Clamp excess digits from $b$. (\textit{mp\_clamp}) \\
3177
12. Return(\textit{MP\_OKAY}). \\
3187
3. $pa \leftarrow 2 \cdot a.used$ \\
3188
4. $\hat W1 \leftarrow 0$ \\
3189
5. for $ix$ from $0$ to $pa - 1$ do \\
3190
\hspace{3mm}5.1 $\_ \hat W \leftarrow 0$ \\
3191
\hspace{3mm}5.2 $ty \leftarrow \mbox{MIN}(a.used - 1, ix)$ \\
3192
\hspace{3mm}5.3 $tx \leftarrow ix - ty$ \\
3193
\hspace{3mm}5.4 $iy \leftarrow \mbox{MIN}(a.used - tx, ty + 1)$ \\
3194
\hspace{3mm}5.5 $iy \leftarrow \mbox{MIN}(iy, \lfloor \left (ty - tx + 1 \right )/2 \rfloor)$ \\
3195
\hspace{3mm}5.6 for $iz$ from $0$ to $iz - 1$ do \\
3196
\hspace{6mm}5.6.1 $\_ \hat W \leftarrow \_ \hat W + a_{tx + iz}a_{ty - iz}$ \\
3197
\hspace{3mm}5.7 $\_ \hat W \leftarrow 2 \cdot \_ \hat W + \hat W1$ \\
3198
\hspace{3mm}5.8 if $ix$ is even then \\
3199
\hspace{6mm}5.8.1 $\_ \hat W \leftarrow \_ \hat W + \left ( a_{\lfloor ix/2 \rfloor}\right )^2$ \\
3200
\hspace{3mm}5.9 $W_{ix} \leftarrow \_ \hat W (\mbox{mod }\beta)$ \\
3201
\hspace{3mm}5.10 $\hat W1 \leftarrow \lfloor \_ \hat W / \beta \rfloor$ \\
3203
6. $oldused \leftarrow b.used$ \\
3204
7. $b.used \leftarrow 2 \cdot a.used$ \\
3205
8. for $ix$ from $0$ to $pa - 1$ do \\
3206
\hspace{3mm}8.1 $b_{ix} \leftarrow W_{ix}$ \\
3207
9. for $ix$ from $pa$ to $oldused - 1$ do \\
3208
\hspace{3mm}9.1 $b_{ix} \leftarrow 0$ \\
3209
10. Clamp excess digits from $b$. (\textit{mp\_clamp}) \\
3210
11. Return(\textit{MP\_OKAY}). \\
3185
3218
\textbf{Algorithm fast\_s\_mp\_sqr.}
3186
This algorithm computes the square of an input using the Comba technique. It is designed to be a replacement for algorithm s\_mp\_sqr when
3187
the number of input digits is less than \textbf{MP\_WARRAY} and less than $\delta \over 2$.
3189
This routine requires two arrays of mp\_words to be placed on the stack. The first array $\hat W$ will hold the double products and the second
3190
array $\hat X$ will hold the squares. Though only at most $MP\_WARRAY \over 2$ words of $\hat X$ are used, it has proven faster on most
3191
processors to simply make it a full size array.
3193
The loop on step 3 will zero the two arrays to prepare them for the squaring step. Step 4.1 computes the squares of the product. Note how
3194
it simply assigns the value into the $\hat X$ array. The nested loop on step 4.2 computes the doubles of the products. This loop
3195
computes the sum of the products for each column. They are not doubled until later.
3197
After the squaring loop, the products stored in $\hat W$ musted be doubled and the carries propagated forwards. It makes sense to do both
3198
operations at the same time. The expression $\hat W_{ix} \leftarrow 2 \hat W_{ix} + \hat {X}_{ix}$ computes the sum of the double product and the
3219
This algorithm computes the square of an input using the Comba technique. It is designed to be a replacement for algorithm
3220
s\_mp\_sqr when the number of input digits is less than \textbf{MP\_WARRAY} and less than $\delta \over 2$.
3221
This algorithm is very similar to the Comba multiplier except with a few key differences we shall make note of.
3223
First, we have an accumulator and carry variables $\_ \hat W$ and $\hat W1$ respectively. This is because the inner loop
3224
products are to be doubled. If we had added the previous carry in we would be doubling too much. Next we perform an
3225
addition MIN condition on $iy$ (step 5.5) to prevent overlapping digits. For example, $a_3 \cdot a_5$ is equal
3226
$a_5 \cdot a_3$. Whereas in the multiplication case we would have $5 < a.used$ and $3 \ge 0$ is maintained since we double the sum
3227
of the products just outside the inner loop we have to avoid doing this. This is also a good thing since we perform
3228
fewer multiplications and the routine ends up being faster.
3230
Finally the last difference is the addition of the ``square'' term outside the inner loop (step 5.8). We add in the square
3231
only to even outputs and it is the square of the term at the $\lfloor ix / 2 \rfloor$ position.
3201
3233
EXAM,bn_fast_s_mp_sqr.c
3235
This implementation is essentially a copy of Comba multiplication with the appropriate changes added to make it faster for
3236
the special case of squaring.
3204
3238
\subsection{Polynomial Basis Squaring}
3205
3239
The same algorithm that performs optimal polynomial basis multiplication can be used to perform polynomial basis squaring. The minor exception
3312
3345
is exactly at the point where Comba squaring can no longer be used (\textit{128 digits}). On slower processors such as the Intel P4
3313
3346
it is actually below the Comba limit (\textit{at 110 digits}).
3315
This routine uses the same error trap coding style as mp\_karatsuba\_sqr. As the temporary variables are initialized errors are redirected to
3316
the error trap higher up. If the algorithm completes without error the error code is set to \textbf{MP\_OKAY} and mp\_clears are executed normally.
3318
\textit{Last paragraph sucks. re-write! -- Tom}
3348
This routine uses the same error trap coding style as mp\_karatsuba\_sqr. As the temporary variables are initialized errors are
3349
redirected to the error trap higher up. If the algorithm completes without error the error code is set to \textbf{MP\_OKAY} and
3350
mp\_clears are executed normally.
3320
3352
\subsection{Toom-Cook Squaring}
3321
3353
The Toom-Cook squaring algorithm mp\_toom\_sqr is heavily based on the algorithm mp\_toom\_mul with the exception that squarings are used
3322
instead of multiplication to find the five relations.. The reader is encouraged to read the description of the latter algorithm and try to
3354
instead of multiplication to find the five relations. The reader is encouraged to read the description of the latter algorithm and try to
3323
3355
derive their own Toom-Cook squaring algorithm.
3325
3357
\subsection{High Level Squaring}
added
removed
libtommath/tommath.src
