ljz@fxgrp.UUCP (Lloyd Zusman) (10/13/87)
Has anyone out there come across any algorithms for doing arithmetic on 64-bit
integers? In other words, we want to define a new data type something like
this ...
typedef struct {
unsigned long low_order;
unsigned long high_order;
} ExtPrecision;
Then, we want to be able to do addition, subtraction, multiplication, and
division on these ExtPrecision things. Furthermore, we need to be able to
convert between these things and ASCII strings so they can be input and
output from a terminal.
These routines need to be in C (i.e., NOT in assembler), need to be totally
machine independent (with the possible exception of byte-ordering
considerations which can be handled via conditional compilation), and need
to be as fast as possible. They MUST use integer arithmetic, NOT floating
point. The routines must look something like ...
ExtAdd(result, a1, a2)
ExtPrecision *result;
ExtPrecision *a1;
ExtPrecision *a2;
{}
ExtSubtract(result, s1, s2)
ExtPrecision *result;
ExtPrecision *s1;
ExtPrecision *s2;
{}
ExtMultiply(result, m1, m2)
ExtPrecision *result;
ExtPrecision *m1;
ExtPrecision *m2;
{}
ExtDivide(result, remainder, d1, d2)
ExtPrecision *result;
ExtPrecision *remainder;
ExtPrecision *d1;
ExtPrecision *d2;
{}
I have already written the routines for adding, subtracting, and multiplying,
and they are very fast. I'm using the "long division" algorithm in Knuth
which is accurate but slow as hell. I need something that's exact (insofar
as it produces an exact quotient and remainder) and fast for division, and
alternatively, I could use an entire 64-bit math library.
Any suggestions as to where I might find such a thing, or where I can find
a fast extended-integer division algorithm, or where I can find someone
who could help me with this project?
Thank you in advance for your quick replies.
Lloyd Zusman
...!ames!fxgrp!ljzhunt@spar.SPAR.SLB.COM (Neil Hunt) (10/13/87)
In article <116@fxgrp.UUCP> ljz@fxgrp.UUCP (Lloyd Zusman) writes: >Has anyone out there come across any algorithms for doing arithmetic on 64-bit >integers? In other words, we want to define a new data type something like >this ... > > typedef struct { > unsigned long low_order; > unsigned long high_order; > } ExtPrecision; > >Then, we want to be able to do addition, subtraction, multiplication, and >division on these ExtPrecision things. Furthermore, we need to be able to >convert between these things and ASCII strings so they can be input and >output from a terminal. If you are on a UNIX machine, I hope that you looked at MP(3) MP(3X) MISCELLANEOUS FUNCTIONS MP(3x) NAME mp, itom, madd, msub, mult, mdiv, mout, pow, gcd, rpow, xtom, mtox, mfree - multiple prevision integer arithmetic. It is written in C (at least all the bits I looked at), and based upon the algorithms in Knuth2, with references. Hope this helps. Neil/.
mkhaw@teknowledge-vaxc.ARPA (Mike Khaw) (10/14/87)
in article <1460@spar.SPAR.SLB.COM>, hunt@spar.SPAR.SLB.COM (Neil Hunt) says: +----------- |>Has anyone out there come across any algorithms for doing arithmetic on 64-bit |>integers? ... ... | If you are on a UNIX machine, I hope that you looked at MP(3) | | MP(3X) MISCELLANEOUS FUNCTIONS MP(3x) | | NAME | mp, itom, madd, msub, mult, mdiv, mout, pow, gcd, rpow, | xtom, mtox, mfree - multiple prevision integer arithmetic. +----------- What kind of Unix machine? Sun 3.2 has it, Ultrix 1.2 doesn't, ... Mike Khaw -- internet: mkhaw@teknowledge-vaxc.arpa usenet: {uunet|sun|ucbvax|decwrl|uw-beaver}!mkhaw%teknowledge-vaxc.arpa USnail: Teknowledge Inc, 1850 Embarcadero Rd, POB 10119, Palo Alto, CA 94303
ljz@fxgrp.UUCP (Lloyd Zusman) (10/14/87)
In article <1460@spar.SPAR.SLB.COM> hunt@spar.UUCP (Neil Hunt) writes: >If you are on a UNIX machine, I hope that you looked at MP(3) > > MP(3X) MISCELLANEOUS FUNCTIONS MP(3x) [ followed by the first few lines of the 'man' page which describes these multiple precision routines already available udner Unix ] Thank you VERY much! I feel a little stupid with this stuff sitting right under my nose. However, I ran these routines in a benchmark and found them to be VERY slow compared to the Add, Subtract, and Multiply routines I had already written (on an order of approximately 20-to-1 or more!). I presume that this is because these MP routines are intended for calculations of arbitrary precision, while my routines are especially written for 64-bit calculations. Does anyone know where I can find the source code to these MP routines or to something similar? I'd like to borrow the division algorithm and hack it to work with my 64-bit thingies. Thank you very much. -- Lloyd Zusman ...!ames!fxgrp!ljz
Alan_T._Cote.OsbuSouth@Xerox.COM (10/17/87)
The following C function should get you started. I once coded the
algorithm in 68000 assembler, and benchmarked on a 68010 as faster than
the corresponding instruction on a 370 mainframe.
Enjoy.
--------------- CUT HERE ------------
typedef struct s_bigint
{
long low;
long high;
} BIGINT;
/* The bigdiv() function performs signed division of a 64-bit signed
* dividend by a 32-bit signed divisor, yielding a 32-bit signed
* quotient and a 32-bit signed remainder. The approach taken is that
* of binary long division. The method is easily extended to cases
* in which all operands are 64-bit, to deal with the situation where
* the quotient must be >32 bits (which can happen -- just try dividing
* 0x100000000000000 by 2, yielding a 32-bit result!).
*
* If the division is successful, bigdiv() returns 0. If quotient
* overflow is detected (eg: division by 0), 1 is returned.
*/
int bigdiv( dividend,
divisor,
quotient,
remainder )
BIGINT *dividend;
long divisor;
long *quotient;
long *remainder;
{
unsigned long dividend_high;
unsigned long dividend_med;
unsigned long dividend_low;
unsigned long temp_quotient;
unsigned long temp_divisor;
int quotient_is_negative;
int bit;
/* Initialize */
temp_quotient = 0;
dividend_high = 0;
quotient_is_negative = 0;
if( dividend->high < 0 )
{
quotient_is_negative = !quotient_is_negative;
dividend_med = ~ dividend->high;
dividend_low = ( ~ dividend->low ) + 1;
if( dividend->low == 0 )
if( ++dividend_med == 0 )
++dividend_high;
}
else
{
dividend_med = dividend->high;
dividend_low = dividend->low;
}
if( divisor < 0 )
{
quotient_is_negative = !quotient_is_negative;
temp_divisor = -divisor;
}
else
temp_divisor = divisor;
/* Do the division */
for( bit = 0; bit < 64; bit++ )
{
/* check for overflow */
if( ( temp_quotient & 0x80000000L ) != 0 )
return( 1 );
/* shift dividend and temporary quotient left one bit */
dividend_high = ( dividend_high << 1 ) + ( ( dividend_med >> 31 ) & 1
);
dividend_med = ( dividend_med << 1 ) + ( ( dividend_low >> 31 ) & 1 );
dividend_low <<= 1;
temp_quotient <<= 1;
/* if dividend_high is divisible by temp_divisor, divide them */
if( dividend_high >= temp_divisor )
{
dividend_high -= temp_divisor;
++temp_quotient;
}
}
/* Establish final results */
if( quotient_is_negative )
*quotient = -temp_quotient;
else
*quotient = temp_quotient;
if( dividend->high < 0 )
*remainder = -dividend_high;
else
*remainder = dividend_high;
return( 0 );
}