bigint.hh

Go to the documentation of this file.

#ifndef CLICK_BIGINT_HH
#define CLICK_BIGINT_HH 1
/*
 * bigint.hh -- multiple-precision arithmetic
 * Eddie Kohler and the authors of GMP
 *
 * Copyright (c) 2008 Meraki, Inc.
 * Derived from the GNU Multiple Precision Arithmetic Library, which is
 * Copyright (c) 2001-2008 Free Software Foundation, Inc.
 *
 * This source code is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * Using the interfaces provided by this file does not create a derived work.
 */
#include <click/integers.hh>
#include <click/type_traits.hh>
#include <click/string.hh>
CLICK_DECLS
 
template <typename L, typename Lhalf = L>
class Bigint { public:
 
    typedef L limb_type;
 
    typedef Lhalf half_limb_type;
 
    enum {
        limb_bits = sizeof(limb_type) * 8,
        half_limb_bits = sizeof(limb_type) * 4
    };
    enum {
        limb_high_bit = limb_type(1) << (limb_bits - 1),
        limb_half = limb_type(1) << half_limb_bits,
        limb_low_mask = limb_half - 1,
        limb_zero = limb_type(0)
    };
 
    static inline half_limb_type low(limb_type x) {
        return x & limb_low_mask;
    }
 
    static inline half_limb_type high(limb_type x) {
        return x >> half_limb_bits;
    }
 
    static void add(limb_type &x1, limb_type &x0,
                    limb_type a1, limb_type a0, limb_type b1, limb_type b0) {
        x0 = a0 + b0;
        x1 = a1 + b1 + (x0 < a0);
    }
 
    static void multiply(limb_type &x1, limb_type &x0,
                         limb_type a, limb_type b) {
        int_multiply(a, b, x0, x1);
    }
 
    static void multiply_half(limb_type &x1, limb_type &x0,
                              limb_type a, half_limb_type b) {
        if (has_fast_int_multiply<limb_type>::value)
            multiply(x1, x0, a, b);
        else {
            half_limb_type al = low(a), ah = high(a);
 
            limb_type r0 = (limb_type) al * b;
 
            limb_type r1 = (limb_type) ah * b;
            limb_type r2 = high(r0) + r1;
 
            x1 = (r2 < r1 ? limb_half : limb_zero) + high(r2);
            x0 = (r2 << half_limb_bits) + low(r0);
        }
    }
 
    static limb_type inverse(limb_type x) {
        limb_type y1 = ~x, y0 = ~limb_zero;
 
        // now calculate [y1y0] / x
        limb_type x1 = high(x), x0 = low(x);
        limb_type q1 = y1 / x1;
        limb_type r1 = y1 - q1 * x1;
        limb_type m = q1 * x0;
        r1 = r1 * limb_half | high(y0);
        if (r1 < m) {
            --q1, r1 += x;
            if (r1 >= x && r1 < m)
                --q1, r1 += x;
        }
        r1 -= m;
 
        limb_type q0 = r1 / x1;
        limb_type r0 = r1 - q0 * x1;
        m = q0 * x0;
        r0 = r0 * limb_half | low(y0);
        if (r0 < m) {
            --q0, r0 += x;
            if (r0 >= x && r0 < m)
                --q0, r0 += x;
        }
 
        return q1 * limb_half | q0;
    }
 
    static click_uintmax_t set(limb_type *a, int n, click_uintmax_t b) {
        while (n > 0) {
            *a++ = b;
            n--;
            b >>= limb_bits;
        }
        return b;
    }
 
    static limb_type multiply_add(limb_type *x, const limb_type *a, int n,
                                  limb_type b) {
        limb_type carry = 0;
        do {
            limb_type x0, x1;
            multiply(x1, x0, *a++, b);
            x0 += carry;
            carry = (x0 < carry) + x1;
            x0 += *x;
            carry += (x0 < *x);
            *x++ = x0;
        } while (--n != 0);
        return carry;
    }
 
    static limb_type multiply(limb_type *x, const limb_type *a, int n,
                              limb_type b, limb_type carry = 0) {
        do {
            limb_type x0, x1;
            multiply(x1, x0, *a++, b);
            x0 += carry;
            carry = (x0 < carry) + x1;
            *x++ = x0;
        } while (--n != 0);
        return carry;
    }
 
    static limb_type multiply_half(limb_type *x, const limb_type *a, int n,
                                   half_limb_type b, limb_type carry = 0) {
        do {
            limb_type x0, x1;
            multiply_half(x1, x0, *a++, b);
            x0 += carry;
            carry = (x0 < carry) + x1;
            *x++ = x0;
        } while (--n != 0);
        return carry;
    }
 
    static limb_type divide(limb_type *x, const limb_type *a, int n,
                            limb_type b) {
        x += n - 1;   /* Make x point at most significant quotient limb */
        a += n - 1;   /* Make a point at most significant dividend limb */
        limb_type r = 0;
 
        if (b & limb_high_bit) {
            /* High quotient limb is 0 or 1, skip a divide step. */
            r = *a;
            *x = (r >= b);
            r -= (b & -*x);
            --x, --a, --n;
 
            /* Multiply-by-inverse, divisor already normalized. */
            limb_type b_inverse = inverse(b);
            while (n > 0) {
                preinverted_divide(*x, r, r, *a, b, b_inverse);
                --x, --a, --n;
            }
 
        } else {
            /* Most significant bit of divisor == 0.  */
            /* Skip a division if high < divisor (high quotient 0).  Testing
               here before normalizing will still skip as often as
               possible.  */
            if (*a < b) {
                r = *a;
                *x = 0;
                --x, --a, --n;
            }
 
            if (n != 0) {
                int norm = ffs_msb(b) - 1; // number of most significant 0 bits
                b <<= norm;
                r <<= norm;
                limb_type b_inverse = inverse(b);
 
                limb_type a1 = *a;
                r |= (a1 >> (limb_bits - norm));
                --a, --n;
                while (n > 0) {
                    limb_type a0 = *a;
                    preinverted_divide(*x, r, r, ((a1 << norm) | (a0 >> (limb_bits - norm))), b, b_inverse);
                    --x, --a, --n;
                    a1 = a0;
                }
                preinverted_divide(*x, r, r, a1 << norm, b, b_inverse);
                r >>= norm;
            }
        }
 
        return r;
    }
 
    static String unparse_clear(limb_type *x, int n, int base = 10, bool uppercase = false) {
        // need d chars, min d s.t. 10^d >= 2^(sizeof(limb_type) * 8 * n)
        // == min d s.t. d >= sizeof(limb_type) * 8 * n / lg 10
        int div = (base >= 16 ? 4 : (base >= 8 ? 3 : 1));
        String s = String::make_uninitialized((n * limb_bits) / div + 1);
        char *q = const_cast<char *>(s.end());
        assert(base >= 2 && base <= 36);
        while (1) {
            while (n > 0 && x[n - 1] == 0)
                --n;
            if (n == 0)
                break;
            int r = divide(x, x, n, base);
            if (r <= 9)
                *--q = '0' + r;
            else if (uppercase)
                *--q = 'A' + r - 10;
            else
                *--q = 'a' + r - 10;
        }
        assert(q >= s.begin());
        if (q == s.end())
            *--q = '0';
        return s.substring(q, s.end());
    }
 
  private:
 
    /* Divide the two-limb number in (a1,a0) by b, with b_inverse being the
       largest limb not larger than (2**(2*limb_bits))/b - (2**limb_bits).  If
       this would yield overflow, b_inverse should be the largest possible
       number (i.e., only ones).  For correct operation, the most significant
       bit of b has to be set.  Put the quotient in q and the remainder in
       r.  */
    /* Like udiv_qrnnd_preinv, but branch-free. */
    static void preinverted_divide(limb_type &q, limb_type &r,
                                   limb_type a1, limb_type a0, limb_type b,
                                   limb_type b_inverse) {
        limb_type a0_mask = (a0 & limb_high_bit ? ~limb_zero : limb_zero);
        limb_type a0_adjusted = a0 + (b & a0_mask);
        limb_type x1, x0;
        multiply(x1, x0, b_inverse, a1 - a0_mask);
        add(x1, x0, x1, x0, a1, a0_adjusted);
        limb_type q1 = ~x1;
        multiply(x1, x0, q1, b);
        add(x1, x0, x1, x0, a1, a0);
        x1 -= b;
        r = x0 + (b & x1);
        q = x1 - q1;
    }
 
};
 
typedef Bigint<uint64_t> bigint;
 
CLICK_ENDDECLS
#endif