libecc/polynomial.h

Go to the documentation of this file.

//
//
// This file is part of the libecc package.
// Copyright (C) 2002 - 2004 by
//
// Carlo Wood, Run on IRC <carlo@alinoe.com>
// RSA-1024 0x624ACAD5 1997-01-26                    Sign & Encrypt
// Fingerprint16 = 32 EC A7 B6 AC DB 65 A6  F6 F6 55 DD 1C DC FF 61
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version.
//
// 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
//

#ifndef LIBECC_POLYNOMIAL_H
#define LIBECC_POLYNOMIAL_H

#include <stdexcept>
#include <libecc/bitset.h>
#include <libecc/debug.h>
#if ECC_DEBUGOUTPUT
#include <libcwd/cwprint.h>
#endif

#if ECC_DEBUG
#define LIBECC_AUGMENTED 1
#define LIBECC_INPLACE (1 || !LIBECC_AUGMENTED)
#define LIBECC_SWAPCOLUMNS (1 || LIBECC_INPLACE)
#else
// Don't change these.
#define LIBECC_AUGMENTED 0
#define LIBECC_INPLACE 1
#define LIBECC_SWAPCOLUMNS 1
#endif

namespacelibecc {

// Forward declarations.
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  classpolynomial;
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  polynomial<m, k, k1, k2> operator*(polynomial<m, k, k1, k2> const&, polynomial<m, k, k1, k2> const&);
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  polynomial<m, k, k1, k2> operator/(polynomial<m, k, k1, k2> const&, polynomial<m, k, k1, k2> const&);
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  bool operator==(polynomial<m, k, k1, k2> const&, polynomial<m, k, k1, k2> const&);
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  bool operator!=(polynomial<m, k, k1, k2> const&, polynomial<m, k, k1, k2> const&);
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  std::ostream& operator<<(std::ostream&, polynomial<m, k, k1, k2> const&);
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  std::ostream& operator<<(std::ostream&, typename polynomial<m, k, k1, k2>::xor_type const&);
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  typename polynomial<m, k, k1, k2>::xor_type operator+(polynomial<m, k, k1, k2> const&, polynomial<m, k, k1, k2> const&);
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  typename polynomial<m, k, k1, k2>::xor_type operator-(polynomial<m, k, k1, k2> const&, polynomial<m, k, k1, k2> const&);

template<unsigned int m, unsigned int k, unsigned int k1 = 0, unsigned int k2 = 0>
  classpolynomial {
    public:
      typedef Operator::bitsetExpression<m, false, false, Operator::bitsetXOR> xor_type;

      // Fix this if you add members in front of M_coefficients.
      static size_t const offsetof_vector = bitset<m>::offsetof_vector;

    private:
      bitset<m> M_coefficients;
      static polynomial<m, k, k1, k2> const one;
      static bool S_normal_initialized;
      static bitset<m> S_normal;

    public:
      static polynomial const& unity(void) { return one; }

    public:
      polynomial(void) { }

      explicit polynomial(bitset_digit_t coefficients) : M_coefficients(coefficients) { }

      polynomial(polynomial const& p) : M_coefficients(p.M_coefficients) { }

      explicit polynomial(bitset<m> const& coefficients) : M_coefficients(coefficients) { }

      polynomial(std::string const& coefficients) : M_coefficients(coefficients) { }

      polynomial(xor_type const& expression) : M_coefficients(expression) { }

      polynomial& operator=(polynomial const& p) { M_coefficients = p.M_coefficients; return *this; }

      polynomial& operator=(bitset<m> const& coefficients) { M_coefficients = coefficients; return *this; }

      polynomial& operator=(xor_type const& expression);

      polynomial(polynomial const& b, polynomial const& c);

      static unsigned int const square_digits = 2 * bitset_base<m>::digits + 4;

      polynomial& square(bitset_digit_t* tmpbuf) const; // tmpbuf must be an array of `square_digits' bitset_digit_t.

      bool sqrt(void);

      // The field arithmetic is implemented in terms of operations on the bits.
      polynomial& operator+=(polynomial const& p) { M_coefficients ^= p.M_coefficients; return *this; }

      polynomial& operator-=(polynomial const& p) { M_coefficients ^= p.M_coefficients; return *this; }

      polynomial& operator*=(polynomial const& p);
#ifdef LIBECC_DOXYGEN
      // Stupid doxygen.
      polynomial& operator*=(typename polynomial<m, k, k1, k2>::xor_type const& expr);
#else
      // The real prototype.
      polynomial& operator*=(xor_type const& expr);
#endif

      polynomial& operator/=(polynomial const& p);
#ifdef LIBECC_DOXYGEN
      // Stupid doxygen.
      polynomial& operator/=(typename polynomial<m, k, k1, k2>::xor_type const& expr);
#else
      // The real prototype.
      polynomial& operator/=(xor_type const& expr);
#endif

      static bitset<m> const& normal(void) { if (!S_normal_initialized) calculate_normal(); return S_normal; }

      int trace(void) const
     {
        // This method was invented by me, so give me credit for it when you use it somewhere. Thank you.
        // Carlo Wood <carlo@alinoe.com> -- 4 December 2004.
        int tr = 0;
        if ((m & 1))
          tr = M_coefficients.template test<0>();
        if (((m - k) & 1))
          tr ^= M_coefficients.template test<m - k>();
        if (k1)
        {
          if (((m - k1) & 1))
            tr ^= M_coefficients.template test<m - k1>();
          if (((m - k2) & 1))
            tr ^= M_coefficients.template test<m - k2>();
        }
        return tr;
      }

      friend xor_type operator+ <>(polynomial const& p1, polynomial const& p2);

      friend xor_type operator- <>(polynomial const& p1, polynomial const& p2);

      friend polynomial operator* <>(polynomial const& p1, polynomial const& p2);
#ifdef LIBECC_DOXYGEN
      // Only added for documentational reasons.
      friend bool operator*(polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2);
      friend bool operator*(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2>::xor_type const& expr);
#endif

      friend polynomial operator/ <>(polynomial const& p1, polynomial const& p2);
#ifdef LIBECC_DOXYGEN
      // Only added for documentational reasons.
      friend bool operator/(polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2);
      friend bool operator/(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2>::xor_type const& expr);
#endif

      friend bool operator== <>(polynomial const& p1, polynomial const& p2);
#ifdef LIBECC_DOXYGEN
      // Only added for documentational reasons.
      friend bool operator==(polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2);
      friend bool operator==(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2>::xor_type const& expr);
#endif

      friend bool operator!= <>(polynomial const& p1, polynomial const& p2);
#ifdef LIBECC_DOXYGEN
      // Only added for documentational reasons.
      friend bool operator!=(polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2);
      friend bool operator!=(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2>::xor_type const& expr);
#endif

      friend std::ostream& operator<< <>(std::ostream& os, polynomial const& p);
#ifdef LIBECC_DOXYGEN
      // Only added for documentational reasons.
      friend std::ostream& operator<<(std::ostream& os, polynomial<m, k, k1, k2>::xor_type const& expr);
#endif

      bitset<m> const& get_bitset(void) const{ return M_coefficients; }

      bitset<m>& get_bitset(void) { return M_coefficients; }

    private:
      static void reduce(bitset_digit_t* buf);
      static bitset_digit_t reducea(bitset_digit_t* a);
      static void calculate_normal(void);

      void multiply_with(polynomial const& p1, bitset<m>& result) const;
#if ECC_DEBUG
#if LIBECC_AUGMENTED
      void print_matrix(bitset<2 * m> const* matrix, bitset<m> const& pivotted);
#else
      void print_matrix(bitset<m> const* matrix, bitset<m> const& pivotted);
#endif
#endif
  };

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  polynomial<m, k, k1, k2> const polynomial<m, k, k1, k2>::one(1);

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  bool polynomial<m, k, k1, k2>::sqrt(void)
  {
    if (!k1)
    {
      bitset<m> highbits;
      highbits.reset();

      // First convert all odd powers into even powers
      if ((m & 1) == 1)
      {
        if ((k & 1) == 1)               // m and k are odd?
        {
          for(unsigned int bit = 1; bit < m; bit += 2)
          {
            if (M_coefficients.test(bit))
            {
              if (bit >= m - k)
                highbits.flip(bit + k - m);
              else
                M_coefficients.flip(bit + k);
              highbits.flip(bit);
            }
          }
        }
        else                    // m is odd and k is even
        {
          for(unsigned int bit = 1; bit < m; bit += 2)
          {
            if (M_coefficients.test(bit))
            {
              if (bit >= m - k)
              {
                M_coefficients.flip(bit + 2 * k - m);
                M_coefficients.flip(bit + k - m);
              }
              else
                M_coefficients.flip(bit + k);
              highbits.flip(bit);
            }
          }
        }
      }
      else if ((k & 1) == 1)    // m is even and k is odd
      {
        for(unsigned int bit = 1; bit < m; bit += 2)
        {
          if (M_coefficients.test(bit))
          {
            if (bit < k)
            {
              M_coefficients.flip(bit + k);
              M_coefficients.flip(bit + m - k);
              highbits.flip(bit + m - k);
            }
            else
            {
              M_coefficients.flip(bit - k);
              highbits.flip(bit - k);
            }
          }
        }
      }
      else                      // m and k are both even (actually, this should never be used as reduction polynomial).
      {
        for(unsigned int bit = 1; bit < m; bit += 2)
          if (M_coefficients.test(bit))
            return false;               // This can't be a square
      }

      // Next handle the remaining even powers
      unsigned int bit_to = 1;
      for(unsigned int bit = 2; bit < m; bit += 2)
      {
        if (M_coefficients.test(bit))
          M_coefficients.set(bit_to);
        else
          M_coefficients.clear(bit_to);
        ++bit_to;
      }
      for(unsigned int bit = m % 2; bit < m; bit += 2)
      {
        if (highbits.test(bit))
          M_coefficients.set(bit_to);
        else
          M_coefficients.clear(bit_to);
        ++bit_to;
      }
    }
    else
    {
      structRoot {
        polynomial<m, k, k1, k2> value;
        Root(polynomial<m, k, k1, k2> const& p) : value(p)
        {
          bitset_digit_t p2buf[libecc::polynomial<m, k, k1, k2>::square_digits];
          polynomial<m, k, k1, k2>& p2 = value.square(p2buf);
          bitset_digit_t p4buf[libecc::polynomial<m, k, k1, k2>::square_digits];
          polynomial<m, k, k1, k2>& p4 = p2.square(p4buf);
          for (unsigned int i = 1; i < m / 2; ++i)
          {
            p4.square(p2buf);
            p2.square(p4buf);
          }
          value = (m % 2 == 0) ? p2 : p4;
        }
      };
      static Root const root_of_t(polynomial<m, k, k1, k2>(2));
      polynomial<m, k, k1, k2> tmp(0);
      bitset<m> tmp2;
      tmp2.reset();
      for(unsigned int bit = 0; bit < m / 2; ++bit)
      {
        if (M_coefficients.test(2 * bit))
          tmp2.set(bit);
        if (M_coefficients.test(2 * bit + 1))
          tmp.get_bitset().set(bit);
      }
      if (m % 2 == 1 && M_coefficients.test(m - 1))
        tmp2.set(m / 2);
      M_coefficients = tmp2;
      *this += tmp * root_of_t.value;
    }
    return true;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>&
  polynomial<m, k, k1, k2>::operator*=(polynomial const& p)
  {
    multiply_with(p, M_coefficients);
    return *this;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>&
  polynomial<m, k, k1, k2>::operator*=(typename polynomial<m, k, k1, k2>::xor_type const& expr)
  {
    return (*this *= polynomial<m, k, k1, k2>(expr));
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>&
  polynomial<m, k, k1, k2>::operator=(xor_type const& expression)
  {
    M_coefficients = expression;
    return *this;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  void
  polynomial<m, k, k1, k2>::multiply_with(polynomial const& p1, bitset<m>& result) const
 {
    bitset_digit_t output[bitset<m>::digits * 2] __attribute__ ((aligned (8)));

    // Find the first non-zero digit in the input polynomial of this object.
    unsigned int digit = 0;
    while(M_coefficients.digit(digit) == 0)             // Still zero?
    {
      output[digit] = 0;                                // That means that the output will end on zero too.
      if (++digit == bitset<m>::digits)
      {
        result.reset();                                 // The whole polynomial is zero, the result will be zero too.
        return;
      }
    }
    unsigned int uninitialized_digit = digit;           // The next digit of `output' that has not yet been initialized.
    // Find the first digit in the input polynomial of this object whose first bit is set.
    for(; digit < bitset<m>::digits; ++digit)
    {
      if ((M_coefficients.digit(digit) & 1))            // Is the first bit set?
      {
        // Set the output to p1 times this bit.
        for (unsigned int d = 0; d < bitset<m>::digits; ++d)
          output[d + digit] = p1.get_bitset().digit(d);
        uninitialized_digit = bitset<m>::digits + digit;
        ++digit;                                        // Set to the next input digit.
        break;
      }
      output[digit] = 0;                                // Initialize this digit of the output to 0.
      ++uninitialized_digit;
    }
    // Set the remaining digits to zero, if any.
    for(unsigned int remaining_digit = uninitialized_digit; remaining_digit < sizeof(output) / sizeof(bitset_digit_t); ++remaining_digit)
      output[remaining_digit] = 0;
    // Find for the remaining input digits the ones that have their first bit set.
    for(; digit < bitset<m>::digits; ++digit)
      if ((M_coefficients.digit(digit) & 1))            // Is the first bit set?
      {
        // Add p1 times this bit to the output.
        for (unsigned int d = 0; d < bitset<m>::digits; ++d)
          output[d + digit] ^= p1.get_bitset().digit(d);
      }
    // Create a bitset that will contain p1, shifted at most bitset_digit_bits - 1 to the left.
    bitset<m + bitset_digit_bits - 1> shifted_p1;
    // Start with having it shifted 1 bit to the left.
    bitset_digit_t carry = 0;
    unsigned int d = 0;
    for(bitset_digit_t const* ptr = p1.get_bitset().digits_ptr(); ptr < p1.get_bitset().digits_ptr() + bitset<m>::digits; ++ptr, ++d)
    {
      shifted_p1.rawdigit(d) = (*ptr << 1) | carry;
      carry = *ptr >> (8 * sizeof(bitset_digit_t) - 1);
    }
    if (d < bitset<m + bitset_digit_bits - 1>::digits)
      shifted_p1.rawdigit(d) = carry;
    for(bitset_digit_t bitmask = 2;;)
    {
      for(unsigned int digit = 0; digit < bitset<m>::digits; ++digit)
        if ((M_coefficients.digit(digit) & bitmask))
        {
          for (unsigned int d = 0; d < shifted_p1.digits; ++d)
            output[d + digit] ^= shifted_p1.digit(d);
        }
      bitmask <<= 1;            // Next bit.
      if (bitmask == 0)         // Done?
        break;
      // Shift p1 one bit further to the left.
      shifted_p1.template shift_op<1, left, assign>(shifted_p1);
    }
    // Reduce the resulting output of the multiplication.
    reduce(output);
    // Copy the reduced output to `result'.
    std::memcpy(result.digits_ptr(), output, bitset<m>::digits * sizeof(bitset_digit_t));
  }

#if ECC_DEBUG
template<unsigned int m>
structdiv_tct {
  bitset_digit_t const* M_p;
  int M_deg;
  int M_low;
  div_tct(bitset_base<m> const& b, int deg, int low) : M_p(b.digits_ptr()), M_deg(deg), M_low(low) { }
  void print_on(std::ostream& os) const
 {
    int lowbit = (M_low >> bitset_digit_bits_log2) * bitset_digit_bits;
    if (lowbit > 0)
      lowbit = 0;
    for (int b = 2 * m - 1; b >= lowbit; --b)
    {
      if (b == M_deg)
        os << "\e[31m";
      int digitoffset = (b >> bitset_digit_bits_log2);
      bitset_digit_t mask = static_cast<bitset_digit_t>(1) << (b & (bitset_digit_bits - 1));
      if (M_p[digitoffset] & mask)
        os << '1';
      else
        os << '0';
      if (b == M_low)
        os << "\e[0m";
      if (b == 0)
        os << '.';
    }
  }
};
#endif

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  polynomial<m, k, k1, k2>&
  polynomial<m, k, k1, k2>::operator/=(polynomial const& p)
  {
    //std::cout << "Entering polynomial<" << m << ", " << k << ", " << k1 << ", " << k2 << ">::operator/=()" << std::endl;
#if ECC_DEBUG
    LibEccDout(dc::polynomial|noprefix_cf, "");
    //LibEccDout(dc::polynomial, "Entering polynomial<" << m << ", " << k << ", " << k1 << ", " << k2 << ">::operator/=()");
    polynomial<m, k, k1, k2> x(p.get_bitset());
    polynomial<m, k, k1, k2> y(M_coefficients);
    LibEccDout(dc::polynomial, "x(t) = " << x);
    LibEccDout(dc::polynomial|flush_cf, "y(t) = " << y);
#endif

    // The following algorithm is based on the algorithm described on
    // page 7 of http://research.sun.com/techrep/2001/smli_tr-2001-95.ps
    // with significant optimization changes by Carlo Wood.
    // The basis of the algorithm is to keep an invariants valid:
    // A * y = U * x  and  B * y = V * x
    // while reducing the size of A and B steadily to 0, until either is 1.
    // The size of A and B is defined as the distance between the lowest
    // and highest 1 in the bitset.

    // Make sure that there is enough space for a full bitset object
    // and align the bitsets on a multiple of bitset_digit_t.
    static unsigned int const digit_offset_UV = ((sizeof(bitset<m>) * 8 - 1) / bitset_digit_bits + 1);
    static unsigned int const offset_UV = digit_offset_UV * bitset_digit_bits;
    // Make room for exponents from at least t^-m till t^2m.
    static unsigned int const digit_size_UV = 3 * digit_offset_UV;
    // Variables A and B do not need this much space.
    static unsigned int const digit_size_AB = bitset<m>::digits;
    // One digit of padding, needed for xor_with_zero_padded.
    static unsigned int const padding_digit_size = 1;
    // 'F' (Floating-point polynomial) will be shifted to the right and
    // is therefore defined to run from t^-2m till t^2m.  This means it will
    // be shifted OVER the other bitsets, but we don't need those anymore anyway.
    static unsigned int const offset_F = 2 * offset_UV;
    static unsigned int const size_F = 2 * m + offset_F;

    // Declare stack space for four variables.
    // This array is going to be used as follows
    // (where padding = padding_digit_size; AB = digit_size_AB and do_UV = digit_offset_UV):
    //
    //          A:->         B:->                U:---->                       V:---->
    // [padding][AB][padding][AB][padding][do_UV][do_UV][do_UV][padding][do_UV][do_UV][do_UV][padding]
    //                                    <-- digit_size_UV -->         <-- digit_size_UV -->
    //                             [do_UV]                       [do_UV]
    //                             F:------------------------>   F:------------------------>
    //                                           R:---->  *or*                 R:---->
    //                             <- offset_F ->                <- offset_F ->
    //
    // This positioning is done to make the total number of bits that we have to clear in
    // as small as possible. R is just an alias for either U or V depending on how the
    // algorithm finishes. Note F overlaps with R and they alias eachother, F might also
    // overlap with B but that variable isn't used anymore then.
    //
    // In order to avoid aliasing problems (see http://dbp-consulting.com/StrictAliasing.pdf),
    // we put the bitsets in a struct together with bitset_digit_t for the padding and
    // use a union to give F a place in all of this. Because objects with constructors
    // aren't allowed in union's we have to use bitset_base.
    structUV_bitpool_t {
      bitset_digit_t offset1[digit_offset_UV];
      union{
        bitset_digit_t offset2[digit_size_UV - digit_offset_UV];
        bitset_base<m> UV;
      };
    };
    structbitpool_t {
      bitset_digit_t padding1[padding_digit_size];
      bitset_base<m> A;
      bitset_digit_t padding2[padding_digit_size];
      bitset_base<m> B;
      bitset_digit_t padding3[padding_digit_size];
      UV_bitpool_t U_bitpool;
      bitset_digit_t padding4[padding_digit_size];
      UV_bitpool_t V_bitpool;
      bitset_digit_t padding5[padding_digit_size];
    };
    structFU_bitpool_t {
      bitset_digit_t padding[3 * padding_digit_size + 2 * digit_size_AB - digit_offset_UV];
      bitset_base<size_F> FU;
    };
    structFV_bitpool_t {
      bitset_digit_t padding[4 * padding_digit_size + 2 * digit_size_AB + 2 * digit_offset_UV];
      bitset_base<size_F> FV;
    };
    unionstrict_aliasing_OK_t {
      bitpool_t ABUV;
      FU_bitpool_t FU;
      FV_bitpool_t FV;
    };
    strict_aliasing_OK_t bitpool __attribute__ ((__aligned__ (ECC_BITS)));
    std::memset((char*)&bitpool, 0, sizeof(bitpool));
    bitset_base<m>& A(bitpool.ABUV.A);
    bitset_base<m>& B(bitpool.ABUV.B);
    bitset_base<m>& U(bitpool.ABUV.U_bitpool.UV);
    bitset_base<m>& V(bitpool.ABUV.V_bitpool.UV);

#if ECC_DEBUG
    assert(sizeof(strict_aliasing_OK_t) == sizeof(bitpool_t));
    assert(sizeof(bitset_base<m>) == digit_size_AB * sizeof(bitset_digit_t));
    assert(sizeof(UV_bitpool_t) == digit_size_UV * sizeof(bitset_digit_t));
    assert(sizeof(bitpool_t) == (5 * padding_digit_size + 2 * digit_size_AB + 2 * digit_size_UV) * sizeof(bitset_digit_t));
    assert((bitset_digit_t*)&bitpool.FU.FU + 2 * digit_offset_UV == (bitset_digit_t*)&U);
    assert((bitset_digit_t*)&bitpool.FV.FV + 2 * digit_offset_UV == (bitset_digit_t*)&V);
#endif

    // The representation of U and V will be done with bitsets of size `digit_size_UV * bitset_digit_bits'.
    // This means that they contain powers of t with a negative exponent.
    // That is not a problem as those are well defined: t^(-n) = 1 / t^n.

    // Let rp = M(t) = t^m + t^k [+ t^k1 + t^k2] + 1.
#if ECC_DEBUG
    bitset<m + 1> rp("1");
    rp.template set<m>();
    rp.template set<k>();
    if (k1)
    {
      rp.template set<k1>();
      rp.template set<k2>();
    }
#endif

    // Let U(t) = y(t) (= M_coefficients).
    LibEccDout(dc::polynomial|flush_cf, "U <- y");
    U = M_coefficients;

    // Guess the maximum and minimum powers to be the possible limits.
    int degU = m - 1;
    int lowU = 0;

    // Let A(t) = x(t).
    LibEccDout(dc::polynomial|flush_cf, "A <- x");
    A = p.get_bitset();

    // Then
    //
    // A(t) * y(t) = U(t) * x(t)  [mod M(t)].

    // Let V(t) = 0
    // Let B = M(t)
    //
    // Then
    //
    // B(t) * y(t) = V(t) * x(t)  [mod M(t)].
    //
    // Let degA be the highest power of t in A.
    typename bitset<m>::const_reverse_iterator degA = A.rbegin();
    degA.find1();
    LibEccDout(dc::polynomial|flush_cf, "deg(A) == " << degA);

    // Let lowA be the lowest power of t in A.
    typename bitset<m>::const_iterator lowA = A.begin();
    lowA.find1();
    LibEccDout(dc::polynomial|flush_cf, "low(A) == " << lowA);

    unsigned int sizeA = degA.get_index() - lowA.get_index();

    // Let n = m - deg(A).
    unsigned int n = m - degA.get_index();
    //
    // Then B'(t) = B(t) - A(t) * t^n will have a degree less than m.
    // And
    //
    // B'(t) * y(t) = B(t) * y(t) - A(t) * y(t) * t^n =
    //              = V(t) * x(t) - U(t) * x(t) * t^n =
    //              = (V(t) - U(t) * t^n) * x(t) =
    //              = V'(t) * x(t)                      [mod M(t)].
    //
    // B <- B'
    LibEccDout(dc::polynomial|flush_cf, "B <- A * t^" << n << " + " << cwprint_using(rp, &bitset<m+1>::base2_print_on));
    B.xor_with_zero_padded(A, lowA.get_index(), degA.get_index(), n);
    B.template flip<m>();
    B.template flip<k>();
    if (k1)
    {
      B.template flip<k1>();
      B.template flip<k2>();
    }
    B.template flip<0>();

    // Let degB be the highest power of t in B.
    typename bitset<m>::const_reverse_iterator degB = B.rbegin();
    degB.find1();
    LibEccDout(dc::polynomial|flush_cf, "deg(B) == " << degB);

    // Let lowB be the lowest power of t in B.
    typename bitset<m>::const_iterator lowB = B.begin();
    lowB.find1();
    LibEccDout(dc::polynomial|flush_cf, "low(B) == " << lowB);

    // V <- V'
    LibEccDout(dc::polynomial|flush_cf, "V <- U * t^" << n <<
        "  [mod " << cwprint_using(rp, &bitset<m + 1>::base2_print_on) << "]");
    V.xor_with_zero_padded(U, 0, m - 1, n);

    int degV = degU + n;
    int lowV = lowU + n;
    
    unsigned int sizeB = degB.get_index() - lowB.get_index();

    if (sizeA > 0 && sizeB > 0)
      for(;;)
      {
        LibEccDout(dc::polynomial|flush_cf, "A = " << cwprint(div_tct<m>(A, degA.get_index(), lowA.get_index())));
        LibEccDout(dc::polynomial|flush_cf, "B = " << cwprint(div_tct<m>(B, degB.get_index(), lowB.get_index())));
        LibEccDout(dc::polynomial|flush_cf, "U = " << cwprint(div_tct<m>(U, degU, lowU)));
        LibEccDout(dc::polynomial|flush_cf, "V = " << cwprint(div_tct<m>(V, degV, lowV)));
        if (sizeA < sizeB)
        {
          int left_shift = lowB.get_index() - lowA.get_index();
          LibEccDout(dc::polynomial|flush_cf, "B <- B + A * t^" << left_shift);
          B.xor_with_zero_padded(A, lowA.get_index(), degA.get_index(), left_shift);
          degB.find1();
          lowB.find1();
          sizeB = degB.get_index() - lowB.get_index();
          LibEccDout(dc::polynomial|flush_cf, "V <- V + U * t^" << left_shift);
          V.xor_with_zero_padded(U, lowU, degU, left_shift);
          degV = std::max(degV, degU + left_shift);
          lowV = std::min(lowV, lowU + left_shift);
          if (sizeB == 0)
            break;
        }
        else
        {
          int left_shift = lowA.get_index() - lowB.get_index();
          LibEccDout(dc::polynomial|flush_cf, "A <- A + B * t^" << left_shift);
          A.xor_with_zero_padded(B, lowB.get_index(), degB.get_index(), left_shift);
          degA.find1();
          lowA.find1();
          sizeA = degA.get_index() - lowA.get_index();
          LibEccDout(dc::polynomial|flush_cf, "U <- U + V * t^" << left_shift);
          U.xor_with_zero_padded(V, lowV, degV, left_shift);
          degU = std::max(degU, degV + left_shift);
          lowU = std::min(lowU, lowV + left_shift);
          if (sizeA == 0)
            break;
        }
      }

    LibEccDout(dc::polynomial|flush_cf, "A = " << cwprint(div_tct<m>(A, degA.get_index(), lowA.get_index())));
    LibEccDout(dc::polynomial|flush_cf, "B = " << cwprint(div_tct<m>(B, degB.get_index(), lowB.get_index())));
    LibEccDout(dc::polynomial|flush_cf, "U = " << cwprint(div_tct<m>(U, degU, lowU)));
    LibEccDout(dc::polynomial|flush_cf, "V = " << cwprint(div_tct<m>(V, degV, lowV)));

    bitset_base<m>* R;
    bitset_base<size_F>* F;
    int low1, lowR;
#if ECC_DEBUG
    int degR;
#endif
    if (sizeA == 0)
    {
      LibEccDout(dc::polynomial|flush_cf, "R = U");
      R = &U;
      F = &bitpool.FU.FU;
      low1 = lowA.get_index();
      lowR = lowU;
#if ECC_DEBUG
      degR = degU;
#endif
    }
    else // sizeB == 0
    {
      LibEccDout(dc::polynomial|flush_cf, "R = V");
      R = &V;
      F = &bitpool.FV.FV;
      low1 = lowB.get_index();
      lowR = lowV;
#if ECC_DEBUG
      degR = degV;
#endif
    }

    *F >>= low1;
    lowR -= low1;
#if ECC_DEBUG
    degR -= low1;
#endif
    // Get rid of negative exponents.
    LibEccDout(dc::polynomial|flush_cf, "lowR = " << lowR);
    LibEccDout(dc::polynomial|flush_cf, "R = " << cwprint(div_tct<m>(*R, degR, lowR)));
    if ((!k1 && k >= bitset_digit_bits) || k2 >= bitset_digit_bits)
    {
      static int const digit_shift_k2 = k2 >> bitset_digit_bits_log2;
      static int const bit_shift_k2 = k2 & (bitset_digit_bits  - 1);
      static int const digit_shift_k1 = k1 >> bitset_digit_bits_log2;
      static int const bit_shift_k1 = k1 & (bitset_digit_bits  - 1);
      static int const digit_shift_k = k >> bitset_digit_bits_log2;
      static int const bit_shift_k = k & (bitset_digit_bits  - 1);
      static int const digit_shift_m = m >> bitset_digit_bits_log2;
      static int const bit_shift_m = m & (bitset_digit_bits  - 1);
      static int const DS_minus_bit_shift_k2_with_compile_warning_evasion = (bitset_digit_bits - bit_shift_k2) & (bitset_digit_bits  - 1);
      static int const DS_minus_bit_shift_k1_with_compile_warning_evasion = (bitset_digit_bits - bit_shift_k1) & (bitset_digit_bits  - 1);
      static int const DS_minus_bit_shift_k_with_compile_warning_evasion = (bitset_digit_bits - bit_shift_k) & (bitset_digit_bits  - 1);
      static int const DS_minus_bit_shift_m_with_compile_warning_evasion = (bitset_digit_bits - bit_shift_m) & (bitset_digit_bits  - 1);
      int first_digit = (lowR + offset_F) >> bitset_digit_bits_log2;
      bitset_digit_t* ptr = F->digits_ptr() + first_digit;
      bitset_digit_t* ptr1 = R->digits_ptr();
      while(ptr < ptr1)
      {
        if (k1)
        {
          ptr[digit_shift_k2] ^= (*ptr) << bit_shift_k2;
          if (bit_shift_k2 != 0)
            ptr[digit_shift_k2 + 1] ^= (*ptr) >> DS_minus_bit_shift_k2_with_compile_warning_evasion;
          ptr[digit_shift_k1] ^= (*ptr) << bit_shift_k1;
          if (bit_shift_k1 != 0)
            ptr[digit_shift_k1 + 1] ^= (*ptr) >> DS_minus_bit_shift_k1_with_compile_warning_evasion;
        }
        ptr[digit_shift_k] ^= (*ptr) << bit_shift_k;
        if (bit_shift_k != 0)
          ptr[digit_shift_k + 1] ^= (*ptr) >> DS_minus_bit_shift_k_with_compile_warning_evasion;
        ptr[digit_shift_m] ^= (*ptr) << bit_shift_m;
        if (bit_shift_m != 0)
          ptr[digit_shift_m + 1] ^= (*ptr) >> DS_minus_bit_shift_m_with_compile_warning_evasion;
        ++ptr;
      }
    }
    else
    {
      for (unsigned int i = lowR + offset_F; i < offset_F; ++i)
      {
        if (F->test(i))
        {
#if ECC_DEBUG
          F->flip(i);           // This is not really needed, but prints nicer output below.
#endif
          if (k1)
          {
            F->flip(i + k2);
            F->flip(i + k1);
          }
          F->flip(i + k);
          F->flip(i + m);
        }
      }
    }
#if ECC_DEBUG
    lowR = 0;
    degR = 2 * m - 1;
#endif
    LibEccDout(dc::polynomial|flush_cf, "R = " << cwprint(div_tct<m>(*R, degR, lowR)));
    reduce(R->digits_ptr());
#if ECC_DEBUG
    degR = m - 1;
#endif
    LibEccDout(dc::polynomial|flush_cf, "R = " << cwprint(div_tct<m>(*R, degR, lowR)));
    *static_cast<bitset_base<m>*>(&M_coefficients) = *R;

    //std::cout << "Leaving polynomial<" << m << ", " << k << ", " << k1 << ", " << k2 << ">::operator/=()" << std::endl;
    return *this;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>&
  polynomial<m, k, k1, k2>::operator/=(typename polynomial<m, k, k1, k2>::xor_type const& expr)
  {
    return (*this /= polynomial<m, k, k1, k2>(expr));
  }

// Solve x^2 + b x = c.
// Assuming that b != 0, there are 2 solutions: x1 and x1 + b.
// This means that during the 'wiping' of the matrix in order
// to solve x, one bit of x will stay undetermined.  We need
// to take special care to make sure that this will be a bit
// for which a bit of 'b' is set, otherwise we'd return a wrong
// value.
//
// If b equals zero, then the solution is the sqrt(c).  Otherwise
// we can devide both sides of the equation by b^2 and solve
// y^2 + y = c/b^2, and set x = b * y.
//
// There will only be a solution to this equation iff 0 = Tr(c/b^2).
// (simply square the equation m-1 times and add them all up).
//
// Note that if y1 is a solution, then so is y1 + 1, hence we
// cannot determine the least significant bit of y.
//
// It is possible to compose a matrix A such that Ax = x^2 + x
// because squaring is an automorphism of the field:
// x is a sum of basis elements, ie x = b1 + b2 + b3 and
// x^2 = b1^2 + b2^2 + b3^2.  Therefore, if there exists a
// matrix S such that Sb_i = b_i^2 for any basis element then
// A = (S + I).  Moreover, such a matrix S must exist because
// there are exactly m basis elements, and a matrix of mxm
// will always be able to satisfy that.

#if ECC_DEBUG
// Debug function.
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  void polynomial<m, k, k1, k2>::print_matrix(
#if LIBECC_AUGMENTED
       bitset<2 * m> const* matrix,
#else
       bitset<m> const* matrix,
#endif
       bitset<m> const& pivotted)
  {
    // Print the matrix.
    for (unsigned int n = 1; n < m; n *= 10)
    {
      LibEccDout(dc::gaussj|continued_cf, "  ");
      for (unsigned int bit = 0; bit < matrix->number_of_bits; ++bit)
      {
        if (bit == m)
          LibEccDout(dc::continued, ' ');
        if ((bit % m) >= 1 && (bit % m) < (m + 1) / 2)
          LibEccDout(dc::continued, "+ ");
        else if (pivotted.test(bit % m))
          LibEccDout(dc::continued, (((bit % m) / n) % 10) << ' ');
        else
          LibEccDout(dc::continued, "  ");
      }
      LibEccDout(dc::finish, "");
    }
    for (unsigned int row = 0; row < m; ++row)
    {
      std::string line;
      if (row >= 1 && row < (m + 1) / 2)
        line = "+ ";
      else if (pivotted.test(row))
        line = "* ";
      else
        line = "  ";
      for (unsigned int bit = 0; bit < matrix->number_of_bits; ++bit)
      {
        if (bit == m)
          line += ' ';
        bool isset = matrix[row].test(bit);
        bool need_color = LIBECC_INPLACE && (matrix->number_of_bits > m) &&
            (((bit % m) >= 1 && (bit % m) < (m + 1) / 2) || pivotted.test(bit % m));
        if (need_color)
        {
          unsigned int corresponding_bit = (bit + m) % (2 * m);
          if (isset == matrix[row].test(corresponding_bit))
            line += "\e[32m";
          else
            line += "\e[31m";
        }
        line += (isset ? '1' : '0');
        if (need_color)
          line += "\e[0m";
        line += ' ';
      }
      LibEccDout(dc::gaussj, line);
    }
    LibEccDout(dc::gaussj|noprefix_cf, "");
  }
#endif

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  polynomial<m, k, k1, k2>::polynomial(polynomial<m, k, k1, k2> const& b, polynomial<m, k, k1, k2> const& c) :
      M_coefficients(0)
  {
    // If b == 0, then x = sqrt(c).
    if (!b.M_coefficients.any())
    {
      M_coefficients = c.M_coefficients;
      sqrt();
      return;
    }

    // Calculate c/b^2.
    bitset_digit_t b2buf[square_digits];
    polynomial<m, k, k1, k2>& b2 = b.square(b2buf);
    polynomial<m, k, k1, k2> cdb2(c);
    cdb2 /= b2;
    if (cdb2.trace() == 1)
      throw std::domain_error("x^2 + bx = c has no solution");

#if LIBECC_AUGMENTED
    typedef bitset<2 * m> matrixrow_type;
#else
    typedef bitset<m> matrixrow_type;
#endif
    static matrixrow_type matrix[m];            // A mx2m or mxm matrix.
    static bool matrix_initialized;
    if (!matrix_initialized)
    {
      std::memset(matrix, 0, sizeof(matrix));
      // Fill this matrix with either the augmented matrix (A|I) or with just A,
      // where A is the matrix such that Ax = x^2 + x.
      for (unsigned int bit = 0; bit < m; ++bit)
      {
        matrix[bit].set(bit);           // The I of A = (S + I).
#if LIBECC_AUGMENTED
        matrix[bit].set(bit + m);               // The I of (A|I).
#endif
      }
      for (unsigned int bit = 0; bit < (m + 1) / 2; ++bit)
        matrix[2 * bit].flip(bit);      // The square of low exponents.
      for (unsigned int bit = (m + 1) / 2; bit < m; ++bit)
        matrix[2 * bit - m].set(bit);   // Reduction with m.
      for (unsigned int bit = (m + 1) / 2; bit < m - k / 2; ++bit)
        matrix[2 * bit - m + k].flip(bit);      // Reduction with m - k.
      if (k1)
      {
        for (unsigned int bit = (m + 1) / 2; bit < m - k1 / 2; ++bit)
          matrix[2 * bit - m + k1].flip(bit);   // Reduction with m - k1.
        for (unsigned int bit = (m + 1) / 2; bit < m - k2 / 2; ++bit)
          matrix[2 * bit - m + k2].flip(bit);   // Reduction with m - k2.
      }
      for (unsigned int bit = m - k / 2; bit < m; ++bit)
      {
        matrix[2 * bit - m + k - m].flip(bit);
        matrix[2 * bit - m + k - m + k].flip(bit);
        if (k1)
        {
          matrix[2 * bit - m + k - m + k1].flip(bit);
          matrix[2 * bit - m + k - m + k2].flip(bit);
        }
      }
      if (k1)
      {
        for (unsigned int bit = m - k1 / 2; bit < m; ++bit)
        {
          matrix[2 * bit - m + k1 - m].flip(bit);
          matrix[2 * bit - m + k1 - m + k].flip(bit);
          matrix[2 * bit - m + k1 - m + k1].flip(bit);
          matrix[2 * bit - m + k1 - m + k2].flip(bit);
        }
        for (unsigned int bit = m - k2 / 2; bit < m; ++bit)
        {
          matrix[2 * bit - m + k2 - m].flip(bit);
          matrix[2 * bit - m + k2 - m + k].flip(bit);
          matrix[2 * bit - m + k2 - m + k1].flip(bit);
          matrix[2 * bit - m + k2 - m + k2].flip(bit);
        }
      }

      bitset<m> pivotted;
      pivotted.reset();

      LibEccDebug(if (dc::gaussj.is_on()) print_matrix(matrix, pivotted));

      // Next, wipe it, so that the left half becomes I.
      // The first half is easy.
      for (unsigned int wipecol = 1; wipecol < (m + 1) / 2; ++wipecol)
      {
        matrix[2 * wipecol] ^= matrix[wipecol];
#if LIBECC_INPLACE
        matrix[2 * wipecol].set(wipecol);               // Store the inverse in-place, destroying the original.
#endif
      }

      // The second half is not.  Use Gauss-Jordan here.
      // Note that pivotting is hardly necessary because our arithmetic is infinitely accurate,
      // but we still need to find a '1' when we encounter a '0' on the main diagonal of course.
      // There will always be at least one '1' in every column, so that partial pivotting suffices
      // (speeding up things obviously), with the exception of the case where that '1' is only
      // found in row 0 (which is our 'singular' row and needs some special attention).
      // The row swapping is done because it is needed if we want to do our work "in-place",
      // reducing the amount of memory needed with a factor of two.

      LibEccDebug(if (dc::gaussj.is_on()) print_matrix(matrix, pivotted));

      unsigned int rowswaps[m];
      rowswaps[0] = 0;
      unsigned int colswaps[m], colswaps_inverse[m];
      for (unsigned int row = 0; row < m; ++row)
      {
        colswaps[row] = row;
        colswaps_inverse[row] = row;
      }

      // Run over all remaining columns and wipe them, immedeately replacing them with the result
      // since once a column is wiped we don't need its contents anymore.  Moreover, while wiping
      // the column it is optionally swapped with another column at the same time.  This, of course,
      // is only done to make the code not understandable anymore for you.
#if LIBECC_SWAPCOLUMNS
      for (unsigned int colcnt = (m + 1) / 2; colcnt < m; ++colcnt)
#else
      for (unsigned int wipecol = (m + 1) / 2; wipecol < m; ++wipecol)
#endif
      {
#if LIBECC_SWAPCOLUMNS
        // Find the next row that wasn't already wiped.
        unsigned int wipecol = colswaps[colcnt];
#if ECC_DEBUG
        LibEccDout(dc::gaussj, "colcnt = " << colcnt);
        for (unsigned int row = 0; row < m; ++row)
        {
          LibEccDout(dc::gaussj, "colswaps[" << row << "] = " << colswaps[row] << "\t\tcolswaps_inverse[" << row << "] = " << colswaps_inverse[row]);
          assert(colswaps[colswaps_inverse[row]] == row);
          assert(colswaps_inverse[colswaps[row]] == row);
        }
        LibEccDout(dc::polynomial|noprefix_cf, "");
#endif
#endif

        // First find a suitable row to wipe with.
        // This searching is called 'pivotting'.
        LibEccDout(dc::gaussj, "Searching for suitable row to wipe with in column " << wipecol);
        unsigned int pivotrow;
        if (!matrix[wipecol].test(wipecol) || pivotted.test(wipecol))
        {
          for (pivotrow = wipecol;;)
          {
            if (++pivotrow == m)
            {
              if (matrix[0].test(wipecol) && !pivotted.template test<0>())
                pivotrow = 0;
              else
              {
                for (pivotrow = (m + 1) / 2; pivotrow < wipecol; ++pivotrow)
                  if (matrix[pivotrow].test(wipecol) && !pivotted.test(pivotrow))
                    break;
              }
              if (pivotrow == wipecol)
              {
                // This happens when we swapped with column 0 (which is all zeroes), for example when m == 14.
                // Just ignore this column.
                pivotrow = m;                   // Flag that we need to continue the main loop.
                pivotted.set(wipecol);
                matrix[wipecol].set(wipecol);   // Copy identity matrix over.
                break;
              }
            }
            if (matrix[pivotrow].test(wipecol) && !pivotted.test(pivotrow))
              break;
          }
          if (pivotrow == m)
            continue;
        }
        else
          pivotrow = wipecol;
        LibEccDout(dc::gaussj, "Using row " << pivotrow << " to wipe column " << wipecol);
        LibEccDout(dc::gaussj, "Before:");
        LibEccDebug(if (dc::gaussj.is_on()) print_matrix(matrix, pivotted));
        pivotted.set(pivotrow);
#if LIBECC_SWAPCOLUMNS
        rowswaps[colcnt] = pivotrow;
        LibEccDout(dc::gaussj, "Setting rowswaps[" << colcnt << "] to " << pivotrow);
#else
        rowswaps[wipecol] = pivotrow;                   // We temporarily use row 'pivotrow' to store row 'wipecol'.
        LibEccDout(dc::gaussj, "Setting rowswaps[" << wipecol << "] to " << pivotrow);
#endif
        if (pivotrow == wipecol)
        {
#if LIBECC_INPLACE
          matrix[pivotrow].set(wipecol);                // Store the inverse in-place, destroying the original.
#endif
          for (unsigned int row = 0; row < m; ++row)
          {
            if (row == pivotrow)
              continue;
            if (matrix[row].test(wipecol))
            {
#if LIBECC_INPLACE
              matrix[row].clear(wipecol);               // Store the inverse in-place, destroying the original.
#endif
              matrix[row] ^= matrix[pivotrow]; 
            }
          }
        }
        else
        {
          // This block contains the main magic.  It's hard to understand I am afraid.
          // Basically this does the same as the code block above, but at the same time
          // swaps the columns 'wipecol' and 'pivotrow'.

#if LIBECC_SWAPCOLUMNS
          // Swap pivot row bits, and set the bit in pivotrow (thats the identity matrix bit).
          if (matrix[pivotrow].test(pivotrow) != matrix[pivotrow].test(wipecol))
          {
            matrix[pivotrow].flip(wipecol);
#if !LIBECC_INPLACE
            matrix[pivotrow].flip(pivotrow);    // No need to flip the 'pivotrow' column when we set it in the next line.
#endif
          }
#endif
#if LIBECC_INPLACE
          matrix[pivotrow].set(pivotrow);               // Store the inverse in-place, destroying the original.
#endif
          for (unsigned int row = 0; row < m; ++row)
          {
            if (row == pivotrow)                        // Don't wipe the row that we use to wipe.
              continue;
            matrixrow_type& mrow = matrix[row];
            if (mrow.test(wipecol))
            {
#if LIBECC_SWAPCOLUMNS
              if (!mrow.test(pivotrow))         // If the value in the two columns differ,
              {
                mrow.clear(wipecol);            // swap the two values.
#if !LIBECC_INPLACE
                mrow.set(pivotrow);             // No need to set pivotrow when it is overwritten in the next line.
#endif
              }
#endif
#if LIBECC_INPLACE
              mrow.clear(pivotrow);             // Store the inverse in-place, destroying the original.
                                                // This represents a 0 from the identity matrix.
#endif
              mrow ^= matrix[pivotrow];         // Ok, now the columns have been swapped and to-be-wiped column
                                                // has been replaced with a clean identity matrix bit. Perform
                                                // the actual wiping.
            }
#if LIBECC_SWAPCOLUMNS
            else if (mrow.test(pivotrow))       // Are the pivotrow and wipecol different?
            {
              mrow.set(wipecol);                // Then flip both, exchanging them effectively. If they
              mrow.clear(pivotrow);             //   were the same, consider them exchanged anyway.
            }
#endif
          }
#if LIBECC_SWAPCOLUMNS
          LibEccDout(dc::gaussj, "Also swapped columns " << pivotrow << " and " << wipecol);
          // Keep colswaps up to date.  We need colswaps_inverse to do that, therefore
          // we need to keep colswaps_inverse up to date too.
          std::swap(colswaps[colswaps_inverse[wipecol]], colswaps[colswaps_inverse[pivotrow]]);
          std::swap(colswaps_inverse[wipecol], colswaps_inverse[pivotrow]);
#endif
        }
        LibEccDout(dc::gaussj, "After:");
        LibEccDebug(if (dc::gaussj.is_on()) print_matrix(matrix, pivotted));
      }

#if ECC_DEBUG
      for (unsigned int i = 0; i < m; ++i)
      {
        if (rowswaps[i] != i)
          LibEccDout(dc::gaussj, i << " : " << rowswaps[i]);
        // Skip the first half of the matrix.
        if (i == 0)
          i = (m + 1) / 2 - 1;
      }
      LibEccDout(dc::gaussj|noprefix_cf, "");
#endif

      if (pivotted.test(0))
      {
        int row0 = (m + 1) / 2;
        while (pivotted.test(row0))
          ++row0;
        rowswaps[0] = row0;
        pivotted.set(row0);
      }

      // Next perform some row rotations, in order to get all rows on their correct places again.
      for (unsigned int i = 0; i < m; ++i)
      {
        if (rowswaps[i] != i)
        {
          unsigned int j = i;
          bitset<2 * m> temp = matrix[j];
          LibEccDout(dc::gaussj|continued_cf, j);
          do
          {
            matrix[j] = matrix[rowswaps[j]];
            LibEccDout(dc::continued, " <-- " << rowswaps[j]);
            j = rowswaps[j];
          }
          while (rowswaps[j] != i);
          matrix[j] = temp;
          LibEccDout(dc::finish, " <-- " << i);
          j = i;
          do
          {
            int pj = j;
            j = rowswaps[pj];
            // Update the administration so that we won't try to rotate them again.
            rowswaps[pj] = pj;
          }
          while (j != i);
        }
        // Skip the first half of the matrix.
        if (i == 0)
          i = (m + 1) / 2 - 1;
      }

      LibEccDebug(if (dc::gaussj.is_on()) print_matrix(matrix, pivotted));
      matrix_initialized = true;
    }

    // Multiply the matrix with cdb2.
    for (unsigned int row = 0; row < m; ++row)
    {
#if LIBECC_AUGMENTED
#if LIBECC_INPLACE
      bitset<m> tmp = matrix[row];
#else
      bitset<2 * m> tmp2;
      matrix[row].template shift_op<m, right, assign>(tmp2);
      bitset<m> tmp = tmp2;
#endif
      tmp &= cdb2.get_bitset();
#else
      bitset<m> tmp = matrix[row] & cdb2.get_bitset();
#endif
      if (tmp.odd())
        M_coefficients.set(row);
    }

    // Finally, multiply with b to get x.
    *this *= b;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline bool
  operator==(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2> const& p2)
  {
    return p1.M_coefficients == p2.M_coefficients;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline bool
  operator==(typename polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2)
  {
    return polynomial<m, k, k1, k2>(expr) == p2;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline bool
  operator==(polynomial<m, k, k1, k2> const& p1, typename polynomial<m, k, k1, k2>::xor_type const& expr)
  {
    return p1 == polynomial<m, k, k1, k2>(expr);
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline bool
  operator!=(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2> const& p2)
  {
    return p1.M_coefficients != p2.M_coefficients;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline bool
  operator!=(typename polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2)
  {
    return polynomial<m, k, k1, k2>(expr) != p2;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline bool
  operator!=(polynomial<m, k, k1, k2> const& p1, typename polynomial<m, k, k1, k2>::xor_type const& expr)
  {
    return p1 != polynomial<m, k, k1, k2>(expr);
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline typename polynomial<m, k, k1, k2>::xor_type
  operator+(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2> const& p2)
  {
    return typename polynomial<m, k, k1, k2>::xor_type(p1.M_coefficients, p2.M_coefficients);
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline typename polynomial<m, k, k1, k2>::xor_type
  operator-(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2> const& p2)
  {
    return typename polynomial<m, k, k1, k2>::xor_type(p1.M_coefficients, p2.M_coefficients);
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>
  operator*(polynomial<m, k, k1, k2> const& p1, polynomial<m, k, k1, k2> const& p2)
  {
    polynomial<m, k, k1, k2> result;
    p1.multiply_with(p2, result.M_coefficients);
    return result;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>
  operator*(typename polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2)
  {
    return polynomial<m, k, k1, k2>(expr) * p2;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>
  operator*(polynomial<m, k, k1, k2> const& p1, typename polynomial<m, k, k1, k2>::xor_type const& expr)
  {
    return p1 * polynomial<m, k, k1, k2>(expr);
  }


template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>
  operator/(polynomial<m, k, k1, k2> const& e1, polynomial<m, k, k1, k2> const& e2)
  {
    polynomial<m, k, k1, k2> tmp(e1);
    tmp /= e2;
    return tmp;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>
  operator/(typename polynomial<m, k, k1, k2>::xor_type const& expr, polynomial<m, k, k1, k2> const& p2)
  {
    return polynomial<m, k, k1, k2>(expr) / p2;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  inline polynomial<m, k, k1, k2>
  operator/(polynomial<m, k, k1, k2> const& p1, typename polynomial<m, k, k1, k2>::xor_type const& expr)
  {
    return p1 / polynomial<m, k, k1, k2>(expr);
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  std::ostream& operator<<(std::ostream& os, polynomial<m, k, k1, k2> const& p)
  {
    p.M_coefficients.base2_print_on(os);
    return os;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  std::ostream& operator<<(std::ostream& os, typename polynomial<m, k, k1, k2>::xor_type const& expr)
  {
    polynomial<m, k, k1, k2> p(expr);
    p.M_coefficients.base2_print_on(os);
    return os;
  }

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  bool polynomial<m, k, k1, k2>::S_normal_initialized;

template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  bitset<m> polynomial<m, k, k1, k2>::S_normal;
 
template<unsigned int m, unsigned int k, unsigned int k1, unsigned int k2>
  void polynomial<m, k, k1, k2>::calculate_normal(void)
  {
#if 0
    bitset<m> single_bit(1);
    polynomial trace;
    bitset_digit_t nextfrob1_buf[square_digits];
    bitset_digit_t nextfrob2_buf[square_digits];
    polynomial* nextfrob1;
    polynomial* nextfrob2;
    for (int bit = 0; bit < m; ++bit)
    {
      trace = single_bit;
      nextfrob1 = &trace.square(nextfrob1_buf);
      for (int i = 0; i < (m - 1) / 2; ++i)
      {
        nextfrob2 = &nextfrob1->square(nextfrob2_buf);
        trace += *nextfrob1 + *nextfrob2;
        if ((m & 1) && i == (m - 3) / 2)
          break;
        nextfrob1 = &nextfrob2->square(nextfrob1_buf);
      }
      if (!(m & 1))
        trace += *nextfrob1;
      if (trace.get_bitset().template test<0>())
        S_normal.set(bit);
      single_bit.template shift_op<1, libecc::left, libecc::assign>(single_bit);
    }
#else
    // We can do that faster... I didn't prove this yet, but it works.
    if ((m & 1))
      S_normal.template set<0>();
    if (((m - k) & 1))
      S_normal.template set<m - k>();
    if (k1)
    {
      if (((m - k1) & 1))
        S_normal.template set<m - k1>();
      if (((m - k2) & 1))
        S_normal.template set<m - k2>();
    }
#endif
    S_normal_initialized = true;
  }

} // namespace libecc

#include <libecc/square.hcc>    // File with different copyright.

#endif // LIBECC_POLYNOMIAL_H