Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

944 lines (870 sloc) 36.989 kb
/* biggishint.c */
/* Biggish integers in this library are arrays of up 32766 16-bit */
/* (unsigned short) integers for arbitrary precision integer */
/* arithmetic up to 524240-bit (16*32765 bit) numbers. A biggish */
/* rational library (biggishrat) using these is also being developed. */
/* The data format for these (fairly) big integers is an array of */
/* short ints allocated on the heap, with the following layout: */
/* +--------------------------------+-------------------------------+ */
/* | required (2 short ints) | optional (up to 16382 times) | */
/* | word 0 | word 1 | word 2 | word 3 etc | */
/* | 14 bits: size | 16 bits: int | 16 bits: int | 16 bits: int | */
/* | 1 bit: overflow | | | | */
/* | 1 bit: sign | | | | */
/* +-----------------+--------------+---------------+---------------+ */
/* Word 0 uses 14 bits for the size in units of pairs of short ints, */
/* 1 bit for overflow/underflow and 1 bit for (minus) sign. */
/* The remaining words are all unsigned short integers, always an odd */
/* number of them because of word 0 and the even array size. Thus */
/* memory allocation occurs in multiples of 4 bytes, a nice alignment */
/* compromise between small and large processor types. */
/* The most compact biggishint uses 4 bytes to count from -65536 to */
/* +65535, the next size (8 bytes) spans -281474976710656 (-2^48) to */
/* +281474976710655 (+2^48-1), and so on to the largest (65532 bytes) */
/* calculating from about -9.22e157811 (-2^524240) to about */
/* +9.22e157811 (+2^524240-1). */
/* The word order is big endian, no matter what the byte order of the */
/* processor may be. There is no reason for this, it just is. */
/* Update: there is a reason to be little endian: the trim() and */
/* shortmultiply() would improve. */
/* All words are unsigned. Negative numbers are stored as positive */
/* numbers and the first word keeps a sign bit. Thus the functions */
/* behave as if each word is simply a digit in a base-65536 number. */
/* The design would use 32-bit words (base-4294967296 numbers) */
/* instead if every possible C compiler had a 64-bit data type to do */
/* carries, but this is not the case. */
/* Functions set the overflow bit if the result does not fit in 32766 */
/* words (including the initial size word). If the overflow bit is */
/* set the other fields change their meaning: */
/* * The sign bit indicates positive overflow or negative overflow. */
/* * The integer value is unusable, and is therefore reduced to the */
/* minimum of 1 word. A design in which the value might represent */
/* positive or negative infinity (as opposed to mere overflow), */
/* or other forms of NaN (Not a Number) is under consideration. */
/* Use biggishint at your risk and without warranty. Give due credit */
/* if you do. Written by Martin Berends. */
/* See also: a much bigger library: http://gmplib.org/manual/ */
/* Donald E. Knuth The Art of Computer Programming Vol 2 */
/* TODO: overflow detection */
/* TODO: change from big endian to little endian */
#include <assert.h> /* assert */
#include <ctype.h> /* tolower */
#include <limits.h> /* USHRT_MAX */
#include <stdio.h> /* printf, only when debugging */
#include <stdlib.h> /* calloc malloc realloc free */
#include <string.h> /* strlen strncmp */
#include "biggishint.h" /* (all externally callable functions) */
#if USHRT_MAX != 65535
#error In this C compiler a short int is not a 16 bit number.
#endif
/* #define BIGGISHINT_TRACE */
/* Internal functions are declared here, their definitions are lower */
/* down. */
unsigned short * biggishint_internal_addsubtract(unsigned short * bi1, unsigned short * bi2, int flipsign2);
int biggishint_internal_bitsize(unsigned long);
unsigned short * biggishint_internal_clone(unsigned short * bi1);
int biggishint_internal_comparemagnitude(unsigned short * bi1, unsigned short * bi2);
unsigned short * biggishint_internal_shiftleft(unsigned short * bi1, unsigned int bitcount);
unsigned short * biggishint_internal_shiftright(unsigned short * bi1, unsigned int bitcount);
void biggishint_internal_shortdivide(unsigned short * bi1, unsigned short * i2);
void biggishint_internal_shortmultiply(unsigned short ** bi1, unsigned short i2);
unsigned short * biggishint_internal_trim(unsigned short ** bi1);
/* --------------------------- Functions ---------------------------- */
/* biggishintAdd */
unsigned short *
biggishintAdd(unsigned short * bi1, unsigned short * bi2)
{
return biggishint_internal_addsubtract(bi1, bi2, 0);
}
/* biggishintCompare */
int
biggishintCompare(unsigned short * bi1, unsigned short * bi2)
{
int sign1, sign2, result;
sign1 = * bi1 & 1;
sign2 = * bi2 & 1;
result = sign1
? ( sign2 ? biggishint_internal_comparemagnitude(bi2, bi1) : -1 )
: (!sign2 ? biggishint_internal_comparemagnitude(bi1, bi2) : 1 );
return result;
}
/* biggishintDivide */
/* see The Art of Computer Programming Vol 2 3rd Ed p270-275 */
unsigned short *
biggishintDivide(unsigned short * bi1, unsigned short * divisor)
{
/* Before starting on the main long division, which is slow, try */
/* to identify divisors that offer an opportunity for a shortcut, */
/* for example shifting right for divisors that are powers of two */
/* or short division for divisors that fit into a small int. */
/* In contast with most of the other routines, this one uses */
/* multiple returns to avoid having many levels of nested */
/* conditionals. */
unsigned short bi1size, divisorsize, dividendword1, dividendword2;
unsigned short * dividend, * pdividend, * pdividendhi, * pdividendlo;
unsigned short * quotient, * pquotient, * pquotienthi, * pquotientlo;
unsigned short * pdivisor, * pdivisorhi, * pdivisorlo;
int sign1, sign2, sign, comparison, divisorshift;
unsigned long dividendcarry, trialdivisor;
unsigned long tempquotient, trialquotientmin, quotientcarry;
bi1size = (* bi1 & 0xfffc) >> 1;
divisorsize = (* divisor & 0xfffc) >> 1;
sign1 = * bi1 & 1;
sign2 = * divisor & 1;
sign = sign1 ^ sign2;
/* Does dividend have fewer words than divisor? */
if (bi1size < divisorsize) { /* quotient becomes 0 */
quotient = (unsigned short *) calloc(2, sizeof(short));
* quotient = 4;
return quotient;
}
/* Is divisor only a 16 bit unsigned number? */
if (divisorsize == 2) { /* use short division instead of long */
unsigned short shortdivisor = divisor[1];
quotient = biggishint_internal_clone(bi1);
biggishint_internal_shortdivide(quotient, & shortdivisor);
* quotient ^= sign2;
return biggishint_internal_trim(&quotient);
}
/* Is dividend less in magnitude than divisor? */
comparison = biggishint_internal_comparemagnitude(bi1, divisor);
if (comparison<0) { /* quotient becomes 0 */
/* Hope the C compiler merges this code with the same code in */
/* (bi1size < bi2size) above. Repeated here because the */
/* earlier case avoids the comparemagnitude function. */
quotient = (unsigned short *) calloc(2, sizeof(short));
* quotient = 4;
return quotient;
}
/* Is dividend equal in magnitude to divisor? */
if (comparison==0) { /* quotient becomes 1 (+ or -) */
quotient = (unsigned short *) calloc(2, sizeof(short));
* quotient = 4 | sign;
quotient[1] = 1;
return quotient;
}
/* Is divisor a power of two or a multiple of a power of two? */
/* ie is there only a single 1 bit or at least one trailing 0 bit? */
if (0) {
; /* TODO: right shift optimization */
}
/* Perform long division because none of the quicker algorithms */
/* above could be used with the given parameters. */
assert( bi1size >= 4 ); assert( divisorsize >= 4 );
/* Initialize dividend with bi1. The loop below will repeatedly */
/* subtract multiples of divisorshifted until the remainder is */
/* less than divisor. */
dividend = biggishint_internal_clone(bi1);
pdividendhi = dividend + (dividend[1] ? 1 : 2);
pdividendlo = dividend + bi1size - 1;
#ifdef BIGGISHINT_TRACE
fprintf(stdout,"dividend %s%04x", (* dividend & 1)?"-":"",
* (pdividend=pdividendhi));
while (++pdividend <= pdividendlo)
fprintf(stdout,"_%04x", * pdividend);
fprintf(stdout," at %p hi %+ld lo %+ld\n",
dividend, pdividendhi-dividend, pdividendlo-dividend);
#endif
/* Point to the divisor's highest and lowest data words */
pdivisorhi = divisor + (divisor[1] ? 1 : 2);
pdivisorlo = divisor + divisorsize - 1;
trialdivisor = (unsigned long) (* pdivisorhi) + 1;
#ifdef BIGGISHINT_TRACE
fprintf(stdout,"divisor %s%04x", (* divisor & 1)?"-":"",
* (pdivisor=pdivisorhi));
while (++pdivisor <= pdivisorlo)
fprintf(stdout,"_%04x", * pdivisor);
fprintf(stdout," at %p hi %+ld lo %+ld\n",
divisor, pdivisorhi-divisor, pdivisorlo-divisor);
#endif
/* Initialize the quotient (result) */
quotient = (unsigned short *) calloc(bi1size, sizeof(short));
* quotient = bi1size << 1 | sign;
/* Work out at which word in quotient the result will begin */
pquotienthi = quotient + (bi1[1] ? 1 : 2) + (pdivisorlo-pdivisorhi);
pquotientlo = quotient + bi1size - 1;
dividendcarry = 0;
/* Calculate one word of the quotient per loop */
while (pquotienthi<=pquotientlo) {
/* To avoid comparing all the words of the divisor, perform a */
/* trial division of dividendcarry by the first word of the */
/* divisor plus one. */
dividendcarry += * pdividendhi;
trialquotientmin = dividendcarry / trialdivisor;
while (trialquotientmin) {
#ifdef BIGGISHINT_TRACE
fprintf(stdout,"carry(%ld) %x dividendcarry %lx/%lx=%lx\n",
pdividendhi-dividend, * pdividendhi, dividendcarry,
trialdivisor, trialquotientmin);
fprintf(stdout," quotient %s%04x", (* quotient & 1)?"-":"",
* (pquotient=quotient));
while (++pquotient <= pquotientlo)
fprintf(stdout,"_%04x", * pquotient);
fprintf(stdout," at %p hi %+ld lo %+ld\n",
quotient, pquotienthi-quotient, pquotientlo-quotient);
#endif
/* Subtract shifted trialquotient*divisor from dividend. */
tempquotient = trialquotientmin;
divisorshift = 0;
while (tempquotient) {
pdividend = pdividendhi + (pdivisorlo - pdivisorhi) - divisorshift++;
#ifdef BIGGISHINT_TRACE
fprintf(stdout," tempquotient %lx pdividend %+ld\n",
tempquotient, pdividend-dividend);
#endif
/* Use tempquotient1 to subtract successive words of */
/* divisor multiplied by tempquotient2 from dividend. */
quotientcarry = 0UL;
for (pdivisor=pdivisorlo; pdivisor>=pdivisorhi; --pdivisor) {
quotientcarry += (* pdivisor) * (tempquotient & 0xffff);
dividendword1 = (* pdividend);
dividendword2 = quotientcarry & 0xffff;
#ifdef BIGGISHINT_TRACE
fprintf(stdout," dividend[%ld] = %04x-%04x = ",
pdividend-dividend, dividendword1, dividendword2);
#endif
if (dividendword1 >= dividendword2) /* just subtract */
dividendword1 -= dividendword2;
else { /* borrow, then subtract */
dividendword1 = 0x10000UL + dividendword1 - dividendword2;
quotientcarry += 0x10000UL;
}
#ifdef BIGGISHINT_TRACE
fprintf(stdout,"%04x tempquotient1 %lx\n",
dividendword1, quotientcarry);
#endif
* pdividend-- = dividendword1;
quotientcarry>>=16;
}
assert( quotientcarry <= 0xffff ); /* no loop required */
if (quotientcarry) {
#ifdef BIGGISHINT_TRACE
fprintf(stdout," dividend[%ld] = %04x-%04lx = ",
pdividend-dividend, * pdividend, quotientcarry);
#endif
(* pdividend) -= (unsigned short) quotientcarry;
#ifdef BIGGISHINT_TRACE
fprintf(stdout,"%04x\n", * pdividend);
#endif
}
tempquotient >>= 16;
} /* Subtracting shifted trialquotient*divisor from dividend */
/* Add shifted trialquotient to quotient */
tempquotient = trialquotientmin;
pquotient = pquotienthi;
while (tempquotient) {
tempquotient += * pquotient;
* pquotient-- = tempquotient & 0xffff;
tempquotient >>= 16;
}
/* Make a new dividendcarry from dividend. */
/* TODO: use assert() to check whether a loop is really needed */
dividendcarry = 0;
for (pdividend=dividend+1; pdividend<=pdividendhi; ++pdividend ) {
dividendcarry = (dividendcarry<<16) + (* pdividend);
}
trialquotientmin = dividendcarry / trialdivisor;
} /* while (trialquotientmin > 0) */
#ifdef BIGGISHINT_TRACE
fprintf(stdout," dividend %s%04x", (* dividend & 1)?"-":"",
* (pdividend=dividend));
while (++pdividend <= pdividendlo)
fprintf(stdout,"_%04x", * pdividend);
fprintf(stdout," at %p hi %+ld lo %+ld\n",
dividend, pdividendhi-dividend, pdividendlo-dividend);
#endif
/* Proceed to the next word in the quotient and dividend */
dividendcarry <<= 16;
++pquotienthi;
++pdividendhi;
}
#ifdef BIGGISHINT_TRACE
fprintf(stdout," quotient %s%04x", (* quotient & 1)?"-":"",
* (pquotient=quotient));
while (++pquotient <= pquotientlo)
fprintf(stdout,"_%04x", * pquotient);
fprintf(stdout," at %p hi %+ld lo %+ld\n",
quotient, pquotienthi-quotient, pquotientlo-quotient);
#endif
/* Subtract (trialquotient2*divisor) from dividend. */
pdividend = pdividendlo;
#ifdef BIGGISHINT_TRACE
fprintf(stdout," trialquotientmin %lx pdividend %+ld\n",
trialquotientmin, pdividend-dividend);
#endif
/* Try to subtract the divisor from the dividend (this is the */
/* only time divisor is used instead of leading word of divisor */
/* + 1. */
dividendcarry = 0;
while (biggishint_internal_comparemagnitude(dividend, divisor) >= 0) {
pdividend = pdividendlo;
dividendword2 = 0; /* used as the borrow word in subtraction */
for (pdivisor=pdivisorlo; pdivisor>=pdivisorhi; --pdivisor) {
dividendword1 = * pdividend;
if ((unsigned long) dividendword1 >= (unsigned long) * pdivisor + dividendword2) {
/* just subtract */
dividendword1 -= * pdivisor + dividendword2;
dividendword2 = 0;
}
else { /* borrow, then subtract */
dividendword1 = (unsigned long) dividendword1 - * pdivisor + dividendword2;
dividendword2 = 1;
}
* pdividend-- = dividendword1;
}
// unsigned short * dividend_temp = dividend;
// dividend = biggishint_internal_addsubtract(dividend_temp, divisor, 1);
// free(dividend_temp);
++dividendcarry;
}
#ifdef BIGGISHINT_TRACE
fprintf(stdout,"final dividendcarry %lx\n", dividendcarry);
#endif
assert( dividendcarry <= 1 );
pquotient = pquotientlo;
while (dividendcarry && (pquotient>quotient)){
dividendcarry += * pquotient;
* pquotient-- = dividendcarry & 0xffff;
dividendcarry >>= 16;
}
free(dividend);
return biggishint_internal_trim(&quotient);
}
/* biggishintFree */
void
biggishintFree(unsigned short * bi1)
{
free(bi1);
}
/* biggishintFromDecimalString */
/* Note: this algorithm is unacceptably inefficient and should be */
/* rewritten, because the two other biggishint functions it calls */
/* cause two malloc() calls per decimal digit. */
unsigned short *
biggishintFromDecimalString(char * str)
{
char * ps, c;
int sign = 0;
unsigned short * bi1, * bi2;
unsigned short digitvalue[2] = {4,0};
ps = str;
if (* ps == '-') { /* Detect a leading minus sign */
sign = 1;
++ps;
}
/* Create the initial biggishint result with a value of 0 */
bi1 = (unsigned short *) calloc(2,2);
* bi1 = 4 | sign;
/* take one digit at a time, convert to binary, accumulate values */
while ( isdigit(c = * ps++) ) {
biggishint_internal_shortmultiply(&bi1, 10);
digitvalue[1] = c - '0';
bi2 = biggishint_internal_addsubtract(bi1, digitvalue, 0);
free(bi1);
bi1 = bi2;
}
return bi1;
}
/* biggishintFromHexadecimalString */
unsigned short *
biggishintFromHexadecimalString(char * str)
{
int hexdigitcount, biggishintwordcount, i, nybble, sign=0;
unsigned short biggishintarraysize, * biggishint, * shortPointer, value;
char character, * strPointer;
strPointer = str;
if (* strPointer == '-') { /* Detect a leading minus sign */
sign = 1;
++strPointer;
}
if (strncmp(strPointer, "0x", 2) == 0) /* skip the '0x' prefix if it exists */
strPointer += 2;
hexdigitcount = strlen(strPointer);
/* The number of short integers holding values must always be odd */
biggishintwordcount = (((hexdigitcount+3)>>3)<<1)+1; /* 1-4=>1, 5-12=>3 etc */
biggishintarraysize = biggishintwordcount + 1;
biggishint = (unsigned short *) calloc(biggishintarraysize, sizeof(unsigned short));
assert( biggishint != NULL );
shortPointer = biggishint;
* shortPointer++ = (biggishintarraysize << 1) | sign;
/* leave one word blank for 5-8 13-16 21-24 digit strings */
if ( (hexdigitcount-1) & 0x4) ++shortPointer;
value = 0;
for (i=hexdigitcount-1; i>=0; --i) {
character = tolower(* strPointer++);
nybble = character - '0' - (character>='a' ? 'a'-'9'-1 : 0);
value = (value<<4) + nybble;
if ((i%4) == 0) {
* shortPointer++ = value;
value = 0;
}
}
return biggishint;
}
/* biggishintFromLong */
unsigned short *
biggishintFromLong(long l)
{
unsigned short * bi1, negative_bit=0;
if (l<0) {
negative_bit = 1;
l = -l;
}
if (l <= USHRT_MAX) {
bi1 = (unsigned short *) calloc(2,sizeof(short));
* bi1 = 4 | negative_bit;
bi1[1] = (unsigned short) l;
}
else {
bi1 = (unsigned short *) calloc(4,sizeof(short));
* bi1 = 8 | negative_bit;
bi1[3] = (unsigned short) l;
bi1[2] = l >> 16;
}
return bi1;
}
/* biggishintMultiply */
unsigned short *
biggishintMultiply(unsigned short * bi1, unsigned short * bi2)
{
unsigned short bi1size, bi2size, res1size, res2size, sign1, sign2;
unsigned short * p1, * p2, * result, * presult;
int i1, i2;
unsigned long resultsize, n1, n2, subtotal, carry;
/* Before starting on the main long multiplication, which is */
/* slow, try to identify multipliers that offer an opportunity */
/* for a shortcut, for example by 0 or 1, shifting left for */
/* multipliers that are multiples of powers of two, or short */
/* multiplication. */
bi1size = (* bi1 & 0xfffc) >> 1;
bi2size = (* bi2 & 0xfffc) >> 1;
if (bi1size == 2) {
result = biggishint_internal_clone(bi2);
* result &= 0xfffe; /* clear the sign bit */
biggishint_internal_shortmultiply(&result, bi1[1]);
}
else {
if (bi2size == 2) {
result = biggishint_internal_clone(bi1);
* result &= 0xfffe; /* clear the sign bit */
biggishint_internal_shortmultiply(&result, bi2[1]);
}
else { /* both bi1 and bi2 are more than 16 bit numbers */
/* Create a result array that is large enough for any */
/* possible product. First calculate the smallest size */
/* according to the contents of bi1 and bi2, regardless */
/* of the need to round up to an even number. */
res1size = bi1size + (bi1[1] ? 0 : -1); /* the first word may be 0 */
res2size = bi2size + (bi2[1] ? 0 : -1);
/* Then add them together and round to an even number */
resultsize = (res1size + res2size + 1) & 0xfffe;
result = (unsigned short *) calloc(resultsize, sizeof(short));
* result = (resultsize << 1);
presult = result + resultsize;
p1 = bi1 + bi1size;
for (i1=1; i1<bi1size; ++i1) {
n1 = * --p1;
presult = result + resultsize - i1;
p2 = bi2 + bi2size;
carry = 0;
for (i2=1; i2<bi2size; ++i2) {
n2 = * --p2;
subtotal = (* presult) + n1 * n2 + carry;
carry = subtotal >> 16;
* presult-- = subtotal & 0xffff;
}
while (carry) {
* presult-- = carry;
carry >>= 16;
}
}
}
}
sign1 = * bi1 & 1;
sign2 = * bi2 & 1;
* result |= sign1 ^ sign2;
return biggishint_internal_trim(&result);
}
/* biggishintShiftLeft */
unsigned short *
biggishintShiftLeft(unsigned short * bi1, unsigned short * bi2)
{ /* TODO: shift counts from 65536 to 524239 and -1 to -524239 bits */
return biggishint_internal_shiftleft(bi1, bi2[1]);
}
unsigned short *
biggishintShiftRight(unsigned short * bi1, unsigned short * bi2)
{
return NULL;
}
/* biggishintSubtract */
unsigned short *
biggishintSubtract(unsigned short * bi1, unsigned short * bi2)
{
return biggishint_internal_addsubtract(bi1, bi2, 1);
}
/* biggishintToDecimalString */
char *
biggishintToDecimalString(unsigned short * bi1)
{
/* The number of decimal digits that will be created is difficult */
/* (or slow) to calculate in advance. This routine initially */
/* over-allocates memory, and then sizes it correctly at the end. */
unsigned short bi1size, * bi2, digit, sign1;
int strsize, leadingzeroes;
char * result, * pdigits, * p1;
/* Calculate the very maximum number of characters that the */
/* resulting string can occupy, including a terminating '\0'. */
/* Each word is '65535' at most, then '\0' */
bi2 = biggishint_internal_clone(bi1);
bi1size = (* bi1 & 0xfffc) >> 1;
sign1 = * bi1 & 1;
strsize = (bi1size-1) * 5 + sign1 + 1;
result = (char *) malloc(strsize);
assert( result != NULL );
pdigits = result;
if (sign1) * pdigits++ = '-';
p1 = result + strsize;
(* --p1) = '\0';
do {
digit = 10;
biggishint_internal_shortdivide(bi2, &digit);
(* --p1) = '0' + digit;
} while ( p1 > pdigits ); /* TODO: find other ways to finish early */
free(bi2);
/* Count how many '0' characters there are at the beginning of */
/* the string, and then move the non '0' characters. */
leadingzeroes = strspn(pdigits, "0");
if (leadingzeroes == strsize - sign1 - 1)
--leadingzeroes;
if (leadingzeroes) {
memmove(pdigits, pdigits+leadingzeroes, strsize-sign1-leadingzeroes); /* (Big Endian)-- ;) */
result = realloc(result, strsize-leadingzeroes);
}
return result;
}
/* biggishintToHexadecimalString */
char *
biggishintToHexadecimalString(unsigned short * bi1)
{
int bi1size, hexstringsize, i, j, value, nybble, emitzero, sign;
char * hexString, * hexPointer;
bi1size = (* bi1 & 0xfffc) >> 1;
sign = * bi1 & 1;
/* Calculate how many characters the hex string needs, including */
/* the "0x" at the beginning and a '\0' at the end */
hexstringsize = (biggishint_internal_bitsize(bi1[1])+3)>>2; /* 0=>0, 1-4=>1, 5-8=>2 */
hexstringsize += (hexstringsize?0:1) /* allow for 0 digit */
+ 3 + sign + ((bi1size-2) << 2); /* '0x' + sign + digits + '\0' */
hexString = (char *) malloc(hexstringsize);
assert( hexString != NULL );
hexPointer = hexString;
if (sign) * hexPointer++ = '-';
* hexPointer++ = '0'; * hexPointer++ = 'x';
emitzero = 0; /* do not emit leading zeroes */
for (i=1; i<bi1size; ++i) {
value = bi1[i];
for (j=3; j>=0; --j) {
nybble = (value >> (j*4)) & 0xf;
if (nybble || emitzero) {
* hexPointer++ = '0' + nybble + ((nybble>9) ? 'a'-'9'-1 : 0);
emitzero = 1;
}
}
}
if (! emitzero)
* hexPointer++ = '0';
* hexPointer = '\0';
return hexString;
}
/* ----------------------- Internal functions ----------------------- */
/* Except for biggishint_internal_trim, the internal functions do not */
/* trim their results, because it costs time, may be redundant, and */
/* increases heap churn. */
/* biggishint_internal_addsubtract */
unsigned short * biggishint_internal_addsubtract(unsigned short * bi1,
unsigned short * bi2, int flipsign2)
{
unsigned short bi1size, bi2size, res1size, res2size, resultsize;
unsigned short * result1, * result2, * larger, * smaller, * p1, * p2;
unsigned int sign1, sign2, sign, carry, i1, i2;
signed long partialresult;
sign1 = * bi1 & 1;
sign2 = (* bi2 & 1) ^ flipsign2;
if (sign1 ^ sign2) { /* different signs, do a subtract */
/* the larger number determines the size and sign of the result */
if (biggishint_internal_comparemagnitude(bi1,bi2) >= 0) {
larger = bi1; smaller = bi2; sign = sign1;
}
else {
smaller = bi1; larger = bi2; sign = sign2;
}
resultsize = (* larger & 0xfffc) >> 1;
result1 = (unsigned short *) calloc(resultsize, sizeof(short));
result2 = result1 + resultsize;
p1 = larger + ((* larger & 0xfffc) >> 1);
p2 = smaller + ((* smaller & 0xfffc) >> 1);
carry = 0;
partialresult = 0;
while (--p1 > larger) {
i1 = * p1;
i2 = (--p2 > smaller) ? * p2 : 0;
partialresult += i1 - i2;
* --result2 = (partialresult >=0)
? (unsigned short) partialresult
: (unsigned short) (partialresult + 65536);
partialresult = (partialresult & 0xffff0000) ? -1 : 0;
}
} /* subtract */
else { /* same signs, do an add */
bi1size = (* bi1 & 0xfffc) >> 1;
bi2size = (* bi2 & 0xfffc) >> 1;
res1size = bi1size + (bi1[1] ? 1 : 0); /* the first word may be 0 */
res2size = bi2size + (bi2[1] ? 1 : 0);
resultsize = ((res1size > res2size ? res1size : res2size) + 1) & 0xfffe;
sign = sign1;
result1 = (unsigned short *) calloc(resultsize, sizeof(short));
assert( result1 != NULL );
/* Initialize pointers to the augend (bi1), addend (bi2) and partialresult */
* result1 = resultsize << 1;
p1 = bi1 + bi1size - 1;
p2 = bi2 + bi2size - 1;
result2 = result1 + resultsize - 1;
carry = 0;
/* Iteratively add words from least significant to most */
while (result2 > result1) {
i1 = i2 = 0;
if (p1 > bi1)
i1 = * p1--;
if (p2 > bi2)
i2 = * p2--;
partialresult = i1 + i2 + carry;
carry = 0;
if (partialresult > USHRT_MAX) {
carry = 1;
partialresult -= (USHRT_MAX + 1);
}
* result2-- = (unsigned short) partialresult;
} /* while */
} /* add */
* result1 = (resultsize << 1) | sign;
return biggishint_internal_trim(&result1);
}
/* biggishint_internal_bitsize */
/* Count how many bits a number uses (0-64), returns 1 + position of first 1 bit */
int
biggishint_internal_bitsize(unsigned long n)
{
int bitsize = 0;
for ( ; n; n >>= 1)
++bitsize;
return bitsize;
}
/* biggishint_internal_clone */
unsigned short *
biggishint_internal_clone(unsigned short * bi1)
{
unsigned short clonebytes, * clone;
clonebytes = * bi1 & 0xfffffffe;
clone = (unsigned short *) malloc(clonebytes);
assert( clone != NULL );
memcpy(clone, bi1, clonebytes);
return clone;
}
/* biggishint_internal_comparemagnitude */
/* returns -1 if bi1<bi2, 0 if bi1==bi2, +1 if bi1>bi2 */
int biggishint_internal_comparemagnitude(unsigned short * bi1, unsigned short * bi2)
{
unsigned short * pi1data, * pi1, * pi2data, * pi2;
unsigned short i1, i2, bi1size, bi2size, loopcount;
int result=0;
/* This function could often be quicker by comparing the sizes of */
/* the two numbers, but that implies trusting the rest of the */
/* code to always trim leading zero words where possible. The */
/* test suite currently lacks the coverage required to enable */
/* that trust. */
bi1size = (* bi1 & 0xfffc) >> 1;
bi2size = (* bi2 & 0xfffc) >> 1;
/* the max number of comparisons is max(bi1size,bi2size)-1 */
loopcount = ((bi1size>bi2size) ? bi1size : bi2size) - 1;
pi1data = bi1 + 1;
pi2data = bi2 + 1;
pi1 = bi1 + bi1size - loopcount;
pi2 = bi2 + bi2size - loopcount;
while (result==0 && loopcount--) {
/* substitute leading zeroes for whichever number is shorter */
i1 = (pi1<pi1data) ? 0 : * pi1; ++pi1;
i2 = (pi2<pi2data) ? 0 : * pi2; ++pi2;
/* compare the two words from each biggishint */
result = (i1<i2) ? -1 : (i1>i2) ? 1 : 0;
}
return result;
}
/* biggishint_internal_shiftleft */
unsigned short *
biggishint_internal_shiftleft(unsigned short * bi1, unsigned int bitcount)
{
unsigned short bi1size, * result, * p1, * p2, carry;
unsigned int inputbitcount, resultbitcount, resultsize, inputloop;
unsigned int shiftleft, shiftright, inputword, resultword;
bi1size = (* bi1 & 0xfffc) >> 1;
/* Calculate the number of data bits needed for the result */
inputbitcount = biggishint_internal_bitsize(bi1[1]) + ((bi1size - 2) << 4);
resultbitcount = inputbitcount + bitcount;
assert( resultbitcount < 524272);
/* Calculate the total number of words to allocate for the result */
resultsize = 1 + ((resultbitcount + 15) >> 4);
result = (unsigned short *) calloc(resultsize, sizeof(short));
* result = resultsize << 1;
/* prepare initial values for the loop below */
shiftleft = bitcount & 0xf;
shiftright = 16 - shiftleft;
p1 = bi1 + 1;
p2 = result + 1;
inputword = * p1++;
inputloop = bi1size;
carry = inputword << shiftleft;
/* Check whether part of the first word of the input needs to */
/* carry over to the second word of the result */
if ( ((inputbitcount-1) & 0xf) > ((resultbitcount-1) & 0xf) ) {
/* Yes, so here store the bits that will not be carried */
resultword = inputword >> shiftright;
* p2++ = resultword;
--inputloop;
}
while (--inputloop) {
inputword = * p1++;
resultword = carry | (inputword >> shiftright);
* p2++ = resultword;
carry = inputword << shiftleft;
}
* p2++ = carry;
return result;
}
/* TODO: biggishint_internal_shiftright */
unsigned short *
biggishint_internal_shiftright(unsigned short * bi1, unsigned int bitcount)
{
unsigned short bi1size, * result;
bi1size = 2;
result = (unsigned short *) malloc(bi1size);
return result;
}
/* biggishint_internal_shortdivide */
/* Short division only (divisor <= 0xffff). */
/* Returns quotient in (* bi1), remainder in (* i2) */
void
biggishint_internal_shortdivide(unsigned short * bi1, unsigned short * i2)
{
unsigned short bi1size, * pi1, * pi2, divisor, remainder, hi, lo;
unsigned long partialdividend, partialquotient;
bi1size = (* bi1 & 0xfffc) >> 1;
divisor = * i2;
remainder = 0;
pi1 = bi1 + 1;
pi2 = bi1 + bi1size;
lo = 0;
while ( pi1 < pi2 ) {
hi = lo;
lo = * pi1;
partialdividend = ((unsigned long)hi << 16) | lo;
partialquotient = partialdividend / divisor;
remainder = partialdividend % divisor;
hi = partialquotient & 0xffff;
lo = remainder;
* pi1++ = hi;
}
* i2 = remainder;
}
/* biggishint_internal_shortmultiply */
void
biggishint_internal_shortmultiply(unsigned short ** bi1, unsigned short multiplier)
{
unsigned short bi1size, productsize, * product, * pi, * pp;
unsigned long productcarry;
if (multiplier == 0) {
* bi1 = realloc(* bi1, sizeof(short)<<1);
(* bi1)[0] = 4; (* bi1)[1] = 0;
}
else {
if (multiplier != 1) {
/* TODO: avoid realloc if possible */
bi1size = (** bi1 & 0xfffc) >> 1;
productsize = bi1size + 2; /* even number of words */
product = (unsigned short *) calloc(productsize, sizeof(short));
* product = productsize << 1;
pi = * bi1 + bi1size - 1;
pp = product + productsize - 1;
productcarry = 0;
do {
productcarry += (* pi--) * (unsigned int)multiplier;
(* pp--) = (unsigned short) productcarry;
productcarry >>= 16;
} while ( pi > * bi1 );
assert( productcarry < 0xffffffff );
while (productcarry) {
(* pp--) = (unsigned short) productcarry;
productcarry >>= 16;
}
free(* bi1);
* bi1 = product;
biggishint_internal_trim(bi1);
}
}
}
/*
#include <stdio.h>
fprintf(stderr,"sign %d\n", sign);
*/
/* biggishint_internal_trim */
/* If possible, remove leading zeroes from the front of the biggishint */
/* Also remove the minus sign from -0 results */
/* Note: it reallocates the data, so it may be at a new address. */
unsigned short *
biggishint_internal_trim(unsigned short ** pbi1)
{
/* For example, change this: */
/* +------+------+------+------+------+------+ */
/* | 000c | 0000 | 0000 | 0000 | 1234 | cdef | */
/* +------+------+------+------+------+------+ */
/* to this: */
/* +------+------+------+------+ */
/* | 0008 | 0000 | 1234 | cdef | */
/* +------+------+------+------+ */
unsigned int sign;
unsigned short bi1size, newsize;
unsigned short * bi1, * pLeft, * pSearch, * pRight, * pAfterZeroes;
bi1 = * pbi1;
/* Count the number of contiguous leading zero words */
bi1size = (* bi1 & 0xfffc) >> 1;
sign = * bi1 & 1;
pLeft = bi1 + 1;
pSearch = bi1;
pRight = bi1 + bi1size; /* just outside the biggishint */
pAfterZeroes = pRight;
/* Try to set pSearch to the address of the first non zero word. */
/* After this loop completes, pAfterZeroes either points to the */
/* first nonzero word, or to the first address after the biggishint. */
while (++pSearch < pAfterZeroes) /* note pAfterZeroes moves! */
if ( * pSearch )
pAfterZeroes = pSearch;
/* If there are leading words filled with zeroes, move the non */
/* zero words to the left to overwrite them */
if (pAfterZeroes > pLeft) {
newsize = (pRight - pAfterZeroes + 2) & 0xfffe; /* always even */
if (newsize < bi1size) {
/* Trim the size of the memory allocation */
/* Bump the destination by 1 if the first non zero word */
/* was at an even subscript */
pLeft += (pAfterZeroes - bi1 + 1) & 1;
/* If the array was little endian, memmove would not happen */
memmove(pLeft, pAfterZeroes, (pRight - pAfterZeroes) << 1);
bi1 = realloc(bi1, newsize << 1);
* bi1 = (newsize << 1) | sign;
* pbi1 = bi1;
}
}
/* Convert the silly case of -0 into 0 */
if (* bi1 == 5 && bi1[1]==0) * bi1 = 4;
return bi1;
}
/* end of biggishint.c */
Jump to Line
Something went wrong with that request. Please try again.