/***
 * This software was written by Nils Maier. No copyright is claimed, and the
 * software is hereby placed in the public domain.
 *
 * In case this attempt to disclaim copyright and place the software in the
 * public domain is deemed null and void in your jurisdiction, then the
 * software is Copyright 2004,2013 Nils Maier and it is hereby released to the
 * general public under the following terms:
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted.
 * There's ABSOLUTELY NO WARRANTY, express or implied.
 */

#ifndef BIGNUM_H
#define BIGNUM_H

#include <cstring>
#include <algorithm>
#include <memory>
#include <stdint.h>

namespace bignum {

  template<size_t dim>
  class ulong {

  public:
    typedef char char_t;
    typedef std::make_unsigned<char_t>::type uchar_t;

  private:
    std::unique_ptr<char_t[]> buf_;

  public:
    inline ulong()
      : buf_(new char_t[dim]())
    {
    }
    inline ulong(size_t t)
      : buf_(new char_t[dim]())
    {
      memcpy(buf_.get(), (char_t*)&t, sizeof(t));
    }
    inline ulong(const ulong<dim>& rhs)
      : buf_(new char_t[dim]())
    {
      memcpy(buf_.get(), rhs.buf_.get(), dim);
    }
    explicit inline ulong(const char_t *data, size_t size)
      : buf_(new char_t[dim]())
    {
      if (size > dim) {
        throw std::bad_alloc();
      }
      memcpy(buf_.get(), data, size);
    }

    virtual ~ulong() {}

    ulong<dim>& operator=(const ulong<dim>& rhs) {
      memcpy(buf_.get(), rhs.buf_.get(), dim);
      return *this;
    }

    bool operator==(const ulong<dim>& rhs) const {
      return memcmp(buf_.get(), rhs.buf_.get(), dim) == 0;
    }
    bool operator!=(const ulong<dim>& rhs) const {
      return memcmp(buf_.get(), rhs.buf_.get(), dim) != 0;
    }
    bool operator>(const ulong<dim>& rhs) const {
      const auto b1 = buf_.get();
      const auto b2 = rhs.buf_.get();
      for (ssize_t i = dim - 1; i >= 0; --i) {
        for (ssize_t j = 1; j >= 0; --j) {
          char_t t = ((uchar_t)(b1[i] << 4 * (1 - j))) >> 4;
          char_t r = ((uchar_t)(b2[i] << 4 * (1 - j))) >> 4;
          if (t != r) {
            return t > r;
          }
        }
      }
      return false;
    }
    bool operator>=(const ulong<dim>& rhs) const {
      return *this == rhs || *this > rhs;
    }
    bool operator<(const ulong<dim>& rhs) const {
      return !(*this >= rhs);
    }
    bool operator<=(const ulong<dim>& rhs) const {
      return *this == rhs || *this < rhs;
    }

    ulong<dim> operator+(const ulong<dim>& rhs) const {
      ulong<dim> rv;
      const auto b1 = buf_.get();
      const auto b2 = rhs.buf_.get();
      const auto rb = rv.buf_.get();
      bool base = false;
      for (size_t i = 0; i < dim; ++i) {
        for (ssize_t j = 0; j < 2; ++j) {
          char_t t = ((uchar_t)(b1[i] << 4 * (1 - j))) >> 4;
          char_t r = ((uchar_t)(b2[i] << 4 * (1 - j))) >> 4;
          if (base) {
            t++;
          }
          if (r + t >= 16) {
            rb[i] += (t + r - 16) << j * 4;
            base = true;
          }
          else {
            rb[i] += (t + r) << j * 4;
            base = false;
          }
        }
      }
      return rv;
    }
    ulong<dim>& operator+=(const ulong<dim>& rhs) {
      *this = *this + rhs;
      return *this;
    }
    ulong<dim>& operator++() {
      *this = *this + 1;
      return *this;
    }
    ulong<dim> operator++(int) {
      ulong<dim> tmp = *this;
      *this = *this + 1;
      return tmp;
    }

    ulong<dim> operator-(const ulong<dim>& rhs) const {
      ulong<dim> rv;
      const auto b1 = buf_.get();
      const auto b2 = rhs.buf_.get();
      const auto rb = rv.buf_.get();
      bool base = false;
      for (size_t i = 0; i < dim; ++i) {
        for (ssize_t j = 0; j < 2; ++j) {
          char_t t = ((uchar_t)(b1[i] << 4 * (1 - j))) >> 4;
          char_t r = ((uchar_t)(b2[i] << 4 * (1 - j))) >> 4;
          if (base) {
            t--;
          }
          if (t >= r) {
            rb[i] += (t - r) << j * 4;
            base = false;
          }
          else {
            rb[i] += (t + 16 - r) << j * 4;
            base = true;
          }
        }
      }
      return rv;
    }
    ulong<dim>& operator-=(const ulong<dim>& rhs) {
      *this = *this - rhs;
      return *this;
    }
    ulong<dim>& operator--() {
      *this = *this - 1;
      return *this;
    }
    ulong<dim> operator--(int) {
      ulong<dim> tmp = *this;
      *this = *this - 1;
      return tmp;
    }

    ulong<dim> operator*(const ulong<dim>& rhs) const {
      ulong<dim> c = rhs, rv;
      const ulong<dim> null;
      size_t cap = c.capacity();
      while (c != null) {
        ulong<dim> tmp = *this;
        tmp.mul(cap - 1);
        rv += tmp;

        ulong<dim> diff(1);
        diff.mul(cap - 1);
        c -= diff;

        cap = c.capacity();
      }
      return rv;
    }
    ulong<dim>& operator*=(const ulong<dim>& rhs) {
      *this = *this * rhs;
      return *this;
    }

    ulong<dim> operator/(const ulong<dim>& rhs) const {
      ulong<dim> quotient, remainder;
      div(rhs, quotient, remainder);
      return quotient;
    }
    ulong<dim>& operator/=(const ulong<dim>& rhs) {
      *this = *this / rhs;
      return *this;
    }

    ulong<dim> operator%(const ulong<dim>& rhs) const {
      ulong<dim> quotient, remainder;
      div(rhs, quotient, remainder);
      return remainder;
    }
    ulong<dim>& operator%=(const ulong<dim>& rhs) {
      *this = *this % rhs;
      return *this;
    }

    ulong<dim> mul_mod(const ulong<dim>& mul, const ulong<dim>& mod) const {
      if (capacity() + mul.capacity() <= dim) {
        return (*this * mul) % mod;
      }
      ulong<dim*2> et(buf_.get(), dim),
        emul(mul.buf_.get(), dim),
        emod(mod.buf_.get(), dim),
        erv = (et * emul) % emod;
      ulong<dim> rv;
      erv.binary(rv.buf_.get(), dim);
      return rv;
    }

    std::unique_ptr<char_t[]> binary() const {
      ulong<dim> c = *this;
      std::unique_ptr<char_t[]> rv;
      rv.swap(c.buf_);
      return rv;
    }
    void binary(char_t *buf, size_t len) const {
      memcpy(buf, buf_.get(), std::min(dim, len));
    }

    size_t length() const {
      return dim;
    }

  private:
    size_t capacity() const {
      size_t rv = dim * 2;
      const auto b = buf_.get();
      for (ssize_t i = dim - 1; i >= 0; --i) {
        char_t f = b[i] >> 4;
        char_t s = (b[i] << 4) >> 4;
        if (!f && !s) {
          rv -= 2;
          continue;
        }
        if (!f) {
          --rv;
        }
        return rv;
      }
      return rv;
    }

    void mul(size_t digits) {
      ulong<dim> tmp = *this;
      auto bt = tmp.buf_.get();
      auto b = buf_.get();
      memset(b, 0, dim);
      const size_t npar = digits % 2;
      const size_t d2 = digits / 2;
      for (size_t i = d2; i < dim; ++i) {
        for (size_t j = 0; j < 2; ++j) {
          char_t c = ((uchar_t)(bt[(dim - 1) - i] << 4 * (1 - j))) >> 4;
          char_t r = c << (npar * (1 - j) * 4 + (1 - npar) * j * 4);
          ssize_t idx = i - d2 - npar * j;
          if (idx >= 0) {
            b[(dim - 1) - idx] += r;
          }
        }
      }
    }

    void div(const ulong<dim>& d, ulong<dim>& q, ulong<dim>& r) const {
      ulong<dim> tmp = d;
      r = *this;
      q = 0;
      size_t cr = r.capacity();
      const size_t cd = d.capacity();
      while (cr > cd) {
        tmp = d;
        tmp.mul(cr - cd - 1);
        ulong<dim> qt(1);
        qt.mul(cr - cd - 1);
        ulong<dim> t = tmp;
        t.mul(1);
        if (r >= t) {
          tmp = t;
          qt.mul(1);
        }
        while (r >= tmp) {
          r -= tmp;
          q += qt;
        }
        cr = r.capacity();
      }
      while (r >= d) {
        r -= d;
        ++q;
      }
    }
  };

} // namespace bignum

#endif
