Multiprecision Arithmetic Part IV

by John Sullivan

mathematics and statistics programming techniques

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 73

Applications of Multiprecision Arithmetic¶

In parts 1 & 2 of this series, I introduced the basic functions for performing multiprecision arithmetic. But why do we need to perform arithmetic to such precision? The day-to-day operations of your bank account, or even of the world’s money markets, require no more precision than that supplied by APL for its ordinary arithmetic operations. Scientific measurements no matter how small or large are usually accompanied by some (small or large) error, so accuracy to more than 4 significant figures is usually unnecessary.

In an article in Scientific American[1], Fred Gruenberger points out that in using arithmetic with only a small number of significant digits, repeated squaring of a number only slightly greater than 1 will soon result in inaccurate results. He then gives the result of squaring 1.0000001 seven times (equivalent to raising it to the power of 128), which results in 1.0000128000812803413… (896 decimal places). Try it with the multiprecision functions. Even on my 486/33 the calculation is quick. Set ⎕PP to 15 and enter 1.0000001*128, and you will see the problem. The result, 1.00001280008129, is wrong in the last displayed place. And as the power gets larger the results become less accurate. Don’t think that this is a silly exercise: some problems involve iterated multiplication like this; can you really trust the results your computer gives you?

Next, Gruenberger discusses the continuing search for larger prime numbers. Before the development of computers, the largest known prime was \(2^{127} - 1\), which has only 39 digits. This activity requires the multiprecision arithmetic of parts 1 and 2 of this series, and the primality testing and factoring of large numbers from part 3. But, having discovered a few 100-digit primes, and found the factors of a few 100-digit non-primes, you may be wondering what sort of practical application there is for these numbers.

Gruenberger mentions some other problems which require high precision arithmetic, such as the solution of simultaneous equations where the determinant of the associated matrix is close to zero, and a couple of mathematical recreations which lead to very large numbers. However, he makes no mention of what is considered by many to be the only commercial application of number theory in the modern world: data encryption and decryption.

This algorithm for public-key cryptography was described in the late 1970s by Rivest, Shamir and Adleman[2], and it is called the RSA system after them. Although all of the number theory behind it has been known for the last couple of centuries, it is the advent of the digital computer that has made the calculations practical. Two numbers, the encryption key \(e\) and the modulus \(n\) are published, and anybody can use them to encode a message, which they then send to you. The idea is that \(n\) cannot be factored without a lot of effort, so the code-breaker cannot easily calculate the encryption key \(d\) to decipher the coded message. There are many books and articles on the RSA algorithm: I obtained my inspiration from chapter 6 of Riesel[3], although I have not used his examples.

Some number theory¶

Before we can go merrily encrypting and decrypting our secret text, we need some number theory. I am not going to prove any of the following statements: look in any good textbook on number theory, such as those mentioned at the end of part 3. These statements are the key to what follows. All of the numbers concerned are, of course, integers.

If \(\operatorname{GCD}(a,n)=1\) then there exists a number \(b\) with \(1 \le b < n\) such that \(ab \equiv 1 \bmod n\). For example \(2 \times 4 \equiv 1 \bmod 7\). We call \(b\) the multiplicative inverse of \(a\) modulo \(n\).

For any integer \(a\) with \(\operatorname{GCD}(a,n)=1\), \(a^{\phi(n)} \equiv 1 \bmod n\), where \(\phi(n)\) is the number of positive integers not exceeding \(n\) that are relatively prime to \(n\) (\(\phi(n)\) is known as Euler’s totient function). Not only that, but if \(n\) is squarefree \(a^{k\phi(n)+1} \equiv a \bmod n\) for all integers \(a\) and \(k\). If \(n\) is prime then \(\phi(n) = n - 1\). If \(\operatorname{GCD}(b,c)=1\) then \(\phi(bc) = \phi(b)\phi(c)\). However, thanks to a theorem due to Carmichael we can do better than this. There is a function, called Carmichael’s function, denoted by \(\lambda(n)\), such that \(\lambda(n)\) divides \(\phi(n)\), and \(a^{\lambda(n)} \equiv 1 \bmod n\), with \(\operatorname{GCD}(a,n)=1\) and \(a^{k\lambda(n)+1} \equiv a \bmod n\) for all integers \(a\) and \(k\) and squarefree \(n\). In general, to calculate Carmichael’s function for a composite number one requires all the prime factors of its argument, which is a beyond the scope of this article; all we need here is that if \(p\) and \(q\) are primes then \(\lambda(pq) = \operatorname{LCM}(p - 1, q - 1)\).

The algorithm¶

The algorithm for what we are about to do is simple. As mentioned above, we need three numbers: the encrypting key, the decrypting key, and the modulus. We represent our plain text as a number by some means, raise it to the power of the encrypting key modulo the modulus, and send the result to the recipient of our secret message. The recipient then raises this number to the power of the decrypting key modulo the modulus, and the result of this is the original message. To make this work, we choose two primes \(p\) and \(q\), and multiply them together to get the modulus \(n\). Choose an encrypting key \(e\) such that it is relatively prime to \(\lambda(n)\). Then the decrypting key \(d\) is just the multiplicative inverse of \(e\) modulo \(\lambda(n)\). (Some books suggest using \(\phi(n)\) here, instead of \(\lambda(n)\), but this just means raising to a larger power and doing more work for the same result.) Effectively, since \(de \equiv 1 \bmod \lambda(n)\), then \(de = k\lambda(n)\), for some \(k\), and we are raising our original number to the power of \(k\lambda(n) + 1 \bmod n\), which, as mentioned above, restores the original number.

The multiplicative inverse function¶

This is put in namespace mp along with the other multiprecision building blocks. It is based on the Euclidean Algorithm: if \(\operatorname{GCD}(a,b)=1\) then there exist \(m,n\) such that \(am + bn = 1\). This can be rearranged as \(am = 1 - bn\), which immediately gives \(am \equiv 1 \bmod n\), thus \(m\) is the multiplicative inverse of \(a\) modulo \(n\). If the GCD is not 1 then this function returns 0, which is an impossible value for the inverse.

As with the GCD function (see part 3) this function is divided into two paths, depending on whether we are processing scalars or multiprecision numbers. The first thing we do is find the GCD of the input (I use \(x\) here to represent \(a\), and \(n\) is the modulus). The Euclidean algorithm says that when we start with \(a\) and \(b\), we obtain the following system of equations:

\[\begin{align*} a &= q_1b + r_1,\ 0 < r_1 < b \\ b &= q_2r_1 + r_2,\ 0 < r_2 < r_1 \\ r_1 &= q_3r_2 + r_3,\ 0 < r_3 < r_2 \\ &\;\;\ldots \\ r_{n-2} &= q_n\,r_{n-1} + r_n,\ 0 < r_n < r_{n-1} \\ r_{n-1} &= q_{n+1}\,r_n + 0 \end{align*}\]

We start off by generating all of these quotients and remainders, storing the quotients in vector \(a\):

    ∇ z←x Inv_mod n;a;b;N;rn;r0;r1
[1]    a←⍬ ⋄ N←n
[2]    →(1∨.≠(⍴x),⍴n)/∆10
[3]    z←0 ⋄ x←N|x
[4]   ∆1:a←a,⊃b←(0,x)⊤n
[5]    n←x ⋄ →(0≠x←1⊃b)/∆1

At this point \(n\) holds the GCD of the original \(n\) and \(x\). If this is not 1 then there is no multiplicative inverse

[6]    →(1≠n)/0

Now we have to go back recursively. We calculate \(r_n\) in terms of \(r_0\) (\(= b\)) and \(r_1\). The recurrence relation we calculate is

\[\begin{equation*} r_i = r_{i-2} - q_i\,r_{i-1} \end{equation*}\]

[7]    r0←1 0 ⋄ r1←0 1
[8]   ∆3:rn←r0-r1×⊃a
[9]    r0←r1 ⋄ r1←rn ⋄ a←1↓a
[10]   →(1<⍴a)/∆3

And our result is the coefficient of \(r_1\) modulo \(n\).

[11]   z←N|1⊃rn ⋄ →0

And now for the multiprecision version. This part of the function is the same as the scalar version, with minor differences to allow for the way the multiprecision functions have been coded.

[12]  ∆10:z←0 0 ⋄ x←x Imod N
[13]  ∆11:a←a,(b←n Idiv x)[0]
[14]   n←x ⋄ →(~0 0 Fequal x←1⊃b)/∆11
[15]   →(~0 1 Fequal n)/0
[16]   r0←(0 1)(0 0) ⋄ r1←(0 0)(0 1)
[17]  ∆13:rn←r0 Fsub¨r1 Fmul¨a[0]
[18]   r0←r1 ⋄ r1←rn ⋄ a←1↓a
[19]   →(1<⍴a)/∆13
[20]   z←(1⊃rn)Imod N
    ∇

Testing the Method¶

For demonstration purposes it is not necessary to search out abstruse prime numbers. At the back of Riesel[3] there are tables of large prime numbers of various types, and I have chosen some of those. You will notice that when using the factoring function I showed in part 3, the factors of \(n\) fall out in no time at all. This is O.K. for demonstrations, but is not supposed to happen in reality!

\[\begin{align*} p &= 9 \times 2^{63} - 1 = 83010348331692982271 \\ q &= 9 \times 2^{63} + 1 = 83010348331692982273 \\ e &= 27 \times 2^{70} - 1 = 31875973759370105192447 \end{align*}\]

Since \(p\) and \(q\) are primes, \(\lambda(pq)\) is just \(\operatorname{LCM}(p - 1, q - 1)\). From this we can calculate the decrypting key:

\[\begin{equation*} d = 2591880212318106415703971813838725054463 \end{equation*}\]

I created another namespace to hide these away, called rsa. It seemed a good idea to use a setup function to generate these numbers:

    ∇ setup;d;e;n;p;q;s
[1]    s←⎕CS'#.mp'
[2]    p←0 ¯1 Fadd 0 9 Fmul 0 2 Fspow 63
[3]    q←p Fadd 0 2
[4]    n←p Fmul q
[5]    e←0 ¯1 Fadd 0 27 Fmul 0 2 Fspow 70
[6]    d←e Inv_mod(0 ¯1 Fadd p)Lcm 0 ¯1 Fadd q
[7]    ⎕CS s
[8]    D←d ⋄ E←e ⋄ N←n ⋄ P←p ⋄ Q←q
    ∇

The LCM function (in namespace mp) is very simple. As with other functions it has two paths, one for scalars and one for multiprecision numbers.

    ∇ z←a Lcm b
[1]    →(1∧.=(⍴a),⍴b)/∆1
[2]    z←⊃(a Fmul b)Idiv a Gcd b
[3]    →0
[4]   ∆1:z←(a×b)÷a Gcd b
    ∇

In order to test the method, I made a function called test. It is sufficient to use a single number, smaller than \(n\), to generate an encrypted version, and to decrypt it again. If the decrypted version equals the original version then the process works.

    ∇ test;M;C
[1]    M←#.mp.Fexec'1234567890123456789012345678901234567890'
[2]    'Original  Message is ',0 #.mp.Ffmt M
[3]    'Encrypted Message is ',0 #.mp.Ffmt C←M #.mp.Impow E N
[4]    'Decrypted Message is ',0 #.mp.Ffmt C #.mp.Impow D N
    ∇

And the following output should surprise nobody.

      test
Original  Message is 1234567890123456789012345678901234567890
Encrypted Message is 6318195541048362681928930668568439899876
Decrypted Message is 1234567890123456789012345678901234567890

Practical Matters¶

How you implement this algorithm is very much up to you and your correspondents. Here is one method that I devised for testing purposes. Again, this is different from the examples in Riesel.

First of all we determine what characters we are going to allow in our messages. I added this line to the end of the setup function, above.

Alf←' ',⎕D,⎕A,'.,:;!?()'

Now you need your normal capitalization function, something like this

    ∇ x←Caps x;a
[1]    a←⎕AV ⋄ a[⎕AV⍳'abcdefghijklmnopqrstuvwxyz']←⎕A
[2]    x←a[⎕AV⍳x]
    ∇

Convert your plaintext to caps, and replace characters not in the "alphabet" with blanks. Generate a set of origin-zero indexes into this alphabet, and treat them as a radix-⍴Alf multiprecision number by catenating a zero at the start (this is a beneficial side-effect of the variable-radix design of the multiprecision suite). Since the first character in our alphabet is a blank this has the effect of stripping leading blanks from the message, which is probably a good idea. Convert this number to a multiprecision number in our default radix, then use encode (⊤) on it to ensure that we break it into parts that are not greater than the modulus \(n\) of our encryption algorithm. Each of the resulting numbers is then encrypted, and the reverse process is carried out on the result, although this time we just use an alphabet of the usual capital letters. Then we split the output into the 5-character blocks, so beloved of espionage systems in the 30’s to early 50’s.

    ∇ z←Encrypt x;a;b
[1]    x←,Caps x
[2]    x[(~x∊Alf)/⍳⍴x]←' '
[3]    z←(#.mp.base,⍴Alf)#.mp.chbase 0,Alf⍳x
[4]    z←N #.mp.enc z
[5]    z←z #.mp.Impow¨⊂E N
[6]    z←N #.mp.dec z
[7]    z←1↓(26,#.mp.base)#.mp.chbase z
[8]    z←⎕A[z]
[9]    z←(5×b←⌈0.2×⍴z)↑z
[10]   z←((6×b)⍴1 1 1 1 1 0)\z
    ∇

The decryption process is just the reverse of this. First we eliminate all non-alphabetical characters from the ciphertext, then we compare it with the ciphertext without the blanks. Obviously, if these are not the same we are working on the wrong text, so we signal an error. Then the process continues as for Encrypt, except that we do not format in 5-character blocks at the end.

    ∇ z←Decrypt x;a;b
[1]    z←x∩⎕A
[2]    ⎕SIGNAL(z≢x~' ')/11
[3]    z←(#.mp.base,26)#.mp.chbase 0,⎕A⍳z
[4]    z←N #.mp.enc z
[5]    z←z #.mp.Impow¨⊂D N
[6]    z←N #.mp.dec z
[7]    z←1↓((⍴Alf),#.mp.base)#.mp.chbase z
[8]    z←Alf[z]
    ∇

And finally, here are the functions dec and enc from namespace mp. dec is similar to the primitive function Decode (⊥) except that the left argument is a single multiprecision number, and the right argument is a vector of multiprecision numbers (so it must be enclosed if there is only one). enc is similar to the primitive function Encode (⊤) except that it only works on one multiprecision number at a time, the left argument is a single multiprecision number which is assumed to repeat as often as required.

    ∇ z←b dec x
[1]    z←0 0
[2]   ∆1:z←(⊃x)Fadd b Fmul z
[3]    →(0≠⍴x←1↓x)/∆1
    ∇

    ∇ z←b enc x
[1]    z←⍬
[2]   ∆1:→(b Fgt x)/∆2
[3]    x←x Idiv b
[4]    z←x[1],z
[5]    x←⊃x
[6]    →∆1
[7]   ∆2:z←(⊂x),z
    ∇

Exercises¶

You haven’t got off scot-free: now it’s your turn to do some work. Using the values of the parameters and the encryption and decryption functions above, decipher the following message.

LFWBQ CXLSF EBSPQ JNROT MPAKI VWHOG URVYB DHGFY LOQBH JXCOR
YVFHL ZGWRT XYPRB GCKUT BOPHN ILFHQ MZSOM OHMJO ILEXK XYVSN
QZXLC JGKED ZKWTR SCHXU GNPHV MYGGE VUXCP EAGTB M

Finally¶

I have enjoyed writing this series of articles, and I hope you have enjoyed reading them. I am always interested in hearing about multiprecision arithmetic and its applications, and if you wish to correspond you can get hold of me through the editor.

The workspace that I used for the examples in this series is available by one of the following methods:

If you are on the Internet you can get a vendor-independent version in Jim Weigang’s APLASCII format by anonymous ftp from archive.uwaterloo.ca, called /languages/apl/workspaces/mult.prec/sullivan/MP

Otherwise send me a blank diskette (you can get hold of me via the Editor). You can either have MP as above, or a Dyalog APL/W 7.1.2 workspace which contains my originals and a lot more.

References¶

[1] Fred Gruenberger, How to handle numbers with thousands of digits, and why one might want to. "Computer Recreations", Scientific American, April 1984.

[2] R L Rivest, A Shamir & L Adleman, A method of obtaining digital signatures and public-key cryptosystems. Communications of the ACM, Vol. 21 No. 2, 1978

[3] Hans Riesel, Prime numbers and computer methods for factorization (Birkhäuser, 1985).

Not specifically referred to in the text, I have quoted from and made use of some of the items in Peter Merritt’s articles in Vector, Vol. 11 No. 1 p.108, and Vol. 11 No. 3 p.119.

In addition to Riesel, there are some other books that contain detailed descriptions of the RSA algorithm, amongst which are the following:

C Bondi (editor), New Applications of Mathematics, Penguin Books, 1991

H Beker & F Piper, Cipher Systems, the Protection of Communications, Northwood Books, London, 1982.

D E Knuth, The Art of Computer Programming, Vol. 2, Seminumerical Algorithms, Addison-Wesley, 2nd edition, 1981