A Simple
Implementation of Extended Precision Integer Arithmetic using the TI-85 Graphing
CalculatorI have a passing interest in basic number theory, and one facet that particularly interests me is long integer arithmetic. I'm not talking about signed 32-bit integers, which is what the term long integer usually means in computer languages like Pascal, but rather numbers which can be much bigger (the programs here can produce output hundreds of digits long), often called "large integers".
What I plan on addressing in this document is a simple system for basic positive long integer math operations. Here we'll cover addition, exponentiation, multiplication, and factorials: the operations that are closed under the set of positive integers.
Note that this is in no way meant to be a guide to implementing large integer math; the method proposed here makes a number of choices that could, in a full program on a full computer, be greatly and easily improved upon. For an excellent overview, see Knuth's Art of Computer Programming: Seminumerical Algorithms.
For starters, we'll need a way to store long integers. For the sake of simplicity, we'll consider a positive integer as a list of digits. For example, a number like 13 is stored will, for our purposes, be stored as the two element list [ 1 3 ].
Now that we have our representation, we're almost ready to start theorizing about possible operations and the algorithms we can use to perform them. Before we really get going, we'll need to understand how our approach works.
One nagging question remains. Why are we allowed to represent our number as a list of digits? This may seem like a silly question, but it isn't.
First off, you say to yourself: "Self... How do I get from my list of digits to my actual number?" It seems like you just take each number from the list and write it down in the order it appears. But what are you doing when in this process?
What we're doing is looking through the list and multiplying each element by a corresponding power of 10. At the last element of the list, which we know as the units digit, we're multiplying by the 0th power of 10, or 1. At the second to last element we multiply by 10¹. At the third from last element we multiply by 10²... and so on and so forth. What we've done is basically used 10 as are radix. (This isn't the method: it is better to use a power of 2 for the radix, but for simplicity and my sanity's sake we'll stick with something easy like 10).
[ an an-1 an-2 an-3 ... a1 ]
From the prior discussion, we know that the above list depicts a number= 10n-1*an
+10n-2*an-2 + ... + 100*a1
=
Does the above sum look familiar? If we let x=10, it's easy to see that our sum is simply a polynomial in x of degree (b-1). The coefficients are given by ab. This means that we can use the rules of algebra to justify and perhaps to find a means to operate on long integers.
The first two operations (addition and multiplication) will operate on long integers, as opposed to just generating long integer results. Because of this we have to worry about input routines. How can make the TI-85 accept long integers as input?
Certainly we can't use its normal storage functions. Entering a long integer on the command line would simply result in a truncated real value being stored to the variable. We could have the user enter a list of digits, which is what we'll need to work with internally, but that would be far too cumbersome. The only other viable method is through the string.
How do we plan on changing string input to a list? The following source code should do the trick:
\START85\ \COMMENT=Program file dated 11/07/96, 21:23 \NAME=ST2L \FILE=C:\LINK\ST2L.85P :lngth ST\->\dimL NL :For(I,1,lngth ST,1) :sub(ST,I,1)\->\P :\St>Eq\(P,G) :G\->\NL(I) :End \STOP85\
What exactly does the above algorithm do? It's essentially a brute force algorithm that runs a for-next loop, finds the current character in the string (ST) and then does something peculiar -- it treats the character it happens to be on as a string to be converted into an equation. The string "3" becomes the equation 3, which is treated almost exactly as the real number 3 (in fact, it becomes a real number after the first operation performed on it)! We then transfer the Ith number in the string to the Ith number in the list, which leaves us with a list of digits!
After we've done our internal work, we'll need some way to output the results in a pretty, usable form. To do this we'll need an algorithm to work backwards. How will we do this?
Previously I had decided to simply go through the list, number by number, each time adding the current number to the string. What I ended up with was a huge progra, to perform a task that should have been relatively simple. The code was neither efficient nor as fast as I would have liked -- luckily, Mark Lybrand stepped in with an excellent replacement algorithm.
His algorithm, which is quite efficient, cycles through the list's elements. When it encounters digit I, it goes to the (I+1)st character of the string, which contains the string equivalent of the digit. For example, if it encounters a 0 it goes to the (0+1), or 1st character of the string, which is "0" -- it's easy to see how and why this new algorithm works out so nicely.
Also note the beginning of the program, which advances the cursor through the list to start at the first non-zero number.
\START85\ \COMMENT=Program file dated 11/19/96, 21:39 \NAME=L2SD \FILE=C:\85\L2SD.85P :"0123456789"\->\DIGITS :" "\->\ST2 :1\->\B :While C(B)==0 :1+B\->\B :End :For(I,B,dimL C,1) :sub(DIGITS,C(I)+1,1)\->\TEMP :ST2+TEMP\->\ST2 :End :sub(ST2,2,lngth ST2-1) \STOP85\
Ok, we know how to get our numbers into the calculator. How do we perform our operations on the resulting lists?
We'll try addition first. If we have [1 3 1 2 4 1] and another [2 7 3 2 1 2 1], we obviously only have to add the lists and, if a number is >10, carry the overflow to the next number. Using the given numbers (lists), we wind up with [3 10 4 4 5 3], which changes into [4 0 4 4 5 3], or, when processed, 404453. Why is this justified?
This is justified because our lists simply represent polynomials in (10^x). If one of its coefficients is >10 (say 10a+b), we get a*10x+1 + b*10x -- we effectively carry the overflow!
Now that we know how, all that's necessary is to write a simple routine. This is easier said than done, but it is possible, once we can navigate through some obstacles.
One thing that's troublesome is that unless the lists have the same dimension, we can't use the addition operator to add the corresponding elements. Having lists with the same number of elements doesn't help as much as we would like anyway: if we have lists like [6 2 1 3 2 1] and [5 3 2 1 2 3], we're going to need an extra digit in our final result. We can't just add them!
To get around this problem, we'll invent a list for our answer that has one more digits than the list with the most elements. For example, if we're adding 321983298123 and 658937211221132, we'd dimension a new list with 16 elements (because 16 is one more than the number of digits in our number with the greatest number of digits).
We'll first work out the special case, when both numbers have the same number of digits -- say, n. Then we dimension a list with 2*n+1 elements.
We then transfer the numbers from the list to their corresponding positions in the new list. For example, if we had the lists [1 3 2 4 2] and [3 2 1 4 5 7], we'd first dimension a list of seven elements. After this we'd add put the second number in the last six positions to have [ 0 3 2 1 4 5 7] and then add the first number to the final positions, to get [0 3 3 4 2 9 9]. All we have to do is check for carries and perform them, if necessary and we have our result. This is the principal used in the following source code:
\START85\ \COMMENT=Program file dated 11/07/96, 21:23 \NAME=ADD \FILE=C:\LINK\ADD.85P : ClLCD :0\->\C :Disp "ADDITION" :Disp "" :InpST ST :ST2L :NL\->\AB :InpST ST :ST2L :NL\->\AC :(dimL AB)*(dimL AB>dimL AC or dimL AB==dimL AC)+(dimL AC)*(dimL AC>dim\#\ L AB)+1\->\dimL C :If dimL AB>dimL AC :Then :AC\->\LL :AB\->\AC :LL\->\AB :""\->\LL :End :If dimL AB==dimL AC :Then :For(P,2,dimL C) :AB(P-1)+AC(P-1)\->\C(P) :End :Else :For(P,2,dimL C) :AC(P-1)\->\C(P) :End :For(P,dimL C-dimL AB+1,dimL C,1) :C(P)+AB(P-dimL C+dimL AB)\->\C(P) :End :End :FIXL :L2SD \STOP85\
Suppose we want to find 3490529510847650949147849619903898133417764638493387843990820577 × 32769132993266709549961988190834461413177642967992942539798288533? Obviously no human (with the possible exception of me and a few of my friends) would be insane enough to perform this computation by hand! Surely there must be a better way!
There is. Programs to perform multiplication are relatively easy to write and can be quite speedy. When multiplied out, we find an amazing number, RSA-129:
1143816257578888676692357799761466120102182967212423625625618429357069 35245733897830597123563958705058989075147599290026879543541
Going backwards from that number to its factorization is quite a bit more difficult (and would have taken 40 quadrillion years using only 1977 methods); we'll leave that to other documents (Mark Janeba has an excellent article on the factorization of RSA-129, for those interested). Instead we'll (what a novel idea) get back on topic...
The fundamental concept we'll use depends upon the fact that since a long integer is nothing more than a polynomial of sorts, we can use an algorithm to multiply the polynomial. Once we have the resulting polynomial and we perform the necessary carries, we'll have our answer.
To multiply a polynomial we first consider a simple case and the operations we would normally perform if working out such a problem by hand. We'll start with the lists [1 2 3] and [4 2] (x²+2x+3 and 4x+2, respectively). To multiply these we note that we go through the first polynomial, multiplying each term by each of the terms of the second and making sure the degree of our x terms are correct, so that we can later combine like terms.
However, in a list there is no telltale sign that tells you what degree x currently is. We can multiply each element of the first list by the elements of the second, but how do we know where each term goes in the final result? The answer is that via clever use of dimL and lngth we can determine what degree each term in the list has and combine them appropriately.
\START85\ \COMMENT=Program file dated 11/07/96, 21:23 \NAME=MULT \FILE=C:\LINK\MULT.85P :0\->\C :ClLCD :Disp "Long Integer" :Disp "Multiplication" :Disp "" :InpST ST :ST2L :NL\->\A :InpST ST :ST2L :NL\->\B :dimL A+dimL B\->\dimL C :For(L,1,dimL A) :For(M,1,dimL B) :C(dimL C-(dimL A-L+dimL B-M))+A(L)*B(M)\->\C(dimL C-(dimL A-L+dimL B-M\#\ )) :End :End :FIXL :L2SD \STOP85\
Work your way through a few simple examples; why the above algorithm works should become clear to you.
Since carries are so often needed in long integer arithmetic, we'll develop a special routine to do exactly this. Our algorithm will simply go from the end of the list to the front, carrying the appropriate excess to the next list position. It looks like this:
\START85\ \COMMENT=Program file dated 11/07/96, 21:23 \NAME=FIXL \FILE=C:\LINK\FIXL.85P : For(I,dimL C,2,\(-)\1) :(C(I)-mod(C(I),10))/10+C(I-1)\->\C(I-1) :mod(C(I),10)\->\C(I) :End \STOP85\
The carries are performed by finding the units digit (via mod), and adding what's left (once divided by 10) to the next position!
First off, we'll define some ground rules. We'll only except integers (not long integers) as input. For factorials, this is obviously not a problem -- for exponents, it's somewhat limiting, but much easier to implement this way.
The algorithm used in PWR and FCTRL is similar: they both rely on the fact that to multiply a number stored this way, all I have to do is multiply each of its elements by the desired quantity and then use the above list cleaning up algorithm on it.
For example, to multiply 3921 by 22, I do something like this: [ 3 9 2 1 ] * 22= [ 66 198 44 22 ]. I then run the previous routine to carry the number, which gives [8 6 2 6 2]. I then just have to run L2SD and I have my number!
PWR just does this a number of times: to calculate 8³, all I do is find out how many digits 8³ has by truncating the log of the answer and adding one. To do this, I note that the logarithm of 8³=3*log(8). This gives me three digits.
[0 0 8] ( I create a list with 3 (or more) places since 8^3 has 3 digits. I then place the number in the last spot and multiply by 8 the remaining (power-1) times -- in this case, 2. ) [0 0 64] ( ×8 and then LONGC here to clean it up ) [0 6 4] ( cleaned up output ) [0 48 32] ( multiplied by 8 again ) [5 1 2] ( cleaned up )
Here's the relevant source:
\START85\ \COMMENT=Program file dated 11/07/96, 21:23 \NAME=PWR \FILE=C:\LINK\PWR.85P :ClLCD :0\->\C :Disp "EXPONENTS" :Disp "" :Input A :Input B :iPart (B*log A)+1\->\dimL C :A\->\C(dimL C) :For(M,1,B-1,1) :If max(C)*A>2\E\13 :Then :FIXL :End :C*A\->\C :End :FIXL :L2SD \STOP85\
FCTRL uses a similar process which you should be able to guess by now... you're also probably wondering why I clean-up the output after each time. I was wondering that as well, and that's why I don't. In the programs below, I only clean up the list when the numbers in them good have significant digits lost when performing the multiplication.
Say I was working with a stupid machine that could only handle 4 digits in a list and I want to take 85. I have
[ 0 0 0 0 8 ] --> [ 0 0 0 0 64 ] [ 0 0 0 0 512 ] [ 0 0 0 0 4096 ] -- Ahhhh! I can't go on any further because if I multiply by 8 again, I'm just going to lose significant digits! I just clean up the list right now to prevent it... [ 0 4 0 9 6 ] [ 0 32 0 72 48 ] (after cleaning up): [ 3 2 7 6 8 ]
This does bring some limitations into it; since the TI-85 only integers up to about 2e13 without loss of significant digits, the base of the number I'm multiplying by can't be any bigger than its (2e13's) square root without risking some danger -- hence the 4.2e6 limitation. Hence, to make 2^300 execute faster, you could use 220 as your number and 15 as your power, but you couldn't use 230 as your number with 10 as your power. Again, the relevant source follows.
\START85\ \COMMENT=Program file dated 11/07/96, 21:23 \NAME=FACT \FILE=C:\LINK\FACT.85P :0\->\C :ClLCD :Disp "N FACTORIAL" :Disp "" :Input N :iPart (log (N!))+1\->\dimL C :N\->\C(dimL C) :For(K,N-1,2,\(-)\1) :C*K\->\C :If max(C)*K>2\E\13 :Then :FIXL :End :End :FIXL :L2SD \STOP85\
There are a few peculiarities of these programs. The concept of losing significant digits (and preventing this via FIXL) may still be vague to some. What I do in these programs is check to see if another step will result in a loss-- if this is true, I clean up (or perform the carries) on the list.
That about wraps it up. Drop me a note with any questions or comments.
There are two implementations currently available. The preferred version (in .85G format), using precisely the algorithms outlined above, is available for download in ZIP format (~1k) from this site.
Older versions of the TI-85 programs outlined here are publicly available. If you have a Graph-Link you can download the program group (~2k) and transfer it to your calculator.
Note that for those lacking the Graph-Link, the source code for the current version of these programs is also available.
![[Wilbur validated]](checked2.gif){kind=small}
Written, maintained and validated by Paul
Pollack