diff options
Diffstat (limited to 'ocaml-binary-annot/otherlibs/num/bng.c')
| -rw-r--r-- | ocaml-binary-annot/otherlibs/num/bng.c | 435 |
1 files changed, 435 insertions, 0 deletions
diff --git a/ocaml-binary-annot/otherlibs/num/bng.c b/ocaml-binary-annot/otherlibs/num/bng.c new file mode 100644 index 0000000000..7450247009 --- /dev/null +++ b/ocaml-binary-annot/otherlibs/num/bng.c @@ -0,0 +1,435 @@ +/***********************************************************************/ +/* */ +/* Objective Caml */ +/* */ +/* Xavier Leroy, projet Cristal, INRIA Rocquencourt */ +/* */ +/* Copyright 2003 Institut National de Recherche en Informatique et */ +/* en Automatique. All rights reserved. This file is distributed */ +/* under the terms of the GNU Library General Public License, with */ +/* the special exception on linking described in file ../../LICENSE. */ +/* */ +/***********************************************************************/ + +/* $Id$ */ + +#include "bng.h" +#include "config.h" + +#if defined(__GNUC__) && BNG_ASM_LEVEL > 0 +#if defined(BNG_ARCH_ia32) +#include "bng_ia32.c" +#elif defined(BNG_ARCH_amd64) +#include "bng_amd64.c" +#elif defined(BNG_ARCH_ppc) +#include "bng_ppc.c" +#elif defined (BNG_ARCH_alpha) +#include "bng_alpha.c" +#elif defined (BNG_ARCH_sparc) +#include "bng_sparc.c" +#elif defined (BNG_ARCH_mips) +#include "bng_mips.c" +#endif +#endif + +#include "bng_digit.c" + +/**** Operations that cannot be overridden ****/ + +/* Return number of leading zero bits in d */ +int bng_leading_zero_bits(bngdigit d) +{ + int n = BNG_BITS_PER_DIGIT; +#ifdef ARCH_SIXTYFOUR + if ((d & 0xFFFFFFFF00000000L) != 0) { n -= 32; d = d >> 32; } +#endif + if ((d & 0xFFFF0000) != 0) { n -= 16; d = d >> 16; } + if ((d & 0xFF00) != 0) { n -= 8; d = d >> 8; } + if ((d & 0xF0) != 0) { n -= 4; d = d >> 4; } + if ((d & 0xC) != 0) { n -= 2; d = d >> 2; } + if ((d & 2) != 0) { n -= 1; d = d >> 1; } + return n - d; +} + +/* Complement the digits of {a,len} */ +void bng_complement(bng a/*[alen]*/, bngsize alen) +{ + for (/**/; alen > 0; alen--, a++) *a = ~*a; +} + +/* Return number of significant digits in {a,alen}. */ +bngsize bng_num_digits(bng a/*[alen]*/, bngsize alen) +{ + while (1) { + if (alen == 0) return 1; + if (a[alen - 1] != 0) return alen; + alen--; + } +} + +/* Return 0 if {a,alen} = {b,blen} + -1 if {a,alen} < {b,blen} + 1 if {a,alen} > {b,blen}. */ +int bng_compare(bng a/*[alen]*/, bngsize alen, + bng b/*[blen]*/, bngsize blen) +{ + bngdigit da, db; + + while (alen > 0 && a[alen-1] == 0) alen--; + while (blen > 0 && b[blen-1] == 0) blen--; + if (alen > blen) return 1; + if (alen < blen) return -1; + while (alen > 0) { + alen--; + da = a[alen]; + db = b[alen]; + if (da > db) return 1; + if (da < db) return -1; + } + return 0; +} + +/**** Generic definitions of the overridable operations ****/ + +/* {a,alen} := {a, alen} + carry. Return carry out. */ +static bngcarry bng_generic_add_carry + (bng a/*[alen]*/, bngsize alen, bngcarry carry) +{ + if (carry == 0 || alen == 0) return carry; + do { + if (++(*a) != 0) return 0; + a++; + } while (--alen); + return 1; +} + +/* {a,alen} := {a,alen} + {b,blen} + carry. Return carry out. + Require alen >= blen. */ +static bngcarry bng_generic_add + (bng a/*[alen]*/, bngsize alen, + bng b/*[blen]*/, bngsize blen, + bngcarry carry) +{ + alen -= blen; + for (/**/; blen > 0; blen--, a++, b++) { + BngAdd2Carry(*a, carry, *a, *b, carry); + } + if (carry == 0 || alen == 0) return carry; + do { + if (++(*a) != 0) return 0; + a++; + } while (--alen); + return 1; +} + +/* {a,alen} := {a, alen} - carry. Return carry out. */ +static bngcarry bng_generic_sub_carry + (bng a/*[alen]*/, bngsize alen, bngcarry carry) +{ + if (carry == 0 || alen == 0) return carry; + do { + if ((*a)-- != 0) return 0; + a++; + } while (--alen); + return 1; +} + +/* {a,alen} := {a,alen} - {b,blen} - carry. Return carry out. + Require alen >= blen. */ +static bngcarry bng_generic_sub + (bng a/*[alen]*/, bngsize alen, + bng b/*[blen]*/, bngsize blen, + bngcarry carry) +{ + alen -= blen; + for (/**/; blen > 0; blen--, a++, b++) { + BngSub2Carry(*a, carry, *a, *b, carry); + } + if (carry == 0 || alen == 0) return carry; + do { + if ((*a)-- != 0) return 0; + a++; + } while (--alen); + return 1; +} + +/* {a,alen} := {a,alen} << shift. + Return the bits shifted out of the most significant digit of a. + Require 0 <= shift < BITS_PER_BNGDIGIT. */ +static bngdigit bng_generic_shift_left + (bng a/*[alen]*/, bngsize alen, + int shift) +{ + int shift2 = BNG_BITS_PER_DIGIT - shift; + bngdigit carry = 0; + if (shift > 0) { + for (/**/; alen > 0; alen--, a++) { + bngdigit d = *a; + *a = (d << shift) | carry; + carry = d >> shift2; + } + } + return carry; +} + +/* {a,alen} := {a,alen} >> shift. + Return the bits shifted out of the least significant digit of a. + Require 0 <= shift < BITS_PER_BNGDIGIT. */ +static bngdigit bng_generic_shift_right + (bng a/*[alen]*/, bngsize alen, + int shift) +{ + int shift2 = BNG_BITS_PER_DIGIT - shift; + bngdigit carry = 0; + if (shift > 0) { + for (a = a + alen - 1; alen > 0; alen--, a--) { + bngdigit d = *a; + *a = (d >> shift) | carry; + carry = d << shift2; + } + } + return carry; +} + +/* {a,alen} := {a,alen} + d * {b,blen}. Return carry out. + Require alen >= blen. */ +static bngdigit bng_generic_mult_add_digit + (bng a/*[alen]*/, bngsize alen, + bng b/*[blen]*/, bngsize blen, + bngdigit d) +{ + bngdigit out, ph, pl; + bngcarry carry; + + alen -= blen; + for (out = 0; blen > 0; blen--, a++, b++) { + bngdigit bd = *b; + /* ph:pl = double-digit product of b's current digit and d */ + BngMult(ph, pl, bd, d); + /* current digit of a += pl + out. Accumulate carries in ph. */ + BngAdd3(*a, ph, *a, pl, out); + /* prepare out for next iteration */ + out = ph; + } + if (alen == 0) return out; + /* current digit of a += out */ + BngAdd2(*a, carry, *a, out); + a++; + alen--; + /* Propagate carry */ + if (carry == 0 || alen == 0) return carry; + do { + if (++(*a) != 0) return 0; + a++; + } while (--alen); + return 1; +} + +/* {a,alen} := {a,alen} - d * {b,blen}. Return carry out. + Require alen >= blen. */ +static bngdigit bng_generic_mult_sub_digit + (bng a/*[alen]*/, bngsize alen, + bng b/*[blen]*/, bngsize blen, + bngdigit d) +{ + bngdigit out, ph, pl; + bngcarry carry; + + alen -= blen; + for (out = 0; blen > 0; blen--, a++, b++) { + bngdigit bd = *b; + /* ph:pl = double-digit product of b's current digit and d */ + BngMult(ph, pl, bd, d); + /* current digit of a -= pl + out. Accumulate carrys in ph. */ + BngSub3(*a, ph, *a, pl, out); + /* prepare out for next iteration */ + out = ph; + } + if (alen == 0) return out; + /* current digit of a -= out */ + BngSub2(*a, carry, *a, out); + a++; + alen--; + /* Propagate carry */ + if (carry == 0 || alen == 0) return carry; + do { + if ((*a)-- != 0) return 0; + a++; + } while (--alen); + return 1; +} + +/* {a,alen} := {a,alen} + {b,blen} * {c,clen}. Return carry out. + Require alen >= blen + clen. */ +static bngcarry bng_generic_mult_add + (bng a/*[alen]*/, bngsize alen, + bng b/*[blen]*/, bngsize blen, + bng c/*[clen]*/, bngsize clen) +{ + bngcarry carry; + for (carry = 0; clen > 0; clen--, c++, alen--, a++) + carry += bng_mult_add_digit(a, alen, b, blen, *c); + return carry; +} + +/* {a,alen} := 2 * {a,alen} + {b,blen}^2. Return carry out. + Require alen >= 2 * blen. */ +static bngcarry bng_generic_square_add + (bng a/*[alen]*/, bngsize alen, + bng b/*[blen]*/, bngsize blen) +{ + bngcarry carry1, carry2; + bngsize i, aofs; + bngdigit ph, pl, d; + + /* Double products */ + for (carry1 = 0, i = 1; i < blen; i++) { + aofs = 2 * i - 1; + carry1 += bng_mult_add_digit(a + aofs, alen - aofs, + b + i, blen - i, b[i - 1]); + } + /* Multiply by two */ + carry1 = (carry1 << 1) | bng_shift_left(a, alen, 1); + /* Add square of digits */ + carry2 = 0; + for (i = 0; i < blen; i++) { + d = b[i]; + BngMult(ph, pl, d, d); + BngAdd2Carry(*a, carry2, *a, pl, carry2); + a++; + BngAdd2Carry(*a, carry2, *a, ph, carry2); + a++; + } + alen -= 2 * blen; + if (alen > 0 && carry2 != 0) { + do { + if (++(*a) != 0) { carry2 = 0; break; } + a++; + } while (--alen); + } + return carry1 + carry2; +} + +/* {a,len-1} := {b,len} / d. Return {b,len} modulo d. + Require MSD of b < d. + If BngDivNeedsNormalization is defined, require d normalized. */ +static bngdigit bng_generic_div_rem_norm_digit + (bng a/*[len-1]*/, bng b/*[len]*/, bngsize len, bngdigit d) +{ + bngdigit topdigit, quo, rem; + intnat i; + + topdigit = b[len - 1]; + for (i = len - 2; i >= 0; i--) { + /* Divide topdigit:current digit of numerator by d */ + BngDiv(quo, rem, topdigit, b[i], d); + /* Quotient is current digit of result */ + a[i] = quo; + /* Iterate with topdigit = remainder */ + topdigit = rem; + } + return topdigit; +} + +#ifdef BngDivNeedsNormalization +/* {a,len-1} := {b,len} / d. Return {b,len} modulo d. + Require MSD of b < d. */ +static bngdigit bng_generic_div_rem_digit + (bng a/*[len-1]*/, bng b/*[len]*/, bngsize len, bngdigit d) +{ + bngdigit rem; + int shift; + + /* Normalize d and b */ + shift = bng_leading_zero_bits(d); + d <<= shift; + bng_shift_left(b, len, shift); + /* Do the division */ + rem = bng_div_rem_norm_digit(a, b, len, d); + /* Undo normalization on b and remainder */ + bng_shift_right(b, len, shift); + return rem >> shift; +} +#endif + +/* {n+dlen, nlen-dlen} := {n,nlen} / {d, dlen}. + {n, dlen} := {n,nlen} modulo {d, dlen}. + Require nlen > dlen and MSD of n < MSD of d. + (This implies MSD of d > 0). */ +static void bng_generic_div_rem + (bng n/*[nlen]*/, bngsize nlen, + bng d/*[dlen]*/, bngsize dlen) +{ + bngdigit topden, quo, rem; + int shift; + bngsize i, j; + + /* Normalize d */ + shift = bng_leading_zero_bits(d[dlen - 1]); + /* Note that no bits of n are lost by the following shift, + since n[nlen-1] < d[dlen-1] */ + bng_shift_left(n, nlen, shift); + bng_shift_left(d, dlen, shift); + /* Special case if d is just one digit */ + if (dlen == 1) { + *n = bng_div_rem_norm_digit(n + 1, n, nlen, *d); + } else { + topden = d[dlen - 1]; + /* Long division */ + for (j = nlen - 1; j >= dlen; j--) { + i = j - dlen; + /* At this point: + - the current numerator is n[j] : ...................... : n[0] + - to be subtracted quo times: d[dlen-1] : ... : d[0] : 0... : 0 + (there are i zeroes at the end) */ + /* Under-estimate the next digit of the quotient (quo) */ + if (topden + 1 == 0) + quo = n[j]; + else + BngDiv(quo, rem, n[j], n[j - 1], topden + 1); + /* Subtract d * quo (shifted i places) from numerator */ + n[j] -= bng_mult_sub_digit(n + i, dlen, d, dlen, quo); + /* Adjust if necessary */ + while (n[j] != 0 || bng_compare(n + i, dlen, d, dlen) >= 0) { + /* Numerator is still bigger than shifted divisor. + Increment quotient and subtract shifted divisor. */ + quo++; + n[j] -= bng_sub(n + i, dlen, d, dlen, 0); + } + /* Store quotient digit */ + n[j] = quo; + } + } + /* Undo normalization on remainder and divisor */ + bng_shift_right(n, dlen, shift); + bng_shift_right(d, dlen, shift); +} + +/**** Construction of the table of operations ****/ + +struct bng_operations bng_ops = { + bng_generic_add_carry, + bng_generic_add, + bng_generic_sub_carry, + bng_generic_sub, + bng_generic_shift_left, + bng_generic_shift_right, + bng_generic_mult_add_digit, + bng_generic_mult_sub_digit, + bng_generic_mult_add, + bng_generic_square_add, + bng_generic_div_rem_norm_digit, +#ifdef BngDivNeedsNormalization + bng_generic_div_rem_digit, +#else + bng_generic_div_rem_norm_digit, +#endif + bng_generic_div_rem +}; + +void bng_init(void) +{ +#ifdef BNG_SETUP_OPS + BNG_SETUP_OPS; +#endif +} |