ossp-pkg/tai/tai_ui64.c
1.3
/*
** OSSP ui64 - 64-Bit Arithmetic
** Copyright (c) 2002-2005 Ralf S. Engelschall <rse@engelschall.com>
** Copyright (c) 2002-2005 The OSSP Project <http://www.ossp.org/>
**
** This file is part of OSSP ui64, a 64-bit arithmetic library
** which can be found at http://www.ossp.org/pkg/lib/ui64/.
**
** Permission to use, copy, modify, and distribute this software for
** any purpose with or without fee is hereby granted, provided that
** the above copyright notice and this permission notice appear in all
** copies.
**
** THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
** WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
** MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
** IN NO EVENT SHALL THE AUTHORS AND COPYRIGHT HOLDERS AND THEIR
** CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
** SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
** LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
** USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
** ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
** OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
** SUCH DAMAGE.
**
** ui64.c: implementation of 64-bit unsigned integer arithmetic
*/
#include <string.h>
#include <ctype.h>
#include "tai_ui64.h"
#define UI64_BASE 256 /* 2^8 */
#define UI64_DIGITS 8 /* 8*8 = 64 bit */
#define UIXX_T(n) struct { unsigned char x[n]; }
/* fill an ui64_t with a sequence of a particular digit */
#define ui64_fill(__x, __n) \
do { int __i; \
for (__i = 0; __i < UI64_DIGITS; __i++) \
(__x).x[__i] = (__n); \
} while(0)
/* the value zero */
ui64_t ui64_zero(void)
{
ui64_t z;
ui64_fill(z, 0);
return z;
}
/* the maximum value */
ui64_t ui64_max(void)
{
ui64_t z;
ui64_fill(z, UI64_BASE-1);
return z;
}
/* convert ISO-C "unsigned long" into internal format */
ui64_t ui64_n2i(unsigned long n)
{
ui64_t z;
int i;
i = 0;
do {
z.x[i++] = (n % UI64_BASE);
} while ((n /= UI64_BASE) > 0 && i < UI64_DIGITS);
for ( ; i < UI64_DIGITS; i++)
z.x[i] = 0;
return z;
}
/* convert internal format into ISO-C "unsigned long";
truncates if sizeof(unsigned long) is less than UI64_DIGITS! */
unsigned long ui64_i2n(ui64_t x)
{
unsigned long n;
int i;
n = 0;
i = (int)sizeof(n);
if (i > UI64_DIGITS)
i = UI64_DIGITS;
while (--i >= 0) {
n = (n * UI64_BASE) + x.x[i];
}
return n;
}
/* convert string representation of arbitrary base into internal format */
ui64_t ui64_s2i(const char *str, char **end, int base)
{
ui64_t z;
const char *cp;
int carry;
static char map[] = {
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, /* 0...9 */
36, 36, 36, 36, 36, 36, 36, /* illegal chars */
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, /* A...M */
23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, /* N...Z */
36, 36, 36, 36, 36, 36, /* illegal chars */
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, /* a...m */
23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35 /* m...z */
};
ui64_fill(z, 0);
if (str == NULL || (base < 2 || base > 36))
return z;
cp = str;
while (*cp != '\0' && isspace((int)(*cp)))
cp++;
while ( *cp != '\0'
&& isalnum((int)(*cp))
&& map[(int)(*cp)-'0'] < base) {
z = ui64_muln(z, base, &carry);
if (carry)
break;
z = ui64_addn(z, map[(int)(*cp)-'0'], &carry);
if (carry)
break;
cp++;
}
if (end != NULL)
*end = (char *)cp;
return z;
}
/* convert internal format into string representation of arbitrary base */
char *ui64_i2s(ui64_t x, char *str, size_t len, int base)
{
static char map[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
char c;
int r;
int n;
int i, j;
if (str == NULL || len < 2 || (base < 2 || base > 36))
return NULL;
n = ui64_len(x);
i = 0;
do {
x = ui64_divn(x, base, &r);
str[i++] = map[r];
while (n > 1 && x.x[n-1] == 0)
n--;
} while (i < ((int)len-1) && (n > 1 || x.x[0] != 0));
str[i] = '\0';
for (j = 0; j < --i; j++) {
c = str[j];
str[j] = str[i];
str[i] = c;
}
return str;
}
/* addition of two ui64_t */
ui64_t ui64_add(ui64_t x, ui64_t y, ui64_t *ov)
{
ui64_t z;
int carry;
int i;
carry = 0;
for (i = 0; i < UI64_DIGITS; i++) {
carry += (x.x[i] + y.x[i]);
z.x[i] = (carry % UI64_BASE);
carry /= UI64_BASE;
}
if (ov != NULL)
*ov = ui64_n2i((unsigned long)carry);
return z;
}
/* addition of an ui64_t and a single digit */
ui64_t ui64_addn(ui64_t x, int y, int *ov)
{
ui64_t z;
int i;
for (i = 0; i < UI64_DIGITS; i++) {
y += x.x[i];
z.x[i] = (y % UI64_BASE);
y /= UI64_BASE;
}
if (ov != NULL)
*ov = y;
return z;
}
/* subtraction of two ui64_t */
ui64_t ui64_sub(ui64_t x, ui64_t y, ui64_t *ov)
{
ui64_t z;
int borrow;
int i;
int d;
borrow = 0;
for (i = 0; i < UI64_DIGITS; i++) {
d = ((x.x[i] + UI64_BASE) - borrow - y.x[i]);
z.x[i] = (d % UI64_BASE);
borrow = (1 - (d/UI64_BASE));
}
if (ov != NULL)
*ov = ui64_n2i((unsigned long)borrow);
return z;
}
/* subtraction of an ui64_t and a single digit */
ui64_t ui64_subn(ui64_t x, int y, int *ov)
{
ui64_t z;
int i;
int d;
for (i = 0; i < UI64_DIGITS; i++) {
d = (x.x[i] + UI64_BASE) - y;
z.x[i] = (d % UI64_BASE);
y = (1 - (d/UI64_BASE));
}
if (ov != NULL)
*ov = y;
return z;
}
/*
7 3 2
* 9 4 2 8
---------
5 8 5 6
+ 1 4 6 4
+ 2 9 2 8
+ 6 5 8 8
---------------
= 6 9 0 1 2 9 6
*/
ui64_t ui64_mul(ui64_t x, ui64_t y, ui64_t *ov)
{
UIXX_T(UI64_DIGITS+UI64_DIGITS) zx;
ui64_t z;
int carry;
int i, j;
/* clear temporary result buffer */
for (i = 0; i < (UI64_DIGITS+UI64_DIGITS); i++)
zx.x[i] = 0;
/* perform multiplication operation */
for (i = 0; i < UI64_DIGITS; i++) {
/* calculate partial product and immediately add to z */
carry = 0;
for (j = 0; j < UI64_DIGITS; j++) {
carry += (x.x[i] * y.x[j]) + zx.x[i+j];
zx.x[i+j] = (carry % UI64_BASE);
carry /= UI64_BASE;
}
/* add carry to remaining digits in z */
for ( ; j < UI64_DIGITS + UI64_DIGITS - i; j++) {
carry += zx.x[i+j];
zx.x[i+j] = (carry % UI64_BASE);
carry /= UI64_BASE;
}
}
/* provide result by splitting zx into z and ov */
memcpy(z.x, zx.x, UI64_DIGITS);
if (ov != NULL)
memcpy((*ov).x, &zx.x[UI64_DIGITS], UI64_DIGITS);
return z;
}
ui64_t ui64_muln(ui64_t x, int y, int *ov)
{
ui64_t z;
int carry;
int i;
carry = 0;
for (i = 0; i < UI64_DIGITS; i++) {
carry += (x.x[i] * y);
z.x[i] = (carry % UI64_BASE);
carry /= UI64_BASE;
}
if (ov != NULL)
*ov = carry;
return z;
}
/*
= 2078 [q]
0615367 [x] : 296 [y]
-0592 [dq]
-----
= 0233
-0000 [dq]
-----
= 2336
-2072 [dq]
-----
= 2647
-2308 [dq]
-----
= 279 [r]
*/
ui64_t ui64_div(ui64_t x, ui64_t y, ui64_t *ov)
{
ui64_t q;
ui64_t r;
int i;
int n, m;
int ovn;
/* determine actual number of involved digits */
n = ui64_len(x);
m = ui64_len(y);
if (m == 1) {
/* simple case #1: reduceable to ui64_divn() */
if (y.x[0] == 0) {
/* error case: division by zero! */
ui64_fill(q, 0);
ui64_fill(r, 0);
}
else {
q = ui64_divn(x, y.x[0], &ovn);
ui64_fill(r, 0);
r.x[0] = (unsigned char)ovn;
}
} else if (n < m) {
/* simple case #2: everything is in the remainder */
ui64_fill(q, 0);
r = x;
} else { /* n >= m, m > 1 */
/* standard case: x[0..n] / y[0..m] */
UIXX_T(UI64_DIGITS+1) rx;
UIXX_T(UI64_DIGITS+1) dq;
ui64_t t;
int km;
int k;
int qk;
unsigned long y2;
unsigned long r3;
int borrow;
int d;
/* rx is x with a leading zero in order to make
sure that n > m and not just n >= m */
memcpy(rx.x, x.x, UI64_DIGITS);
rx.x[UI64_DIGITS] = 0;
for (k = n - m; k >= 0; k--) {
/* efficiently compute qk by guessing
qk := rx[k+m-2...k+m]/y[m-2...m-1] */
km = k + m;
y2 = (y.x[m-1]*UI64_BASE) + y.x[m-2];
r3 = (rx.x[km]*(UI64_BASE*UI64_BASE)) +
(rx.x[km-1]*UI64_BASE) + rx.x[km-2];
qk = r3 / y2;
if (qk >= UI64_BASE)
qk = UI64_BASE - 1;
/* dq := y*qk (post-adjust qk if guessed incorrectly) */
t = ui64_muln(y, qk, &ovn);
memcpy(dq.x, t.x, UI64_DIGITS);
dq.x[m] = (unsigned char)ovn;
for (i = m; i > 0; i--)
if (rx.x[i+k] != dq.x[i])
break;
if (rx.x[i+k] < dq.x[i]) {
t = ui64_muln(y, --qk, &ovn);
memcpy(dq.x, t.x, UI64_DIGITS);
dq.x[m] = (unsigned char)ovn;
}
/* store qk */
q.x[k] = (unsigned char)qk;
/* rx := rx - dq*(b^k) */
borrow = 0;
for (i = 0; i < m+1; i++) {
d = ((rx.x[k+i] + UI64_BASE) - borrow - dq.x[i]);
rx.x[k+i] = (d % UI64_BASE);
borrow = (1 - (d/UI64_BASE));
}
}
memcpy(r.x, rx.x, m);
/* fill out results with leading zeros */
for (i = n-m+1; i < UI64_DIGITS; i++)
q.x[i] = 0;
for (i = m; i < UI64_DIGITS; i++)
r.x[i] = 0;
}
/* provide results */
if (ov != NULL)
*ov = r;
return q;
}
ui64_t ui64_divn(ui64_t x, int y, int *ov)
{
ui64_t z;
unsigned int carry;
int i;
carry = 0;
for (i = (UI64_DIGITS - 1); i >= 0; i--) {
carry = (carry * UI64_BASE) + x.x[i];
z.x[i] = (carry / y);
carry %= y;
}
if (ov != NULL)
*ov = carry;
return z;
}
ui64_t ui64_and(ui64_t x, ui64_t y)
{
ui64_t z;
int i;
for (i = 0; i < UI64_DIGITS; i++)
z.x[i] = (x.x[i] & y.x[i]);
return z;
}
ui64_t ui64_or(ui64_t x, ui64_t y)
{
ui64_t z;
int i;
for (i = 0; i < UI64_DIGITS; i++)
z.x[i] = (x.x[i] | y.x[i]);
return z;
}
ui64_t ui64_xor(ui64_t x, ui64_t y)
{
ui64_t z;
int i;
for (i = 0; i < UI64_DIGITS; i++)
z.x[i] = ((x.x[i] & ~(y.x[i])) | (~(x.x[i]) & (y.x[i])));
return z;
}
ui64_t ui64_not(ui64_t x)
{
ui64_t z;
int i;
for (i = 0; i < UI64_DIGITS; i++)
z.x[i] = ~(x.x[i]);
return z;
}
ui64_t ui64_rol(ui64_t x, int s, ui64_t *ov)
{
UIXX_T(UI64_DIGITS+UI64_DIGITS) zx;
ui64_t z;
int i;
int carry;
if (s <= 0) {
/* no shift at all */
if (ov != NULL)
*ov = ui64_zero();
return x;
}
else if (s > 64) {
/* too large shift */
if (ov != NULL)
*ov = ui64_zero();
return ui64_zero();
}
else if (s == 64) {
/* maximum shift */
if (ov != NULL)
*ov = x;
return ui64_zero();
}
else { /* regular shift */
/* shift (logically) left by s/8 bytes */
for (i = 0; i < UI64_DIGITS+UI64_DIGITS; i++)
zx.x[i] = 0;
for (i = 0; i < UI64_DIGITS; i++)
zx.x[i+(s/8)] = x.x[i];
/* shift (logically) left by remaining s%8 bits */
s %= 8;
if (s > 0) {
carry = 0;
for (i = 0; i < UI64_DIGITS+UI64_DIGITS; i++) {
carry += (zx.x[i] * (1 << s));
zx.x[i] = (carry % UI64_BASE);
carry /= UI64_BASE;
}
}
memcpy(z.x, zx.x, UI64_DIGITS);
if (ov != NULL)
memcpy((*ov).x, &zx.x[UI64_DIGITS], UI64_DIGITS);
}
return z;
}
ui64_t ui64_ror(ui64_t x, int s, ui64_t *ov)
{
UIXX_T(UI64_DIGITS+UI64_DIGITS) zx;
ui64_t z;
int i;
int carry;
if (s <= 0) {
/* no shift at all */
if (ov != NULL)
*ov = ui64_zero();
return x;
}
else if (s > 64) {
/* too large shift */
if (ov != NULL)
*ov = ui64_zero();
return ui64_zero();
}
else if (s == 64) {
/* maximum shift */
if (ov != NULL)
*ov = x;
return ui64_zero();
}
else { /* regular shift */
/* shift (logically) right by s/8 bytes */
for (i = 0; i < UI64_DIGITS+UI64_DIGITS; i++)
zx.x[i] = 0;
for (i = 0; i < UI64_DIGITS; i++)
zx.x[UI64_DIGITS+i-(s/8)] = x.x[i];
/* shift (logically) right by remaining s%8 bits */
s %= 8;
if (s > 0) {
carry = 0;
for (i = (UI64_DIGITS+UI64_DIGITS - 1); i >= 0; i--) {
carry = (carry * UI64_BASE) + zx.x[i];
zx.x[i] = (carry / (1 << s));
carry %= (1 << s);
}
}
memcpy(z.x, &zx.x[UI64_DIGITS], UI64_DIGITS);
if (ov != NULL)
memcpy((*ov).x, zx.x, UI64_DIGITS);
}
return z;
}
int ui64_cmp(ui64_t x, ui64_t y)
{
int i;
i = UI64_DIGITS - 1;
while (i > 0 && x.x[i] == y.x[i])
i--;
return (x.x[i] - y.x[i]);
}
int ui64_len(ui64_t x)
{
int i;
for (i = UI64_DIGITS; i > 1 && x.x[i-1] == 0; i--)
;
return i;
}
