It was thus said that the Great Jerome H. Fine once stated: What language? In assembly addition and subtraction are pretty trivial if you have a carry bit you can use (not all RISC based CPUs have a carry bit). Just make sure you keep track of endian issues (given an array A[n] is A[0] the Most Significant or Least Significant word?). So, in Intel x86, the addition routine may look like (A[0] = LSW) multi_byte_add: mov ebx,[esp+8] ; source A value mov esi,[esp+12] ; source B value mov edi,[esp+16] ; destination mov ecx,[esp+20] ; size (in 32-bit words) clc ; clear carry bit .mba_loop: mov eax,[ebx] adc eax,[esi] mov [edi],eax lea ebx,[ebx+4] lea esi,[esi+4] lea edi,[edi+4] dec ecx jnz .mba_loop ret Subtraction works simularly. A full multiplication routine is a bit more daunting---if you know that you are only going to multiply a big number with a word, then it's a bit easier, but still, it comes down do how one does multiplication by hand: A B C D W X Y Z --------------------- AZ BZ CZ DZ AY BY CY DY AX BX CX DX AW BW CW DW (and don't forget the carries!) One trick you might want to try is to take advantage of the fact that these are binary values so multiplication becomes shifts and adds but with larger values like these, doing the multibyte rotates may not be as effective as using the actual MUL instructions. Division is nasty, but again, works like long division. Doing this stuff in C is possible, but you loose out in speed big time. A multibyte addition will look like (warning---untested code! and A[0] = LSW): #include <limits.h> void multi_byte_add(unsigned short *s1,unsigned short *s2,unsigned short *dest,size_t size) { int result; int carry; for (carry = 0 ; size ; size-- , s1++ , s2++ , dest++) { result = *s1 + *s2 + carry; if ((carry = (result > USHRT_MAX))) result -= SHRT_MAX; *dest = (unsigned short)result; } } I have some C code that does some of this---one is a C implementation of a techique Woz used in calculating e to 100,000 digits (as described in Byte magazine of either Dec 1980 or Dec 1981) and another one that does division in C, if you are interested. -spc (Good luck)

Jerome Fine replies: While the logic will probably be in FORTRAN 77 under RT-11 on an emulated PDP-11, I might actually check that it all works on a real PDP-11. However, since FORTRAN 77 is rather nasty using more than 32 bit integers, I will probably code the routines in MACRO-11 for 64/128 bit add and compare after I first write the program for 32 bit integers which can calculate up to 2 billion. As for endian issues, I can define the most significant word to be at either end, but I tend to think that I will place the most significant at the high end of the word since this allows FORTRAN 77 to use octal to output the value correctly with 64 bit floating point values. After that, if the run times are reasonable, I can extend the code to allow 64 bits - although I am not sure why since even running the code to 20 billion takes more than 10 times 2 billion and doesn't really prove anything. Looking at google, I have noted that PI(N) has been successfully calculated up to: N = 40,000,000,000,000,000,000,000 which is already over 70 bits, although it was not done by finding the prime numbers, but by calculating Li(N) and then the error between PI(N) and Li(N). For the PDP-11, I tend to understand that MUL is only for signed 16 bit numbers, so I don't think MUL can be useful with unsigned values which are required for a 64 bit multiple as per the above example. CAN SOMEONE PLEASE CORRECT ME IF I AM WRONG??? That said, in MACRO-11, it should not be too difficult to do multiple shift and add operations. However, maybe someone already had some code handy? I am interested!!!!!!! Thank you for your reply! Sincerely yours, Jerome Fine -- If you attempted to send a reply and the original e-mail address has been discontinued due a high volume of junk e-mail, then the semi-permanent e-mail address can be obtained by replacing the four characters preceding the 'at' with the four digits of the current year.

Jerome Fine replies: I suspect that there is gross confusion over pi (to how many figures of accuracy you choose) and PI(N) which is the NUMBER of primes between 2 and the number N. NOTE that in both cases, the GREEK letter "pi" is used (I do not know how to type GREEK letters under Netscape!!). Why the GREEK letter "pi" is used for PI(N) is beyond my comprehension, but if you google "prime numbers" and look at a few of the references, you will find it there. For example: N PI(N) 10 4 (the primes are 2, 3, 7, 9) 100 25 1,000 168 10,000 1,229 100,000 9,592 By the way, I just learned this use for "pi" this year even though I have been curious about prime numbers for over 20 years. I finally requested some books from the library and then did a few googles. That is too easy a solution - anything a bit more difficult? More seriously, if I am using all 16 bits as unsigned integers, then the high order bit can't be easily eliminated. Of course another method is to keep all the values below 10,000 decimal which can then allow me to use FORTRAN much more easily. But I really prefer to use the full 64 bits. On the other hand, if I use the sieve method to find the primes, I don't need either multiplication or division, just 64 bit adding and 64 bit compares, probably increments by 2 and 4 and of course conversion of 64 bit values to decimal output. The last can also be done by repeated subtraction of powers of 10, so again no division. Thank you for the files. I have not had a chance to understand them as yet, but they will prove to be very useful. Sincerely yours, Jerome Fine -- If you attempted to send a reply and the original e-mail address has been discontinued due a high volume of junk e-mail, then the semi-permanent e-mail address can be obtained by replacing the four characters preceding the 'at' with the four digits of the current year.

... * The comment says 7 bytes long; the code actually does 8 (7..0) * Why not use movzbw instead of clrw-and-movb? * It's canonically bicw, not bitcw. * You need bicw #ff00,r6, not #00ff. * movzbw r6,r6 does what your bi[t]cw tries to, and is shorter. * If you'd store the numbevrs little-endian, the VAX way, instead of big-endian, you could do it much more cheaply with emul - especially if you do use a multiple of 4 bytes. Provided you have integer overflow traps turned off, the difference doesn't matter; the only difference between signed and unsigned multiply is what constitutes arithmetic overflow. Here's what I'd do, which rolls all the above suggestions together: ; All numbers are 8 bytes long, unsigned, LSB-first (native-endian) ; Overflow is ignored (high 64 bits of the 128-bit product are lost) ; r0 - points to first byte of multiplier A ; r1 - points to first byte of multiplier B ; r2 - points to first byte of result ; r3,r4,r5 - ignored on input and garbaged on output emul (r0),(r1),$0,r3 emul (r0),4(r1),r4,r4 emul 4(r0),(r1),r4,r4 movq r3,(r2) If you want the whole 128-bit product, then ; r0 - input A, 8 bytes ; r1 - input B, 8 bytes ; r2 - output, 16 bytes ; r3-r9 - scratch emul (r0),(r1),$0,r3 emul (r0),4(r1),r4,r4 emul 4(r0),4(r1),r5,r5 emul 4(r0),(r1),r4,r8 addl2 r9,r5 adwc $0,r6 movl r3,(r2) movl r8,4(r2) movq r5,8(r2) What VAX doesn't? Except you don't. You probably want EDIV, which does give both. /~\ The ASCII der Mouse \ / Ribbon Campaign X Against HTML mouse at rodents.montreal.qc.ca / \ Email! 7D C8 61 52 5D E7 2D 39 4E F1 31 3E E8 B3 27 4B

Jerome> If the code could be written in FORTRAN, it would be fairly Jerome> simple. Then common source code could be used. But, even if Jerome> possible, FORTRAN is simply not structured to easily allow Jerome> multiple operations on succeeding values which are required Jerome> when adding carry bits between even two words, let alone 8 * Jerome> 16 bit words needed for 128 bit arithmetic. Which means Jerome> coding in assembler! Actually, carry is easy. There's a simple trick I learned on the Alpha, but it applies anywhere. Given UNSIGNED variables: c = a + b; if (c < a) { /* you got a carry */ } Tweak as needed for signed. Jerome> Since I do NOT find that to be a problem, I am not very Jerome> concerned. But in case there was already a library around, I Jerome> thought I would ask. Jerome> By the way, are there any standard algorithms for the 4 basic Jerome> operations (add, subtract, multiply and divide) for 128 bit Jerome> numbers which are composed of 8 * 16 bit words? As per your Jerome> suggestion, I would probably use: CHARACTER * 16 ARRAY ( nnn Jerome> ) Knuth vol. 1 may help. Reading GNU MP may also help, though that is likely to be optimized for much bigger integers (1024 bits or more). paul

Jerome Fine replies: Fortunately on the PDP-11, the carry is automatic in assembler, just "ADC" to the next word to be included in the sum before adding the next pair of words. Likewise for subtraction. Multiplication and addition are a bit more complicated, of course. Sincerely yours, Jerome Fine -- If you attempted to send a reply and the original e-mail address has been discontinued due a high volume of junk e-mail, then the semi-permanent e-mail address can be obtained by replacing the four characters preceding the 'at' with the four digits of the current year.

Jerome> What I am complaining about is not the internal operations, Jerome> but the conversion routines. I went to the trouble of Jerome> printing out the values in OCTAL so that I could figure out Jerome> the problem. A 56 bit mantissa should allow a value of Jerome> 36,028,797,018,963,968 when the leading (invisible) bit is Jerome> one with the remaining bits zero. The displayed value is Jerome> 36,028,797,018,963,965 instead or less by 3. ... That would be a bug. Not surprising. Floating point arithmetic should be, and can be, accurate to the last bit. But very often it is not, because it's *hard* to do that. DEC had a team of specialists who worked on the numeric algorithms for VAX. The PDP-11 libraries didn't get the benefit of that work, so it is not surprising to find errors in the last few bits. As I said, if you want exact answers, DO NOT use float. Jerome> Any idea where I could download the source code? Given that Jerome> FORTRAN is not really suitable for shifting and using carry Jerome> bits, I would port to MACRO-11. It's right in the GNU software library listing: http://directory.fsf.org/GNU/gnump.html Jerome> In general, sieve methods for prime numbers don't ever need Jerome> division and rarely need multiplication. Addition and Jerome> comparing are the primary operations. Increment is also Jerome> needed (with increment by 2 and 4 being a big time saver). Jerome> So the only other thing needed are accurate conversion Jerome> routines to allow output in decimal instead of octal. I Jerome> presume that the latter is best done using subtraction of Jerome> multiple powers of TEN which are stored in a table unless Jerome> there are better methods of converting the output? Any Jerome> suggestions? I'm sure I/O is part of GnuMP. Decimal? Probably, I haven't looked. Multiplication is very definitely in there, it has to be, that's a critical function for crypto. Division I don't know about either. It certainly should be there if GnuMP was aiming to be a general arithmetic package. If all else fails, it certainly has modular arithmetic (again, because of crypto) from which you can synthesize division if it isn't already there. Oh yes, I'm quite sure it also includes primality tests -- presumably the Miller-Rabin test. paul

Jerome> I am not familiar with the Miller-Rabin test? If it has Jerome> anything to do with prime numbers, I would be very Jerome> interested. Depending on how deep you want to dig, if you like primes, you should get yourself some cryptography textbooks. A good primer, without a lot of serious math, is Bruce Schneier's "Applied Cryptography". If you want to understand the math in depth, get "Handbook of Applied Cryptography" (I think that's the correct title) by Menezes et al. I have both, though I'm stuck at around chapter 2 in the latter because the math goes beyond what I know. I think Knuth has a discussion of the M-R test -- surely you have Knuth, right? If not, you need to fix that. The M-R test is a "probabilistic primality test" -- a fast way of testing a number for being prime that has a certain probability of giving the wrong answer. You can do the test repeatedly, and if the tests are done right so they are independent, the test will become more and more accurate. This test is generally used in crypto programs because it is fast. There are also tests for primality that are exact rather than probabilistic -- those are MUCH slower. For multi-thousand-bit numbers, people use M-R unless they absolutely, positively require a guarantee of primality and can afford the hours or days of CPU time investment. An example: the Diffie-Hellman parameters for the IKE key exchange protocol (part of IPsec) were subjected to the full primality test, not just M-R. paul