Dividing 512 bit number by 128 bit number in C program
vaidehi...@gmail.com
I have a huge problem finding " a mod m " using C program, when m is
128 bit number and a is 512 bit number. Can anybody suggest me a
solution.
Thanks, in advance.
Le Chaud Lapin
Well, I can tell you what not to do: that which I did.
I did what I did first because I was hopeful to just get something
working quickly, and then out of curiosity to see just how bad my
performance would be. The performance was really bad. I will direct
you to what you should do instead after talking a bit about what I did,
which is probably what many do who find themselves on the path of big
integer division.
WHAT YOU SHOULD NOT DO:
There might be a curious urge to use the elementary school algorithm of
repeatedly aligning the 128-bit value under the far-left end of the
512-bit value, and subtracting, bit by bit, until you can subtract no
more.
For each iteration, you would find the highest-order bit of the 512-bit
value, and place the highest-order bit of the 128 bit value just under
it, then begin the subtraction starting at the lowest-order bit of the
modulus, borrowing as necessary. Before you start your subtraction,
you would check to make sure that simply lining up the dividend
(ultimate remainder) with the modulus does not lead you into a trap -
it might turn out that you need to shift the remainder one extra
position to the right so that that from which you are about to subtract
is actually >= to the modulus. In any case, you are finished when you
cannot perform any more subtractions because the dividend (remainder)
is less than the divisor (modulus). Unfortunately, this algorithm is
very low and cannot be used for any practical purposes other than
teaching.
WHAT YOU SHOULD DO:
Use Knuth's Algorithm D from 4.3.1. I have not any of his excellent
works. I am only going on what I have read. If you search for "Knuth
4.3.1 D hacker's delight" on Google, you can find a section on
multi-word division.
WHAT YOU MIGHT DO:
If just want something that works, I can send you my C++ function that
does division. It works, and could probably be understood and changed
to C easily, but is too slow to be of any practical use.
-Le Chaud Lapin-
vaidehi...@gmail.com
to look at it.
Thanks..
Vaidehi
Le Chaud Lapin
> Thanks a lot this helps. Could you send me your C++ code. I would like
> to look at it.
I decided to post here for benefit of others who would make the same
mistake I made.
Description of the algorithm:
Division is naturally the most complicated of the basic big-integer
operations. The following is a C++ static member function of class
Integer that does big integer division. It takes dividend and divisor
by const reference, and quotient and remainder by reference,
calculating the quotient and the remainder before returning. There is
no return value.
Division by 0 results in an Integer::DIVISION_BY_ZERO exception.
Division by 1 results in remainder being set to 0 and quotient being
set to the dividend.
The function consists mostly of a large outer loop.
In one of the inner loops, the divisor is positioned so that its most
significant bit is just beneath the most significant bit of the
dividend. A bit-wise subtraction of divisor from dividend is performed
from right to left, starting with the low-order bit of the divisor,
borrowing bits if necessary. The positioning of the divisor is
important - it is not enough to simply align the high-order bits and
forget, as it might happen that the bit sequence of dividend involved
in the subtraction yields a number whose magnitude is insufficient to
accommodate a subtraction by the divisor. A small snippet of code
checks for this situation before allowing the subtraction to occur by
scanning the bits of both the dividend and the divisor from left to
right, and if it happens that a dividend bit is found to be 0 while the
corresponding divisor bit is 1, then naturally, subtraction is
impossible. In this case, a simple shift of divisor one bit to the
right will remedy, but if the divisor cannot be shifted because it is
already as far to the right as it can go (the divisor_offset is 0)
then naturally we are done with algorithm.
As we perform the subtraction, we pay careful attention to borrowing,
taking all tri-state cases into consideration. Those familiar with the
subtraction logic of an ALU will find these cases familiar. I decided
to use IF statements to handle the case instead a truth table and
direct computation for clarity, since the algorithm is already to slow
for practical use.
Each time a subtraction occurs, we set a corresponding bit in the
quotient. This bit just so happens to have the same offset as
divsor_offset.
After each subtraction, if we find that there are fewer bits in the
dividend than divisor, we are done.
The function ::adjust_extent of class Integer is simply present to trim
the memory use of the Integer object.
Bit testing for 1 or 0 is done using array index into the buffer of the
Integer object and C/C++ logical operators.
The sign of the quotient is based on the sign of the dividend and the
divisor. The sign of the remainder, for lack of a better scheme, uses
the exact same logic to compute the sign of the quotient.
Again, this function should not be used in commercial code because it
is too slow. The only use would be to instruct students how *not* to
design a division function in C/C++. An efficient division algorithm
can be found in "D. Knuth, The Art of Computer Programming Vo.2,
Section 4.3.1"
-Le Chaud Lapin-
void Integer::divide (const Integer ÷nd, const Integer &divisor,
Integer "ient, Integer &remainder)
{
// Division by zero is not permitted:
if (divisor.is_zero())
throw Integer::DIVISION_BY_ZERO;
if (divisor.is_one())
{
quotient = dividend;
remainder = 0;
return;
}
quotient.nullify();
quotient.adjust_extent(dividend.extent());
// We copy the dividend into the remainder since the remainder will be
computed by successively reducing it, beginning with remainder ==
dividend:
remainder = dividend;
// Determine the initial number of bits in the remainder:
unsigned int remainder_length_bits = remainder.length_bits();
// Determine the initial number of bits in the divisor:
const unsigned int divisor_length_bits = divisor.length_bits();
// Keep trying to reduce remainder by dividing by divisor until there
are no more spots to shift divisor into (divisor_offset == 0):
while (remainder_length_bits >= divisor_length_bits)
{
unsigned int divisor_offset = remainder_length_bits -
divisor_length_bits;
// Check if the substraction of divisor can be performed:
bool subtractable = true;
unsigned int divisor_bit_index = divisor_length_bits - 1;
do // DEFECT: The high-order bits of both remainder and divisor are
always set so we can skip them.
{
register unsigned int remainder_bit_index = divisor_bit_index +
divisor_offset;
if (remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] & (1 <<
(remainder_bit_index % WORD_LENGTH_BITS)))
{
// Remainder bit is 1.
if (!(divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS))))
// Divisor bit is 0.
break; // ...subtraction can be performed.
// Else, both remainder bit and divisor bit are equal to 1.
}
else
{
// Remainder bit is 0.
if (divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS)))
{
// Divisor bit is 1.
// If we cannot shift the divisor to the right one position...
if (divisor_offset == 0)
subtractable = false; // ...subtraction cannot be performed.
// Otherwise...
else
--divisor_offset; // ...subtraction can be performed after
reduction of divisor_offset by 1:
break; // In any case, we break because we have remainder_bit ==
0, divisor_bit == 1, which is always an interesting situation.
}
// Else, both remainder bit and divisor bit are equal to 0.
}
}
while (divisor_bit_index--);
// If subtraction cannot be performed, we are finished:
if (!subtractable)
break;
// Perform a subtraction:
bool borrow = false;
for (divisor_bit_index = 0; divisor_bit_index < divisor_length_bits;
++divisor_bit_index)
{
register unsigned int remainder_bit_index = divisor_bit_index +
divisor_offset;
if (remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] & (1 <<
(remainder_bit_index % WORD_LENGTH_BITS)))
{
// Remainder bit is 1.
if (divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS)))
{
// Divisor bit is 1.
if (!borrow)
// Set remainder bit to 0 subtract (1-1) and satisfy borrow of 0:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] &= ~(1
<< (remainder_bit_index % WORD_LENGTH_BITS));
// Else leave remainder bit == 1 because borrow made it 0, but to
subtract divisor it need its own borrow, so (2-1) == 1
}
else
{
// Divisor bit is 0.
// In this case, we only disturb remainder bit if borrow needs to
be satisfied.
if (borrow)
{
// Set remainder bit to 0 to satisfy borrow:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] &= ~(1
<< (remainder_bit_index % WORD_LENGTH_BITS));
borrow = false;
}
}
}
else
{
// Remainder bit is 0.
if (divisor.buffer[divisor_bit_index / WORD_LENGTH_BITS] & (1 <<
(divisor_bit_index % WORD_LENGTH_BITS)))
{
// Divisor bit is 1.
if (!borrow)
{
// Set remainder bit to 1 to subtract (2-1). Indicate necessity
to borrow:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] |= (1 <<
(remainder_bit_index % WORD_LENGTH_BITS));
borrow = true;
}
// Else remainder bit should stay 0 because continuing borrow gave
us 2, borrow took 1 and divisor took 1.
}
else
{
// Divisor bit is 0.
if (borrow)
// Set remainder bit to 1 because continuing borrow gave us 2,
borrow took 1:
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] |= (1 <<
(remainder_bit_index % WORD_LENGTH_BITS));
// Else there is nothing to do because all 3 states are 0.
}
}
}
if (borrow)
{
unsigned int remainder_bit_index = divisor_offset +
divisor_length_bits;
remainder.buffer[remainder_bit_index / WORD_LENGTH_BITS] &= ~(1 <<
(remainder_bit_index % WORD_LENGTH_BITS));
}
// Set appropriate quotient bit since subtraction was performed:
quotient.buffer[divisor_offset / WORD_LENGTH_BITS] |= (1 <<
(divisor_offset % WORD_LENGTH_BITS));
while (remainder_length_bits >= divisor_length_bits)
{
if (remainder.buffer[(remainder_length_bits - 1)/ WORD_LENGTH_BITS]
& (Integer::Word(1) << ((remainder_length_bits - 1)%
WORD_LENGTH_BITS)))
break;
--remainder_length_bits;
}
}
quotient.adjust_extent();
quotient.sign_ = static_cast<Integer::Sign>(dividend.sign_ *
divisor.sign_);
remainder.adjust_extent();
remainder.sign_ = static_cast<Integer::Sign>(dividend.sign_ *
divisor.sign_);
}
Christian Siebert
> 128 bit number and a is 512 bit number. Can anybody suggest me a
> solution.
example:
/* gcc -O2 -lgmp a_mod_m.c -o a_mod_m */
#include <stdio.h>
#include <gmp.h>
int main(int argc, char *argv[])
{
mpz_t a, m, r;
int base = 10;
if (argc == 3) {
mpz_init_set_str(a, argv[1], base);
mpz_init_set_str(m, argv[2], base);
mpz_init(r);
mpz_mod(r, a, m);
gmp_printf("%Zd\n", r);
}
else puts("usage: a_mod_m <integer 'a'> <integer 'm'>");
return 0;
}
Phil Carmody
Do not post the same post in multiple newsgroups, cross-post instead.
That way people won't needlessly respond with suggestions that you've
already seen in other newsgroups.
Phil
--
"Home taping is killing big business profits. We left this side blank
so you can help." -- Dead Kennedys, written upon the B-side of tapes of
/In God We Trust, Inc./.
Le Chaud Lapin
> vaidehi...@gmail.com wrote:
> > Thanks a lot this helps. Could you send me your C++ code. I would like
> > to look at it.
>
> I decided to post here for benefit of others who would make the same
> mistake I made.
I would like to make a retraction of the statement in the line above.
I realized last night just before going to bed that my algorithm is not
as useless as I thought. My gut feeling is that it can be made to run
much faster than it is in the implementation I gave.
Just how fast remains a matter of exploration. I will explore.
-Le Chaud Lapin-
Unruh
>Hi all,
>I have a huge problem finding " a mod m " using C program, when m is
>128 bit number and a is 512 bit number. Can anybody suggest me a
>solution.
Get a bignum package. There are a number around-- in pgp, in gpg, in tom's
libcrypt (libtomcrypt?)
>Thanks, in advance.