Below is the file 'tommath.tex' from this revision. You can also download the file.
\documentclass[b5paper]{book} \usepackage{hyperref} \usepackage{makeidx} \usepackage{amssymb} \usepackage{color} \usepackage{alltt} \usepackage{graphicx} \usepackage{layout} \def\union{\cup} \def\intersect{\cap} \def\getsrandom{\stackrel{\rm R}{\gets}} \def\cross{\times} \def\cat{\hspace{0.5em} \| \hspace{0.5em}} \def\catn{$\|$} \def\divides{\hspace{0.3em} | \hspace{0.3em}} \def\nequiv{\not\equiv} \def\approx{\raisebox{0.2ex}{\mbox{\small $\sim$}}} \def\lcm{{\rm lcm}} \def\gcd{{\rm gcd}} \def\log{{\rm log}} \def\ord{{\rm ord}} \def\abs{{\mathit abs}} \def\rep{{\mathit rep}} \def\mod{{\mathit\ mod\ }} \renewcommand{\pmod}[1]{\ ({\rm mod\ }{#1})} \newcommand{\floor}[1]{\left\lfloor{#1}\right\rfloor} \newcommand{\ceil}[1]{\left\lceil{#1}\right\rceil} \def\Or{{\rm\ or\ }} \def\And{{\rm\ and\ }} \def\iff{\hspace{1em}\Longleftrightarrow\hspace{1em}} \def\implies{\Rightarrow} \def\undefined{{\rm ``undefined"}} \def\Proof{\vspace{1ex}\noindent {\bf Proof:}\hspace{1em}} \let\oldphi\phi \def\phi{\varphi} \def\Pr{{\rm Pr}} \newcommand{\str}[1]{{\mathbf{#1}}} \def\F{{\mathbb F}} \def\N{{\mathbb N}} \def\Z{{\mathbb Z}} \def\R{{\mathbb R}} \def\C{{\mathbb C}} \def\Q{{\mathbb Q}} \definecolor{DGray}{gray}{0.5} \newcommand{\emailaddr}[1]{\mbox{$<${#1}$>$}} \def\twiddle{\raisebox{0.3ex}{\mbox{\tiny $\sim$}}} \def\gap{\vspace{0.5ex}} \makeindex \begin{document} \frontmatter \pagestyle{empty} \title{Implementing Multiple Precision Arithmetic \\ ~ \\ Draft Edition } \author{\mbox{ %\begin{small} \begin{tabular}{c} Tom St Denis \\ Algonquin College \\ \\ Mads Rasmussen \\ Open Communications Security \\ \\ Greg Rose \\ QUALCOMM Australia \\ \end{tabular} %\end{small} } } \maketitle This text has been placed in the public domain. This text corresponds to the v0.30 release of the LibTomMath project. \begin{alltt} Tom St Denis 111 Banning Rd Ottawa, Ontario K2L 1C3 Canada Phone: 1-613-836-3160 Email: tomstdenis@iahu.ca \end{alltt} This text is formatted to the international B5 paper size of 176mm wide by 250mm tall using the \LaTeX{} {\em book} macro package and the Perl {\em booker} package. \tableofcontents \listoffigures \chapter*{Prefaces to the Draft Edition} I started this text in April 2003 to complement my LibTomMath library. That is, explain how to implement the functions contained in LibTomMath. The goal is to have a textbook that any Computer Science student can use when implementing their own multiple precision arithmetic. The plan I wanted to follow was flesh out all the ideas and concepts I had floating around in my head and then work on it afterwards refining a little bit at a time. Chance 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 text. Choosing to not waste any time I dove right into the project even before my spring semester was finished. I wrote a bit off and on at first. The moment my exams were finished I jumped into long 12 to 16 hour days. The result after only a couple of months was a ten chapter, three hundred page draft that I quickly had distributed to anyone who wanted 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 managed to grab a certain level of attention having people from around the world ask me for copies of the text was certain rewarding. 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. 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 finished yet. I haven't given up on the project, only had some setbacks. First O'Reilly declined to publish the text then 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 onto finishing the book not securing a contract. 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. Even the simplest introductory material has to be lined with references and figures. A lot of the text has to be re-written from point form to prose form to ensure an easier read. Why am I doing all this work for free then? Simple. My philosophy is quite simply ``Open Source. Open Academia. Open Minds'' which means that to achieve a goal of open minds, that is, people willing to accept new ideas and explore the unknown you have to make available material they can access freely without hinderance. I've been writing free software since I was about sixteen but only recently have I hit upon software that people have come to depend upon. I started LibTomCrypt in December 2001 and now several major companies use it as integral portions of their software. Several educational institutions use it as a matter of course and many freelance developers use it as part of their projects. To further my contributions I started the LibTomMath project in December 2002 aimed at providing multiple precision arithmetic routines that students could learn from. That is write routines that are not only easy to understand and follow but provide quite impressive performance considering they are all in standard portable ISO C. The second leg of my philosophy is ``Open Academia'' which is where this textbook comes in. In the end, when all is said and done the text will be useable by educational institutions as a reference on multiple precision arithmetic. At this time I feel I should share a little information about myself. The most common question I was asked at Crypto'03, perhaps just out of professional courtesy, was which school I either taught at or attended. The unfortunate truth is that I neither teach at or attend a school of academic reputation. I'm currently at Algonquin College which is what I'd like to call ``somewhat academic but mostly vocational'' college. In otherwords, job training. I'm a 21 year old computer science student mostly self-taught in the areas I am aware of (which includes a half-dozen computer science fields, a few fields of mathematics and some English). I look forward to teaching someday but I am still far off from that goal. Now it would be improper for me to not introduce the rest of the texts co-authors. While they are only contributing corrections and editorial feedback their support has been tremendously helpful in presenting the concepts laid out in the text so far. Greg has always been there for me. He has tracked my LibTom projects since their inception and even sent cheques to help pay tuition from time to time. His background has provided a wonderful source to bounce ideas off of and improve the quality of my writing. Mads is another fellow who has just ``been there''. I don't even recall what his interest in the LibTom projects is but I'm definitely glad he has been around. His ability to catch logical errors in my written English have saved me on several occasions to say the least. 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 been getting the feeling that people are starting to use my text and I owe them some updated material. My current tenative plan is to edit one chapter every two weeks starting January 4th. It seems insane but my lower course load at college 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 people who will take it. \begin{flushright} Tom St Denis \end{flushright} \newpage I found the opportunity to work with Tom appealing for several reasons, not only could I broaden my own horizons, but also contribute to educate others facing the problem of having to handle big number mathematical calculations. This book is Tom's child and he has been caring and fostering the project ever since the beginning with a clear mind of how he wanted the project to turn out. I have helped by proofreading the text and we have had several discussions about the layout and language used. I hold a masters degree in cryptography from the University of Southern Denmark and have always been interested in the practical aspects of cryptography. Having worked in the security consultancy business for several years in S\~{a}o Paulo, Brazil, I have been in touch with a great deal of work in which multiple precision mathematics was needed. Understanding the possibilities for speeding up multiple precision calculations is often very important since we deal with outdated machine architecture where modular reductions, for example, become painfully slow. This text is for people who stop and wonder when first examining algorithms such as RSA for the first time and asks themselves, ``You tell me this is only secure for large numbers, fine; but how do you implement these numbers?'' \begin{flushright} Mads Rasmussen S\~{a}o Paulo - SP Brazil \end{flushright} \newpage It's all because I broke my leg. That just happened to be at about the same time that Tom asked for someone to review the section of the book about Karatsuba multiplication. I was laid up, alone and immobile, and thought ``Why not?'' I vaguely knew what Karatsuba multiplication was, but not really, so I thought I could help, learn, and stop myself from watching daytime cable TV, all at once. At the time of writing this, I've still not met Tom or Mads in meatspace. I've been following Tom's progress since his first splash on the sci.crypt Usenet news group. I watched him go from a clueless newbie, to the cryptographic equivalent of a reformed smoker, to a real contributor to the field, over a period of about two years. I've been impressed with his obvious intelligence, and astounded by his productivity. Of course, he's young enough to be my own child, so he doesn't have my problems with staying awake. When I reviewed that single section of the book, in its very earliest form, I was very pleasantly surprised. So I decided to collaborate more fully, and at least review all of it, and perhaps write some bits too. There's still a long way to go with it, and I have watched a number of close friends go through the mill of publication, so I think that the way to go is longer than Tom thinks it is. Nevertheless, it's a good effort, and I'm pleased to be involved with it. \begin{flushright} Greg Rose, Sydney, Australia, June 2003. \end{flushright} \mainmatter \pagestyle{headings} \chapter{Introduction} \section{Multiple Precision Arithmetic} \subsection{What is Multiple Precision Arithmetic?} When we think of long-hand arithmetic such as addition or multiplication we rarely consider the fact that we instinctively raise or lower the precision of the numbers we are dealing with. For example, in decimal we almost immediate can reason that $7$ times $6$ is $42$. However, $42$ has two digits of precision as opposed to one digit we started with. Further multiplications of say $3$ result in a larger precision result $126$. In these few examples we have multiple precisions for the numbers we are working with. Despite the various levels of precision a single subset\footnote{With the occasional optimization.} of algorithms can be designed to accomodate them. By way of comparison a fixed or single precision operation would lose precision on various operations. For example, in the decimal system with fixed precision $6 \cdot 7 = 2$. Essentially at the heart of computer based multiple precision arithmetic are the same long-hand algorithms taught in schools to manually add, subtract, multiply and divide. \subsection{The Need for Multiple Precision Arithmetic} The most prevalent need for multiple precision arithmetic, often referred to as ``bignum'' math, is within the implementation of public-key cryptography algorithms. Algorithms such as RSA \cite{RSAREF} and Diffie-Hellman \cite{DHREF} require integers of significant magnitude to resist known cryptanalytic attacks. For example, at the time of this writing a typical RSA modulus would be at least greater than $10^{309}$. However, modern programming languages such as ISO C \cite{ISOC} and Java \cite{JAVA} only provide instrinsic support for integers which are relatively small and single precision. \begin{figure}[!here] \begin{center} \begin{tabular}{|r|c|} \hline \textbf{Data Type} & \textbf{Range} \\ \hline char & $-128 \ldots 127$ \\ \hline short & $-32768 \ldots 32767$ \\ \hline long & $-2147483648 \ldots 2147483647$ \\ \hline long long & $-9223372036854775808 \ldots 9223372036854775807$ \\ \hline \end{tabular} \end{center} \caption{Typical Data Types for the C Programming Language} \label{fig:ISOC} \end{figure} The largest data type guaranteed to be provided by the ISO C programming language\footnote{As per the ISO C standard. However, each compiler vendor is allowed to augment the precision as they see fit.} can only represent values up to $10^{19}$ as shown in figure \ref{fig:ISOC}. On its own the C language is insufficient to accomodate the magnitude required for the problem at hand. An RSA modulus of magnitude $10^{19}$ could be trivially factored\footnote{A Pollard-Rho factoring would take only $2^{16}$ time.} on the average desktop computer, rendering any protocol based on the algorithm insecure. Multiple precision algorithms solve this very problem by extending the range of representable integers while using single precision data types. Most advancements in fast multiple precision arithmetic stem from the need for faster and more efficient cryptographic primitives. Faster modular reduction and exponentiation algorithms such as Barrett's algorithm, which have appeared in various cryptographic journals, can render algorithms such as RSA and Diffie-Hellman more efficient. In fact, several major companies such as RSA Security, Certicom and Entrust have built entire product lines on the implementation and deployment of efficient algorithms. However, cryptography is not the only field of study that can benefit from fast multiple precision integer routines. Another auxiliary use of multiple precision integers is high precision floating point data types. The basic IEEE \cite{IEEE} standard floating point type is made up of an integer mantissa $q$, an exponent $e$ and a sign bit $s$. Numbers are given in the form $n = q \cdot b^e \cdot -1^s$ where $b = 2$ is the most common base for IEEE. Since IEEE floating point is meant to be implemented in hardware the precision of the mantissa is often fairly small (\textit{23, 48 and 64 bits}). The mantissa is merely an integer and a multiple precision integer could be used to create a mantissa of much larger precision than hardware alone can efficiently support. This approach could be useful where scientific applications must minimize the total output error over long calculations. Another use for large integers is within arithmetic on polynomials of large characteristic (i.e. $GF(p)[x]$ for large $p$). In fact the library discussed within this text has already been used to form a polynomial basis library\footnote{See \url{http://poly.libtomcrypt.org} for more details.}. \subsection{Benefits of Multiple Precision Arithmetic} \index{precision} The benefit of multiple precision representations over single or fixed precision representations is that no precision is lost while representing the result of an operation which requires excess precision. For example, the product of two $n$-bit integers requires at least $2n$ bits of precision to be represented faithfully. A multiple precision algorithm would augment the precision of the destination to accomodate the result while a single precision system would truncate excess bits to maintain a fixed level of precision. It is possible to implement algorithms which require large integers with fixed precision algorithms. For example, elliptic curve cryptography (\textit{ECC}) is often implemented on smartcards by fixing the precision of the integers to the maximum size the system will ever need. Such an approach can lead to vastly simpler algorithms which can accomodate the integers required even if the host platform cannot natively accomodate them\footnote{For example, the average smartcard processor has an 8 bit accumulator.}. However, as efficient as such an approach may be, the resulting source code is not normally very flexible. It cannot, at runtime, accomodate inputs of higher magnitude than the designer anticipated. Multiple precision algorithms have the most overhead of any style of arithmetic. For the the most part the overhead can be kept to a minimum with careful planning, but overall, it is not well suited for most memory starved platforms. However, multiple precision algorithms do offer the most flexibility in terms of the magnitude of the inputs. That is, the same algorithms based on multiple precision integers can accomodate any reasonable size input without the designer's explicit forethought. This leads to lower cost of ownership for the code as it only has to be written and tested once. \section{Purpose of This Text} The purpose of this text is to instruct the reader regarding how to implement efficient multiple precision algorithms. That is to not only explain a limited subset of the core theory behind the algorithms but also the various ``house keeping'' elements that are neglected by authors of other texts on the subject. Several well reknowned texts \cite{TAOCPV2,HAC} give considerably detailed explanations of the theoretical aspects of algorithms and often very little information regarding the practical implementation aspects. In most cases how an algorithm is explained and how it is actually implemented are two very different concepts. For example, the Handbook of Applied Cryptography (\textit{HAC}), algorithm 14.7 on page 594, gives a relatively simple algorithm for performing multiple precision integer addition. However, the description lacks any discussion concerning the fact that the two integer inputs may be of differing magnitudes. As a result the implementation is not as simple as the text would lead people to believe. Similarly the division routine (\textit{algorithm 14.20, pp. 598}) does not discuss how to handle sign or handle the dividend's decreasing magnitude in the main loop (\textit{step \#3}). Both texts also do not discuss several key optimal algorithms required such as ``Comba'' and Karatsuba multipliers and fast modular inversion, which we consider practical oversights. These optimal algorithms are vital to achieve any form of useful performance in non-trivial applications. To solve this problem the focus of this text is on the practical aspects of implementing a multiple precision integer package. As a case study the ``LibTomMath''\footnote{Available at \url{http://math.libtomcrypt.org}} package is used to demonstrate algorithms with real implementations\footnote{In the ISO C programming language.} that have been field tested and work very well. The LibTomMath library is freely available on the Internet for all uses and this text discusses a very large portion of the inner workings of the library. The algorithms that are presented will always include at least one ``pseudo-code'' description followed by the actual C source code that implements the algorithm. The pseudo-code can be used to implement the same algorithm in other programming languages as the reader sees fit. This text shall also serve as a walkthrough of the creation of multiple precision algorithms from scratch. Showing the reader how the algorithms fit together as well as where to start on various taskings. \section{Discussion and Notation} \subsection{Notation} A multiple precision integer of $n$-digits shall be denoted as $x = (x_{n-1} ... x_1 x_0)_{ \beta }$ and represent the integer $x \equiv \sum_{i=0}^{n-1} x_i\beta^i$. The elements of the array $x$ are said to be the radix $\beta$ digits of the integer. For example, $x = (1,2,3)_{10}$ would represent the integer $1\cdot 10^2 + 2\cdot10^1 + 3\cdot10^0 = 123$. \index{mp\_int} The term ``mp\_int'' shall refer to a composite structure which contains the digits of the integer it represents, as well as auxilary data required to manipulate the data. These additional members are discussed further in section \ref{sec:MPINT}. For the purposes of this text a ``multiple precision integer'' and an ``mp\_int'' are assumed to be synonymous. When an algorithm is specified to accept an mp\_int variable it is assumed the various auxliary data members are present as well. An expression of the type \textit{variablename.item} implies that it should evaluate to the member named ``item'' of the variable. For example, a string of characters may have a member ``length'' which would evaluate to the number of characters in the string. If the string $a$ equals ``hello'' then it follows that $a.length = 5$. For certain discussions more generic algorithms are presented to help the reader understand the final algorithm used to solve a given problem. When an algorithm is described as accepting an integer input it is assumed the input is a plain integer with no additional multiple-precision members. That is, algorithms that use integers as opposed to mp\_ints as inputs do not concern themselves with the housekeeping operations required such as memory management. These algorithms will be used to establish the relevant theory which will subsequently be used to describe a multiple precision algorithm to solve the same problem. \subsection{Precision Notation} For the purposes of this text a single precision variable must be able to represent integers in the range $0 \le x < q \beta$ while a double precision variable must be able to represent integers in the range $0 \le x < q \beta^2$. The variable $\beta$ represents the radix of a single digit of a multiple precision integer and must be of the form $q^p$ for $q, p \in \Z^+$. The extra radix-$q$ factor allows additions and subtractions to proceed without truncation of the carry. Since all modern computers are binary, it is assumed that $q$ is two, for all intents and purposes. \index{mp\_digit} \index{mp\_word} Within the source code that will be presented for each algorithm, the data type \textbf{mp\_digit} will represent a single precision integer type, while, the data type \textbf{mp\_word} will represent a double precision integer type. In several algorithms (notably the Comba routines) temporary results will be stored in arrays of double precision mp\_words. For the purposes of this text $x_j$ will refer to the $j$'th digit of a single precision array and $\hat x_j$ will refer to the $j$'th digit of a double precision array. Whenever an expression is to be assigned to a double precision variable it is assumed that all single precision variables are promoted to double precision during the evaluation. Expressions that are assigned to a single precision variable are truncated to fit within the precision of a single precision data type. For example, if $\beta = 10^2$ a single precision data type may represent a value in the range $0 \le x < 10^3$, while a double precision data type may represent a value in the range $0 \le x < 10^5$. Let $a = 23$ and $b = 49$ represent two single precision variables. The single precision product shall be written as $c \leftarrow a \cdot b$ while the double precision product shall be written as $\hat c \leftarrow a \cdot b$. In this particular case, $\hat c = 1127$ and $c = 127$. The most significant digit of the product would not fit in a single precision data type and as a result $c \ne \hat c$. \subsection{Algorithm Inputs and Outputs} Within the algorithm descriptions all variables are assumed to be scalars of either single or double precision as indicated. The only exception to this rule is when variables have been indicated to be of type mp\_int. This distinction is important as scalars are often used as array indicies and various other counters. \subsection{Mathematical Expressions} The $\lfloor \mbox{ } \rfloor$ brackets imply an expression truncated to an integer not greater than the expression itself. For example, $\lfloor 5.7 \rfloor = 5$. Similarly the $\lceil \mbox{ } \rceil$ brackets imply an expression rounded to an integer not less than the expression itself. For example, $\lceil 5.1 \rceil = 6$. Typically when the $/$ division symbol is used the intention is to perform an integer division with truncation. For example, $5/2 = 2$ which will often be written as $\lfloor 5/2 \rfloor = 2$ for clarity. When an expression is written as a fraction a real value division is implied, for example ${5 \over 2} = 2.5$. The norm of a multiple precision integer, for example, $\vert \vert x \vert \vert$ will be used to represent the number of digits in the representation of the integer. For example, $\vert \vert 123 \vert \vert = 3$ and $\vert \vert 79452 \vert \vert = 5$. \subsection{Work Effort} \index{big-Oh} To measure the efficiency of the specified algorithms, a modified big-Oh notation is used. In this system all single precision operations are considered to have the same cost\footnote{Except where explicitly noted.}. That is a single precision addition, multiplication and division are assumed to take the same time to complete. While this is generally not true in practice, it will simplify the discussions considerably. Some algorithms have slight advantages over others which is why some constants will not be removed in the notation. For example, a normal baseline multiplication (section \ref{sec:basemult}) requires $O(n^2)$ work while a baseline squaring (section \ref{sec:basesquare}) requires $O({{n^2 + n}\over 2})$ work. In standard big-Oh notation these would both be said to be equivalent to $O(n^2)$. However, in the context of the this text this is not the case as the magnitude of the inputs will typically be rather small. As a result small constant factors in the work effort will make an observable difference in algorithm efficiency. All of the algorithms presented in this text have a polynomial time work level. That is, of the form $O(n^k)$ for $n, k \in \Z^{+}$. This will help make useful comparisons in terms of the speed of the algorithms and how various optimizations will help pay off in the long run. \section{Exercises} Within the more advanced chapters a section will be set aside to give the reader some challenging exercises related to the discussion at hand. These exercises are not designed to be prize winning problems, but instead to be thought provoking. Wherever possible the problems are forward minded, stating problems that will be answered in subsequent chapters. The reader is encouraged to finish the exercises as they appear to get a better understanding of the subject material. That being said, the problems are designed to affirm knowledge of a particular subject matter. Students in particular are encouraged to verify they can answer the problems correctly before moving on. Similar to the exercises of \cite[pp. ix]{TAOCPV2} these exercises are given a scoring system based on the difficulty of the problem. However, unlike \cite{TAOCPV2} the problems do not get nearly as hard. The scoring of these exercises ranges from one (the easiest) to five (the hardest). The following table sumarizes the scoring system used. \begin{figure}[here] \begin{center} \begin{small} \begin{tabular}{|c|l|} \hline $\left [ 1 \right ]$ & An easy problem that should only take the reader a manner of \\ & minutes to solve. Usually does not involve much computer time \\ & to solve. \\ \hline $\left [ 2 \right ]$ & An easy problem that involves a marginal amount of computer \\ & time usage. Usually requires a program to be written to \\ & solve the problem. \\ \hline $\left [ 3 \right ]$ & A moderately hard problem that requires a non-trivial amount \\ & of work. Usually involves trivial research and development of \\ & new theory from the perspective of a student. \\ \hline $\left [ 4 \right ]$ & A moderately hard problem that involves a non-trivial amount \\ & of work and research, the solution to which will demonstrate \\ & a higher mastery of the subject matter. \\ \hline $\left [ 5 \right ]$ & A hard problem that involves concepts that are difficult for a \\ & novice to solve. Solutions to these problems will demonstrate a \\ & complete mastery of the given subject. \\ \hline \end{tabular} \end{small} \end{center} \caption{Exercise Scoring System} \end{figure} Problems at the first level are meant to be simple questions that the reader can answer quickly without programming a solution or devising new theory. These problems are quick tests to see if the material is understood. Problems at the second level are also designed to be easy but will require a program or algorithm to be implemented to arrive at the answer. These two levels are essentially entry level questions. Problems at the third level are meant to be a bit more difficult than the first two levels. The answer is often fairly obvious but arriving at an exacting solution requires some thought and skill. These problems will almost always involve devising a new algorithm or implementing a variation of another algorithm previously presented. Readers who can answer these questions will feel comfortable with the concepts behind the topic at hand. Problems at the fourth level are meant to be similar to those of the level three questions except they will require additional research to be completed. The reader will most likely not know the answer right away, nor will the text provide the exact details of the answer until a subsequent chapter. Problems at the fifth level are meant to be the hardest problems relative to all the other problems in the chapter. People who can correctly answer fifth level problems have a mastery of the subject matter at hand. Often problems will be tied together. The purpose of this is to start a chain of thought that will be discussed in future chapters. The reader is encouraged to answer the follow-up problems and try to draw the relevance of problems. \section{Introduction to LibTomMath} \subsection{What is LibTomMath?} LibTomMath is a free and open source multiple precision integer library written entirely in portable ISO C. By portable it is meant that the library does not contain any code that is computer platform dependent or otherwise problematic to use on any given platform. The library has been successfully tested under numerous operating systems including Unix\footnote{All of these trademarks belong to their respective rightful owners.}, MacOS, Windows, Linux, PalmOS and on standalone hardware such as the Gameboy Advance. The library is designed to contain enough functionality to be able to develop applications such as public key cryptosystems and still maintain a relatively small footprint. \subsection{Goals of LibTomMath} Libraries which obtain the most efficiency are rarely written in a high level programming language such as C. However, even though this library is written entirely in ISO C, considerable care has been taken to optimize the algorithm implementations within the library. Specifically the code has been written to work well with the GNU C Compiler (\textit{GCC}) on both x86 and ARM processors. Wherever possible, highly efficient algorithms, such as Karatsuba multiplication, sliding window exponentiation and Montgomery reduction have been provided to make the library more efficient. Even with the nearly optimal and specialized algorithms that have been included the Application Programing Interface (\textit{API}) has been kept as simple as possible. Often generic place holder routines will make use of specialized algorithms automatically without the developer's specific attention. One such example is the generic multiplication algorithm \textbf{mp\_mul()} which will automatically use Toom--Cook, Karatsuba, Comba or baseline multiplication based on the magnitude of the inputs and the configuration of the library. Making LibTomMath as efficient as possible is not the only goal of the LibTomMath project. Ideally the library should be source compatible with another popular library which makes it more attractive for developers to use. In this case the MPI library was used as a API template for all the basic functions. MPI was chosen because it is another library that fits in the same niche as LibTomMath. Even though LibTomMath uses MPI as the template for the function names and argument passing conventions, it has been written from scratch by Tom St Denis. The project is also meant to act as a learning tool for students, the logic being that no easy-to-follow ``bignum'' library exists which can be used to teach computer science students how to perform fast and reliable multiple precision integer arithmetic. To this end the source code has been given quite a few comments and algorithm discussion points. \section{Choice of LibTomMath} LibTomMath was chosen as the case study of this text not only because the author of both projects is one and the same but for more worthy reasons. Other libraries such as GMP \cite{GMP}, MPI \cite{MPI}, LIP \cite{LIP} and OpenSSL \cite{OPENSSL} have multiple precision integer arithmetic routines but would not be ideal for this text for reasons that will be explained in the following sub-sections. \subsection{Code Base} The LibTomMath code base is all portable ISO C source code. This means that there are no platform dependent conditional segments of code littered throughout the source. This clean and uncluttered approach to the library means that a developer can more readily discern the true intent of a given section of source code without trying to keep track of what conditional code will be used. The code base of LibTomMath is well organized. Each function is in its own separate source code file which allows the reader to find a given function very quickly. On average there are $76$ lines of code per source file which makes the source very easily to follow. By comparison MPI and LIP are single file projects making code tracing very hard. GMP has many conditional code segments which also hinder tracing. When compiled with GCC for the x86 processor and optimized for speed the entire library is approximately $100$KiB\footnote{The notation ``KiB'' means $2^{10}$ octets, similarly ``MiB'' means $2^{20}$ octets.} which is fairly small compared to GMP (over $250$KiB). LibTomMath is slightly larger than MPI (which compiles to about $50$KiB) but LibTomMath is also much faster and more complete than MPI. \subsection{API Simplicity} LibTomMath is designed after the MPI library and shares the API design. Quite often programs that use MPI will build with LibTomMath without change. The function names correlate directly to the action they perform. Almost all of the functions share the same parameter passing convention. The learning curve is fairly shallow with the API provided which is an extremely valuable benefit for the student and developer alike. The LIP library is an example of a library with an API that is awkward to work with. LIP uses function names that are often ``compressed'' to illegible short hand. LibTomMath does not share this characteristic. The GMP library also does not return error codes. Instead it uses a POSIX.1 \cite{POSIX1} signal system where errors are signaled to the host application. This happens to be the fastest approach but definitely not the most versatile. In effect a math error (i.e. invalid input, heap error, etc) can cause a program to stop functioning which is definitely undersireable in many situations. \subsection{Optimizations} While LibTomMath is certainly not the fastest library (GMP often beats LibTomMath by a factor of two) it does feature a set of optimal algorithms for tasks such as modular reduction, exponentiation, multiplication and squaring. GMP and LIP also feature such optimizations while MPI only uses baseline algorithms with no optimizations. GMP lacks a few of the additional modular reduction optimizations that LibTomMath features\footnote{At the time of this writing GMP only had Barrett and Montgomery modular reduction algorithms.}. LibTomMath is almost always an order of magnitude faster than the MPI library at computationally expensive tasks such as modular exponentiation. In the grand scheme of ``bignum'' libraries LibTomMath is faster than the average library and usually slower than the best libraries such as GMP and OpenSSL by only a small factor. \subsection{Portability and Stability} LibTomMath will build ``out of the box'' on any platform equipped with a modern version of the GNU C Compiler (\textit{GCC}). This means that without changes the library will build without configuration or setting up any variables. LIP and MPI will build ``out of the box'' as well but have numerous known bugs. Most notably the author of MPI has recently stopped working on his library and LIP has long since been discontinued. GMP requires a configuration script to run and will not build out of the box. GMP and LibTomMath are still in active development and are very stable across a variety of platforms. \subsection{Choice} LibTomMath is a relatively compact, well documented, highly optimized and portable library which seems only natural for the case study of this text. Various source files from the LibTomMath project will be included within the text. However, the reader is encouraged to download their own copy of the library to actually be able to work with the library. \chapter{Getting Started} \section{Library Basics} The trick to writing any useful library of source code is to build a solid foundation and work outwards from it. First, a problem along with allowable solution parameters should be identified and analyzed. In this particular case the inability to accomodate multiple precision integers is the problem. Futhermore, the solution must be written as portable source code that is reasonably efficient across several different computer platforms. After a foundation is formed the remainder of the library can be designed and implemented in a hierarchical fashion. That is, to implement the lowest level dependencies first and work towards the most abstract functions last. For example, before implementing a modular exponentiation algorithm one would implement a modular reduction algorithm. By building outwards from a base foundation instead of using a parallel design methodology the resulting project is highly modular. Being highly modular is a desirable property of any project as it often means the resulting product has a small footprint and updates are easy to perform. Usually when I start a project I will begin with the header file. I define the data types I think I will need and prototype the initial functions that are not dependent on other functions (within the library). After I implement these base functions I prototype more dependent functions and implement them. The process repeats until I implement all of the functions I require. For example, in the case of LibTomMath I implemented functions such as mp\_init() well before I implemented mp\_mul() and even further before I implemented mp\_exptmod(). As an example as to why this design works note that the Karatsuba and Toom-Cook multipliers were written \textit{after} the dependent function mp\_exptmod() was written. Adding the new multiplication algorithms did not require changes to the mp\_exptmod() function itself and lowered the total cost of ownership (\textit{so to speak}) and of development for new algorithms. This methodology allows new algorithms to be tested in a complete framework with relative ease. \begin{center} \begin{figure}[here] \includegraphics{pics/design_process.ps} \caption{Design Flow of the First Few Original LibTomMath Functions.} \label{pic:design_process} \end{figure} \end{center} Only after the majority of the functions were in place did I pursue a less hierarchical approach to auditing and optimizing the source code. For example, one day I may audit the multipliers and the next day the polynomial basis functions. It only makes sense to begin the text with the preliminary data types and support algorithms required as well. This chapter discusses the core algorithms of the library which are the dependents for every other algorithm. \section{What is a Multiple Precision Integer?} Recall that most programming languages, in particular ISO C \cite{ISOC}, only have fixed precision data types that on their own cannot be used to represent values larger than their precision will allow. The purpose of multiple precision algorithms is to use fixed precision data types to create and manipulate multiple precision integers which may represent values that are very large. As a well known analogy, school children are taught how to form numbers larger than nine by prepending more radix ten digits. In the decimal system the largest single digit value is $9$. However, by concatenating digits together larger numbers may be represented. Newly prepended digits (\textit{to the left}) are said to be in a different power of ten column. That is, the number $123$ can be described as having a $1$ in the hundreds column, $2$ in the tens column and $3$ in the ones column. Or more formally $123 = 1 \cdot 10^2 + 2 \cdot 10^1 + 3 \cdot 10^0$. Computer based multiple precision arithmetic is essentially the same concept. Larger integers are represented by adjoining fixed precision computer words with the exception that a different radix is used. What most people probably do not think about explicitly are the various other attributes that describe a multiple precision integer. For example, the integer $154_{10}$ has two immediately obvious properties. First, the integer is positive, that is the sign of this particular integer is positive as opposed to negative. Second, the integer has three digits in its representation. There is an additional property that the integer posesses that does not concern pencil-and-paper arithmetic. The third property is how many digits placeholders are available to hold the integer. The human analogy of this third property is ensuring there is enough space on the paper to write the integer. For example, if one starts writing a large number too far to the right on a piece of paper they will have to erase it and move left. Similarly, computer algorithms must maintain strict control over memory usage to ensure that the digits of an integer will not exceed the allowed boundaries. These three properties make up what is known as a multiple precision integer or mp\_int for short. \subsection{The mp\_int Structure} \label{sec:MPINT} The mp\_int structure is the ISO C based manifestation of what represents a multiple precision integer. The ISO C standard does not provide for any such data type but it does provide for making composite data types known as structures. The following is the structure definition used within LibTomMath. \index{mp\_int} \begin{verbatim} typedef struct { int used, alloc, sign; mp_digit *dp; } mp_int; \end{verbatim} The mp\_int structure can be broken down as follows. \begin{enumerate} \item The \textbf{used} parameter denotes how many digits of the array \textbf{dp} contain the digits used to represent a given integer. The \textbf{used} count must be positive (or zero) and may not exceed the \textbf{alloc} count. \item The \textbf{alloc} parameter denotes how many digits are available in the array to use by functions before it has to increase in size. When the \textbf{used} count of a result would exceed the \textbf{alloc} count all of the algorithms will automatically increase the size of the array to accommodate the precision of the result. \item The pointer \textbf{dp} points to a dynamically allocated array of digits that represent the given multiple precision integer. It is padded with $(\textbf{alloc} - \textbf{used})$ zero digits. The array is maintained in a least significant digit order. As a pencil and paper analogy the array is organized such that the right most digits are stored first starting at the location indexed by zero\footnote{In C all arrays begin at zero.} in the array. For example, if \textbf{dp} contains $\lbrace a, b, c, \ldots \rbrace$ where \textbf{dp}$_0 = a$, \textbf{dp}$_1 = b$, \textbf{dp}$_2 = c$, $\ldots$ then it would represent the integer $a + b\beta + c\beta^2 + \ldots$ \index{MP\_ZPOS} \index{MP\_NEG} \item The \textbf{sign} parameter denotes the sign as either zero/positive (\textbf{MP\_ZPOS}) or negative (\textbf{MP\_NEG}). \end{enumerate} \subsubsection{Valid mp\_int Structures} Several rules are placed on the state of an mp\_int structure and are assumed to be followed for reasons of efficiency. The only exceptions are when the structure is passed to initialization functions such as mp\_init() and mp\_init\_copy(). \begin{enumerate} \item The value of \textbf{alloc} may not be less than one. That is \textbf{dp} always points to a previously allocated array of digits. \item The value of \textbf{used} may not exceed \textbf{alloc} and must be greater than or equal to zero. \item The value of \textbf{used} implies the digit at index $(used - 1)$ of the \textbf{dp} array is non-zero. That is, leading zero digits in the most significant positions must be trimmed. \begin{enumerate} \item Digits in the \textbf{dp} array at and above the \textbf{used} location must be zero. \end{enumerate} \item The value of \textbf{sign} must be \textbf{MP\_ZPOS} if \textbf{used} is zero; this represents the mp\_int value of zero. \end{enumerate} \section{Argument Passing} A convention of argument passing must be adopted early on in the development of any library. Making the function prototypes consistent will help eliminate many headaches in the future as the library grows to significant complexity. In LibTomMath the multiple precision integer functions accept parameters from left to right as pointers to mp\_int structures. That means that the source (input) operands are placed on the left and the destination (output) on the right. Consider the following examples. \begin{verbatim} mp_mul(&a, &b, &c); /* c = a * b */ mp_add(&a, &b, &a); /* a = a + b */ mp_sqr(&a, &b); /* b = a * a */ \end{verbatim} The left to right order is a fairly natural way to implement the functions since it lets the developer read aloud the functions and make sense of them. For example, the first function would read ``multiply a and b and store in c''. Certain libraries (\textit{LIP by Lenstra for instance}) accept parameters the other way around, to mimic the order of assignment expressions. That is, the destination (output) is on the left and arguments (inputs) are on the right. In truth, it is entirely a matter of preference. In the case of LibTomMath the convention from the MPI library has been adopted. Another very useful design consideration, provided for in LibTomMath, is whether to allow argument sources to also be a destination. For example, the second example (\textit{mp\_add}) adds $a$ to $b$ and stores in $a$. This is an important feature to implement since it allows the calling functions to cut down on the number of variables it must maintain. However, to implement this feature specific care has to be given to ensure the destination is not modified before the source is fully read. \section{Return Values} A well implemented application, no matter what its purpose, should trap as many runtime errors as possible and return them to the caller. By catching runtime errors a library can be guaranteed to prevent undefined behaviour. However, the end developer can still manage to cause a library to crash. For example, by passing an invalid pointer an application may fault by dereferencing memory not owned by the application. In the case of LibTomMath the only errors that are checked for are related to inappropriate inputs (division by zero for instance) and memory allocation errors. It will not check that the mp\_int passed to any function is valid nor will it check pointers for validity. Any function that can cause a runtime error will return an error code as an \textbf{int} data type with one of the following values. \index{MP\_OKAY} \index{MP\_VAL} \index{MP\_MEM} \begin{center} \begin{tabular}{|l|l|} \hline \textbf{Value} & \textbf{Meaning} \\ \hline \textbf{MP\_OKAY} & The function was successful \\ \hline \textbf{MP\_VAL} & One of the input value(s) was invalid \\ \hline \textbf{MP\_MEM} & The function ran out of heap memory \\ \hline \end{tabular} \end{center} When an error is detected within a function it should free any memory it allocated, often during the initialization of temporary mp\_ints, and return as soon as possible. The goal is to leave the system in the same state it was when the function was called. Error checking with this style of API is fairly simple. \begin{verbatim} int err; if ((err = mp_add(&a, &b, &c)) != MP_OKAY) { printf("Error: %s\n", mp_error_to_string(err)); exit(EXIT_FAILURE); } \end{verbatim} The GMP \cite{GMP} library uses C style \textit{signals} to flag errors which is of questionable use. Not all errors are fatal and it was not deemed ideal by the author of LibTomMath to force developers to have signal handlers for such cases. \section{Initialization and Clearing} The logical starting point when actually writing multiple precision integer functions is the initialization and clearing of the mp\_int structures. These two algorithms will be used by the majority of the higher level algorithms. Given the basic mp\_int structure an initialization routine must first allocate memory to hold the digits of the integer. Often it is optimal to allocate a sufficiently large pre-set number of digits even though the initial integer will represent zero. If only a single digit were allocated quite a few subsequent re-allocations would occur when operations are performed on the integers. There is a tradeoff between how many default digits to allocate and how many re-allocations are tolerable. Obviously allocating an excessive amount of digits initially will waste memory and become unmanageable. If the memory for the digits has been successfully allocated then the rest of the members of the structure must be initialized. Since the initial state of an mp\_int is to represent the zero integer, the allocated digits must be set to zero. The \textbf{used} count set to zero and \textbf{sign} set to \textbf{MP\_ZPOS}. \subsection{Initializing an mp\_int} An mp\_int is said to be initialized if it is set to a valid, preferably default, state such that all of the members of the structure are set to valid values. The mp\_init algorithm will perform such an action. \begin{figure}[here] \begin{center} \begin{tabular}{l} \hline Algorithm \textbf{mp\_init}. \\ \textbf{Input}. An mp\_int $a$ \\ \textbf{Output}. Allocate memory and initialize $a$ to a known valid mp\_int state. \\ \hline \\ 1. Allocate memory for \textbf{MP\_PREC} digits. \\ 2. If the allocation failed return(\textit{MP\_MEM}) \\ 3. for $n$ from $0$ to $MP\_PREC - 1$ do \\ \hspace{3mm}3.1 $a_n \leftarrow 0$\\ 4. $a.sign \leftarrow MP\_ZPOS$\\ 5. $a.used \leftarrow 0$\\ 6. $a.alloc \leftarrow MP\_PREC$\\ 7. Return(\textit{MP\_OKAY})\\ \hline \end{tabular} \end{center} \caption{Algorithm mp\_init} \end{figure} \textbf{Algorithm mp\_init.} The \textbf{MP\_PREC} name represents a constant\footnote{Defined in the ``tommath.h'' header file within LibTomMath.} used to dictate the minimum precision of allocated mp\_int integers. Ideally, it is at least equal to $32$ since for most purposes that will be more than enough. Memory for the default number of digits is allocated first. If the allocation fails the algorithm returns immediately with the \textbf{MP\_MEM} error code. If the allocation succeeds the remaining members of the mp\_int structure must be initialized to reflect the default initial state. The allocated digits are all set to zero (step three) to ensure they are in a known state. The \textbf{sign}, \textbf{used} and \textbf{alloc} are subsequently initialized to represent the zero integer. By step seven the algorithm returns a success code and the mp\_int $a$ has been successfully initialized to a valid state representing the integer zero. \textbf{Remark.} This function introduces the idiosyncrasy that all iterative loops, commonly initiated with the ``for'' keyword, iterate incrementally when the ``to'' keyword is placed between two expressions. For example, ``for $a$ from $b$ to $c$ do'' means that a subsequent expression (or body of expressions) are to be evaluated upto $c - b$ times so long as $b \le c$. In each iteration the variable $a$ is substituted for a new integer that lies inclusively between $b$ and $c$. If $b > c$ occured the loop would not iterate. By contrast if the ``downto'' keyword were used in place of ``to'' the loop would iterate decrementally. \vspace{+3mm}\begin{small} \hspace{-5.1mm}{\bf File}: bn\_mp\_init.c \vspace{-3mm} \begin{alltt} 016 017 /* init a new bigint */ 018 int mp_init (mp_int * a) 019 \{ 020 /* allocate memory required and clear it */ 021 a->dp = OPT_CAST(mp_digit) XCALLOC (sizeof (mp_digit), MP_PREC); 022 if (a->dp == NULL) \{ 023 return MP_MEM; 024 \} 025 026 /* set the used to zero, allocated digits to the default precision 027 * and sign to positive */ 028 a->used = 0; 029 a->alloc = MP_PREC; 030 a->sign = MP_ZPOS; 031 032 return MP_OKAY; 033 \} \end{alltt} \end{small} One immediate observation of this initializtion function is that it does not return a pointer to a mp\_int structure. It is assumed that the caller has already allocated memory for the mp\_int structure, typically on the application stack. The call to mp\_init() is used only to initialize the members of the structure to a known default state. Before any of the other members of the structure are initialized memory from the application heap is allocated with the calloc() function (line @22,calloc@). The size of the allocated memory is large enough to hold \textbf{MP\_PREC} mp\_digit variables. The calloc() function is used instead\footnote{calloc() will allocate memory in the same manner as malloc() except that it also sets the contents to zero upon successfully allocating the memory.} of malloc() since digits have to be set to zero for the function to finish correctly. The \textbf{OPT\_CAST} token is a macro definition which will turn into a cast from void * to mp\_digit * for C++ compilers. It is not required for C compilers. After the memory has been successfully allocated the remainder of the members are initialized (lines 28 through 30) to their respective default states. At this point the algorithm has succeeded and a success code is returned to the calling function. If this function returns \textbf{MP\_OKAY} it is safe to assume the mp\_int structure has been properly initialized and is safe to use with other functions within the library. \subsection{Clearing an mp\_int} When an mp\_int is no longer required by the application, the memory that has been allocated for its digits must be returned to the application's memory pool with the mp\_clear algorithm. \begin{figure}[here] \begin{center} \begin{tabular}{l} \hline Algorithm \textbf{mp\_clear}. \\ \textbf{Input}. An mp\_int $a$ \\ \textbf{Output}. The memory for $a$ is freed for reuse. \\ \hline \\ 1. If $a$ has been previously freed then return(\textit{MP\_OKAY}). \\ 2. for $n$ from 0 to $a.used - 1$ do \\ \hspace{3mm}2.1 $a_n \leftarrow 0$ \\ 3. Free the memory allocated for the digits of $a$. \\ 4. $a.used \leftarrow 0$ \\ 5. $a.alloc \leftarrow 0$ \\ 6. $a.sign \leftarrow MP\_ZPOS$ \\ 7. Return(\textit{MP\_OKAY}). \\ \hline \end{tabular} \end{center} \caption{Algorithm mp\_clear} \end{figure} \textbf{Algorithm mp\_clear.} This algorithm releases the memory allocated for an mp\_int back into the memory pool for reuse. It is designed such that a given mp\_int structure can be cleared multiple times between initializations without attempting to free the memory twice\footnote{In ISO C for example, calling free() twice on the same memory block causes undefinied behaviour.}. The first step determines if the mp\_int structure has been marked as free already. If it has, the algorithm returns success immediately as no further actions are required. Otherwise, the algorithm will proceed to put the structure in a known empty and otherwise invalid state. First the digits of the mp\_int are set to zero. The memory that has been allocated for the digits is then freed. The \textbf{used} and \textbf{alloc} counts are both set to zero and the \textbf{sign} set to \textbf{MP\_ZPOS}. This known fixed state for cleared mp\_int structures will make debuging easier for the end developer. That is, if they spot (via their debugger) an mp\_int they are using that is in this state it will be obvious that they erroneously and prematurely cleared the mp\_int structure. Note that once an mp\_int has been cleared the mp\_int structure is no longer in a valid state for any other algorithm with the exception of algorithms mp\_init, mp\_init\_copy, mp\_init\_size and mp\_clear. \vspace{+3mm}\begin{small} \hspace{-5.1mm}{\bf File}: bn\_mp\_clear.c \vspace{-3mm} \begin{alltt} 016 017 /* clear one (frees) */ 018 void 019 mp_clear (mp_int * a) 020 \{ 021 /* only do anything if a hasn't been freed previously */ 022 if (a->dp != NULL) \{ 023 /* first zero the digits */ 024 memset (a->dp, 0, sizeof (mp_digit) * a->used); 025 026 /* free ram */ 027 XFREE(a->dp); 028 029 /* reset members to make debugging easier */ 030 a->dp = NULL; 031 a->alloc = a->used = 0; 032 a->sign = MP_ZPOS; 033 \} 034 \} \end{alltt} \end{small} The ``if'' statement (line 22) prevents the heap from being corrupted if a user double-frees an mp\_int. This is because once the memory is freed the pointer is set to \textbf{NULL} (line 30). Without the check, code that accidentally calls mp\_clear twice for a given mp\_int structure would try to free the memory allocated for the digits twice. This may cause some C libraries to signal a fault. By setting the pointer to \textbf{NULL} it helps debug code that may inadvertently free the mp\_int before it is truly not needed, because attempts to reference digits should fail immediately. The allocated digits are set to zero before being freed (line 24). This is ideal for cryptographic situations where the integer that the mp\_int represents might need to be kept a secret. \section{Maintenance Algorithms} The previous sections describes how to initialize and clear an mp\_int structure. To further support operations that are to be performed on mp\_int structures (such as addition and multiplication) the dependent algorithms must be able to augment the precision of an mp\_int and initialize mp\_ints with differing initial conditions. These algorithms complete the set of low level algorithms required to work with mp\_int structures in the higher level algorithms such as addition, multiplication and modular exponentiation. \subsection{Augmenting an mp\_int's Precision} When storing a value in an mp\_int structure, a sufficient number of digits must be available to accomodate the entire result of an operation without loss of precision. Quite often the size of the array given by the \textbf{alloc} member is large enough to simply increase the \textbf{used} digit count. However, when the size of the array is too small it must be re-sized appropriately to accomodate the result. The mp\_grow algorithm will provide this functionality. \newpage\begin{figure}[here] \begin{center} \begin{tabular}{l} \hline Algorithm \textbf{mp\_grow}. \\ \textbf{Input}. An mp\_int $a$ and an integer $b$. \\ \textbf{Output}. $a$ is expanded to accomodate $b$ digits. \\ \hline \\ 1. if $a.alloc \ge b$ then return(\textit{MP\_OKAY}) \\ 2. $u \leftarrow b\mbox{ (mod }MP\_PREC\mbox{)}$ \\ 3. $v \leftarrow b + 2 \cdot MP\_PREC - u$ \\ 4. Re-Allocate the array of digits $a$ to size $v$ \\ 5. If the allocation failed then return(\textit{MP\_MEM}). \\ 6. for n from a.alloc to $v - 1$ do \\ \hspace{+3mm}6.1 $a_n \leftarrow 0$ \\ 7. $a.alloc \leftarrow v$ \\ 8. Return(\textit{MP\_OKAY}) \\ \hline \end{tabular} \end{center} \caption{Algorithm mp\_grow} \end{figure} \textbf{Algorithm mp\_grow.} It is ideal to prevent re-allocations from being performed if they are not required (step one). This is useful to prevent mp\_ints from growing excessively in code that erroneously calls mp\_grow. The requested digit count is padded up to next multiple of \textbf{MP\_PREC} plus an additional \textbf{MP\_PREC} (steps two and three). This helps prevent many trivial reallocations that would grow an mp\_int by trivially small values. It is assumed that the reallocation (step four) leaves the lower $a.alloc$ digits of the mp\_int intact. This is much akin to how the \textit{realloc} function from the standard C library works. Since the newly allocated digits are assumed to contain undefined values they are initially set to zero. \vspace{+3mm}\begin{small} \hspace{-5.1mm}{\bf File}: bn\_mp\_grow.c \vspace{-3mm} \begin{alltt} 016 017 /* grow as required */ 018 int mp_grow (mp_int * a, int size) 019 \{ 020 int i; 021 mp_digit *tmp; 022 023 /* if the alloc size is smaller alloc more ram */ 024 if (a->alloc < size) \{ 025 /* ensure there are always at least MP_PREC digits extra on top */ 026 size += (MP_PREC * 2) - (size % MP_PREC); 027 028 /* reallocate the array a->dp 029 * 030 * We store the return in a temporary variable 031 * in case the operation failed we don't want 032 * to overwrite the dp member of a. 033 */ 034 tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * size); 035 if (tmp == NULL) \{ 036 /* reallocation failed but "a" is still valid [can be freed] */ 037 return MP_MEM; 038 \} 039 040 /* reallocation succeeded so set a->dp */ 041 a->dp = tmp; 042 043 /* zero excess digits */ 044 i = a->alloc; 045 a->alloc = size; 046 for (; i < a->alloc; i++) \{ 047 a->dp[i] = 0; 048 \} 049 \} 050 return MP_OKAY; 051 \} \end{alltt} \end{small} The first step is to see if we actually need to perform a re-allocation at all (line 24). If a reallocation must occur the digit count is padded upwards to help prevent many trivial reallocations (line 26). Next the reallocation is performed and the return of realloc() is stored in a temporary pointer named $tmp$ (line 36). The return is stored in a temporary instead of $a.dp$ to prevent the code from losing the original pointer in case the reallocation fails. Had the return been stored in $a.dp$ instead there would be no way to reclaim the heap originally used. If the reallocation fails the function will return \textbf{MP\_MEM} (line 37), otherwise, the value of $tmp$ is assigned to the pointer $a.dp$ and the function continues. A simple for loop from line 46 to line 51 will zero all digits that were above the old \textbf{alloc} limit to make sure the integer is in a known state. \subsection{Initializing Variable Precision mp\_ints} Occasionally the number of digits required will be known in advance of an initialization, based on, for example, the size of input mp\_ints to a given algorithm. The purpose of algorithm mp\_init\_size is similar to mp\_init except that it will allocate \textit{at least} a specified number of digits. \begin{figure}[here] \begin{small} \begin{center} \begin{tabular}{l} \hline Algorithm \textbf{mp\_init\_size}. \\ \textbf{Input}. An mp\_int $a$ and the requested number of digits $b$. \\ \textbf{Output}. $a$ is initialized to hold at least $b$ digits. \\ \hline \\ 1. $u \leftarrow b \mbox{ (mod }MP\_PREC\mbox{)}$ \\ 2. $v \leftarrow b + 2 \cdot MP\_PREC - u$ \\ 3. Allocate $v$ digits. \\ 4. for $n$ from $0$ to $v - 1$ do \\ \hspace{3mm}4.1 $a_n \leftarrow 0$ \\ 5. $a.sign \leftarrow MP\_ZPOS$\\ 6. $a.used \leftarrow 0$\\ 7. $a.alloc \leftarrow v$\\ 8. Return(\textit{MP\_OKAY})\\ \hline \end{tabular} \end{center} \end{small} \caption{Algorithm mp\_init\_size} \end{figure} \textbf{Algorithm mp\_init\_size.} This algorithm will initialize an mp\_int structure $a$ like algorithm mp\_init with the exception that the number of digits allocated can be controlled by the second input argument $b$. The input size is padded upwards so it is a multiple of \textbf{MP\_PREC} plus an additional \textbf{MP\_PREC} digits. This padding is used to prevent trivial allocations from becoming a bottleneck in the rest of the algorithms. Like algorithm mp\_init, the mp\_int structure is initialized to a default state representing the integer zero. This particular algorithm is useful if it is known ahead of time the approximate size of the input. If the approximation is correct no further memory re-allocations are required to work with the mp\_int. \vspace{+3mm}\begin{small} \hspace{-5.1mm}{\bf File}: bn\_mp\_init\_size.c \vspace{-3mm} \begin{alltt} 016 017 /* init an mp_init for a given size */ 018 int mp_init_size (mp_int * a, int size) 019 \{ 020 /* pad size so there are always extra digits */ 021 size += (MP_PREC * 2) - (size % MP_PREC); 022 023 /* alloc mem */ 024 a->dp = OPT_CAST(mp_digit) XCALLOC (sizeof (mp_digit), size); 025 if (a->dp == NULL) \{ 026 return MP_MEM; 027 \} 028 a->used = 0; 029 a->alloc = size; 030 a->sign = MP_ZPOS; 031 032 return MP_OKAY; 033 \} \end{alltt} \end{small} The number of digits $b$ requested is padded (line 21) by first augmenting it to the next multiple of \textbf{MP\_PREC} and then adding \textbf{MP\_PREC} to the result. If the memory can be successfully allocated the mp\_int is placed in a default state representing the integer zero. Otherwise, the error code \textbf{MP\_MEM} will be returned (line 26). The digits are allocated and set to zero at the same time with the calloc() function (line @25,calloc@). The \textbf{used} count is set to zero, the \textbf{alloc} count set to the padded digit count and the \textbf{sign} flag set to \textbf{MP\_ZPOS} to achieve a default valid mp\_int state (lines 28, 29 and 30). If the function returns succesfully then it is correct to assume that the mp\_int structure is in a valid state for the remainder of the functions to work with. \subsection{Multiple Integer Initializations and Clearings} Occasionally a function will require a series of mp\_int data types to be made available simultaneously. The purpose of algorithm mp\_init\_multi is to initialize a variable length array of mp\_int structures in a single statement. It is essentially a shortcut to multiple initializations. \newpage\begin{figure}[here] \begin{center} \begin{tabular}{l} \hline Algorithm \textbf{mp\_init\_multi}. \\ \textbf{Input}. Variable length array $V_k$ of mp\_int variables of length $k$. \\ \textbf{Output}. The array is initialized such that each mp\_int of $V_k$ is ready to use. \\ \hline \\ 1. for $n$ from 0 to $k - 1$ do \\ \hspace{+3mm}1.1. Initialize the mp\_int $V_n$ (\textit{mp\_init}) \\ \hspace{+3mm}1.2. If initialization failed then do \\ \hspace{+6mm}1.2.1. for $j$ from $0$ to $n$ do \\ \hspace{+9mm}1.2.1.1. Free the mp\_int $V_j$ (\textit{mp\_clear}) \\ \hspace{+6mm}1.2.2. Return(\textit{MP\_MEM}) \\ 2. Return(\textit{MP\_OKAY}) \\ \hline \end{tabular} \end{center} \caption{Algorithm mp\_init\_multi} \end{figure} \textbf{Algorithm mp\_init\_multi.} The algorithm will initialize the array of mp\_int variables one at a time. If a runtime error has been detected (\textit{step 1.2}) all of the previously initialized variables are cleared. The goal is an ``all or nothing'' initialization which allows for quick recovery from runtime errors. \vspace{+3mm}\begin{small} \hspace{-5.1mm}{\bf File}: bn\_mp\_init\_multi.c \vspace{-3mm} \begin{alltt} 016 #include <stdarg.h> 017 018 int mp_init_multi(mp_int *mp, ...) 019 \{ 020 mp_err res = MP_OKAY; /* Assume ok until proven otherwise */ 021 int n = 0; /* Number of ok inits */ 022 mp_int* cur_arg = mp; 023 va_list args; 024 025 va_start(args, mp); /* init args to next argument from caller */ 026 while (cur_arg != NULL) \{ 027 if (mp_init(cur_arg) != MP_OKAY) \{ 028 /* Oops - error! Back-track and mp_clear what we already 029 succeeded in init-ing, then return error. 030 */ 031 va_list clean_args; 032 033 /* end the current list */ 034 va_end(args); 035 036 /* now start cleaning up */ 037 cur_arg = mp; 038 va_start(clean_args, mp); 039 while (n--) \{ 040 mp_clear(cur_arg); 041 cur_arg = va_arg(clean_args, mp_int*); 042 \} 043 va_end(clean_args); 044 res = MP_MEM; 045 break; 046 \} 047 n++; 048 cur_arg = va_arg(args, mp_int*); 049 \} 050 va_end(args); 051 return res; /* Assumed ok, if error flagged above. */ 052 \} 053 \end{alltt} \end{small} This function intializes a variable length list of mp\_int structure pointers. However, instead of having the mp\_int structures in an actual C array they are simply passed as arguments to the function. This function makes use of the ``...'' argument syntax of the C programming language. The list is terminated with a final \textbf{NULL} argument appended on the right. The function uses the ``stdarg.h'' \textit{va} functions to step portably through the arguments to the function. A count $n$ of succesfully initialized mp\_int structures is maintained (line 47) such that if a failure does occur, the algorithm can backtrack and free the previously initialized structures (lines 27 to 46). \subsection{Clamping Excess Digits} When a function anticipates a result will be $n$ digits it is simpler to assume this is true within the body of the function instead of checking during the computation. For example, a multiplication of a $i$ digit number by a $j$ digit produces a result of at most $i + j$ digits. It is entirely possible that the result is $i + j - 1$ though, with no final carry into the last position. However, suppose the destination had to be first expanded (\textit{via mp\_grow}) to accomodate $i + j - 1$ digits than further expanded to accomodate the final carry. That would be a considerable waste of time since heap operations are relatively slow. The ideal solution is to always assume the result is $i + j$ and fix up the \textbf{used} count after the function terminates. This way a single heap operation (\textit{at most}) is required. However, if the result was not checked there would be an excess high order zero digit. For example, suppose the product of two integers was $x_n = (0x_{n-1}x_{n-2}...x_0)_{\beta}$. The leading zero digit will not contribute to the precision of the result. In fact, through subsequent operations more leading zero digits would accumulate to the point the size of the integer would be prohibitive. As a result even though the precision is very low the representation is excessively large. The mp\_clamp algorithm is designed to solve this very problem. It will trim high-order zeros by decrementing the \textbf{used} count until a non-zero most significant digit is found. Also in this system, zero is considered to be a positive number which means that if the \textbf{used} count is decremented to zero, the sign must be set to \textbf{MP\_ZPOS}. \begin{figure}[here] \begin{center} \begin{tabular}{l} \hline Algorithm \textbf{mp\_clamp}. \\ \textbf{Input}. An mp\_int $a$ \\ \textbf{Output}. Any excess leading zero digits of $a$ are removed \\ \hline \\ 1. while $a.used > 0$ and $a_{a.used - 1} = 0$ do \\ \hspace{+3mm}1.1 $a.used \leftarrow a.used - 1$ \\ 2. if $a.used = 0$ then do \\ \hspace{+3mm}2.1 $a.sign \leftarrow MP\_ZPOS$ \\ \hline \\ \end{tabular} \end{center} \caption{Algorithm mp\_clamp} \end{figure} \textbf{Algorithm mp\_clamp.} As can be expected this algorithm is very simple. The loop on step one is expected to iterate only once or twice at the most. For example, this will happen in cases where there is not a carry to fill the last position. Step two fixes the sign for when all of the digits are zero to ensure that the mp\_int is valid at all times. \vspace{+3mm}\begin{small} \hspace{-5.1mm}{\bf File}: bn\_mp\_clamp.c \vspace{-3mm} \begin{alltt} 016 017 /* trim unused digits 018 * 019 * This is used to ensure that leading zero digits are 020 * trimed and the leading "used" digit will be non-zero 021 * Typically very fast. Also fixes the sign if there 022 * are no more leading digits 023 */ 024 void 025 mp_clamp (mp_int * a) 026 \{ 027 /* decrease used while the most significant digit is 028 * zero. 029 */ 030 while (a->used > 0 && a->dp[a->used - 1] == 0) \{ 031 --(a->used); 032 \} 033 034 /* reset the sign flag if used == 0 */ 035 if (a->used == 0) \{ 036 a->sign = MP_ZPOS; 037 \} 038 \} \end{alltt} \end{small} Note on line 27 how to test for the \textbf{used} count is made on the left of the \&\& operator. In the C programming language the terms to \&\& are evaluated left to right with a boolean short-circuit if any condition fails. This is important since if the \textbf{used} is zero the test on the right would fetch below the array. That is obviously undesirable. The parenthesis on line 30 is used to make sure the \textbf{used} count is decremented and not the pointer ``a''. \section*{Exercises} \begin{tabular}{cl} $\left [ 1 \right ]$ & Discuss the relevance of the \textbf{used} member of the mp\_int structure. \\ & \\ $\left [ 1 \right ]$ & Discuss the consequences of not using padding when performing allocations. \\ & \\ $\left [ 2 \right ]$ & Estimate an ideal value for \textbf{MP\_PREC} when performing 1024-bit RSA \\ & encryption when $\beta = 2^{28}$. \\ & \\ $\left [ 1 \right ]$ & Discuss the relevance of the algorithm mp\_clamp. What does it prevent? \\ & \\ $\left [ 1 \right ]$ & Give an example of when the algorithm mp\_init\_copy might be useful. \\ & \\ \end{tabular} %%% % CHAPTER FOUR %%% \chapter{Basic Operations} \section{Introduction} In the previous chapter a series of low level algorithms were established that dealt with initializing and maintaining mp\_int structures. This chapter will discuss another set of seemingly non-algebraic algorithms which will form the low level basis of the entire library. While these algorithm are relatively trivial it is important to understand how they work before proceeding since these algorithms will be used almost intrinsically in the following chapters. The algorithms in this chapter deal primarily with more ``programmer'' related tasks such as creating copies of mp\_int structures, assigning small values to mp\_int structures and comparisons of the values mp\_int structures represent. \section{Assigning Values to mp\_int Structures} \subsection{Copying an mp\_int} Assigning the value that a given mp\_int structure represents to another mp\_int structure shall be known as making a copy for the purposes of this text. The copy of the mp\_int will be a separate entity that represents the same value as the mp\_int it was copied from. The mp\_copy algorithm provides this functionality. \newpage\begin{figure}[here] \begin{center} \begin{tabular}{l} \hline Algorithm \textbf{mp\_copy}. \\ \textbf{Input}. An mp\_int $a$ and $b$. \\ \textbf{Output}. Store a copy of $a$ in $b$. \\ \hline \\ 1. If $b.alloc < a.used$ then grow $b$ to $a.used$ digits. (\textit{mp\_grow}) \\ 2. for $n$ from 0 to $a.used - 1$ do \\ \hspace{3mm}2.1 $b_{n} \leftarrow a_{n}$ \\ 3. for $n$ from $a.used$ to $b.used - 1$ do \\ \hspace{3mm}3.1 $b_{n} \leftarrow 0$ \\ 4. $b.used \leftarrow a.used$ \\ 5. $b.sign \leftarrow a.sign$ \\ 6. return(\textit{MP\_OKAY}) \\ \hline \end{tabular} \end{center} \caption{Algorithm mp\_copy} \end{figure} \textbf{Algorithm mp\_copy.} This algorithm copies the mp\_int $a$ such that upon succesful termination of the algorithm the mp\_int $b$ will represent the same integer as the mp\_int $a$. The mp\_int $b$ shall be a complete and distinct copy of the mp\_int $a$ meaing that the mp\_int $a$ can be modified and it shall not affect the value of the mp\_int $b$. If $b$ does not have enough room for the digits of $a$ it must first have its precision augmented via the mp\_grow algorithm. The digits of $a$ are copied over the digits of $b$ and any excess digits of $b$ are set to zero (step two and three). The \textbf{used} and \textbf{sign} members of $a$ are finally copied over the respective members of $b$. \textbf{Remark.} This algorithm also introduces a new idiosyncrasy that will be used throughout the rest of the text. The error return codes of other algorithms are not explicitly checked in the pseudo-code presented. For example, in step one of the mp\_copy algorithm the return of mp\_grow is not explicitly checked to ensure it succeeded. Text space is limited so it is assumed that if a algorithm fails it will clear all temporarily allocated mp\_ints and return the error code itself. However, the C code presented will demonstrate all of the error handling logic required to implement the pseudo-code. \vspace{+3mm}\begin{small} \hspace{-5.1mm}{\bf File}: bn\_mp\_copy.c \vspace{-3mm} \begin{alltt} 016 017 /* copy, b = a */ 018 int 019 mp_copy (mp_int * a, mp_int * b) 020 \{ 021 int res, n; 022 023 /* if dst == src do nothing */ 024 if (a == b) \{ 025 return MP_OKAY; 026 \} 027 028 /* grow dest */ 029 if (b->alloc < a->used) \{ 030 if ((res = mp_grow (b, a->used)) != MP_OKAY) \{ 031 return res; 032 \} 033 \} 034 035 /* zero b and copy the parameters over */ 036 \{ 037 register mp_digit *tmpa, *tmpb; 038 039 /* pointer aliases */ 040 041 /* source */ 042 tmpa = a->dp; 043 044 /* destination */ 045 tmpb = b->dp; 046 047 /* copy all the digits */ 048 for (n = 0; n < a->used; n++) \{ 049 *tmpb++ = *tmpa++; 050 \} 051 052 /* clear high digits */ 053 for (; n < b->used; n++) \{ 054 *tmpb++ = 0; 055 \} 056 \} 057 058 /* copy used count and sign */ 059 b->used = a->used; 060 b->sign = a->sign; 061 return MP_OKAY; 062 \} \end{alltt} \end{small} Occasionally a dependent algorithm may copy an mp\_int effectively into itself such as when the input and output mp\_int structures passed to a function are one and the same. For this case it is optimal to return immediately without copying digits (line 24). The mp\_int $b$ must have enough digits to accomodate the used digits of the mp\_int $a$. If $b.alloc$ is less than $a.used$ the algorithm mp\_grow is used to augment the precision of $b$ (lines 29 to 33). In order to simplify the inner loop that copies the digits from $a$ to $b$, two aliases $tmpa$ and $tmpb$ point directly at the digits of the mp\_ints $a$ and $b$ respectively. These aliases (lines 42 and 45) allow the compiler to access the digits without first dereferencing the mp\_int pointers and then subsequently the pointer to the digits. After the aliases are established the digits from $a$ are copied into $b$ (lines 48 to 50) and then the excess digits of $b$ are set to zero (lines 53 to 55). Both ``for'' loops make use of the pointer aliases and in fact the alias for $b$ is carried through into the second ``for'' loop to clear the excess digits. This optimization allows the alias to stay in a machine register fairly easy between the two loops. \textbf{Remarks.} The use of pointer aliases is an implementation methodology first introduced in this function that will be used considerably in other functions. Technically, a pointer alias is simply a short hand alias used to lower the number of pointer dereferencing operations required to access data. For example, a for loop may resemble \begin{alltt} for (x = 0; x < 100; x++) \{ a->num[4]->dp[x] = 0; \} \end{