libecc/fieldmath.h

Go to the documentation of this file.

//
//
// This file is part of the libecc package.
// Copyright (C) 2002, 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_FIELDMATH_H
#define LIBECC_FIELDMATH_H

#include <gmpxx.h>
#include <libecc/bitset.h>

namespacelibecc {

template<class field_type>
  classmultiplicative_square {
    public:
      void operator()(field_type& field_element) const{ field_element *= field_element; }
  };

template<class field_type, class square_functor>
  field_type exponentiation(field_type const& base, mpz_class const& exponent,
      square_functor const& do_square = multiplicative_square<field_type>())
  {
    field_type result(mpz_tstbit(exponent.get_mpz_t(), 0) ? base : field_type::unity());
    field_type current_factor(base);
    unsigned long next_bit = 1;
    for(;;)
    {
      unsigned long current_bit = mpz_scan1(exponent.get_mpz_t(), next_bit);
      if (current_bit == (unsigned long)-1)
        break;
      for(int cnt = current_bit - next_bit; cnt >= 0; --cnt)
        do_square(current_factor);
      result *= current_factor;
      next_bit = current_bit + 1;
    }
    return result;
  }

template<unsigned int m>
  bitset<m>&
  gcd(bitset<m>& polynomial_coefficients0, bitset<m>& polynomial_coefficients1)
  {
    bitset<m>* polycoef[2];
    polycoef[0] = &polynomial_coefficients0;
    polycoef[1] = &polynomial_coefficients1;
    typename bitset<m>::const_reverse_iterator iter[2];

    // Determine the highest degree bit that is set, in the two polynomial_coefficients.
    for (int p = 0; p < 2; ++p)
    {
      iter[p] = polycoef[p]->rbegin();
      iter[p].find1();
      if (iter[p] == polycoef[p]->rend())
        return *polycoef[1 - p];                // The other polynomials is zero.
    }
    int distance = (iter[0] - iter[1]);
    if (distance < 0)
      distance = -distance;
    int smallest = (iter[0] > iter[1]) ? 0 : 1;
    for (int largest = 1 - smallest;; largest = 1 - largest, smallest = 1 - smallest)
    {
      bitset<m> temp(*polycoef[smallest]);
      temp <<= distance;
      int last_shift = distance;
      *polycoef[largest] ^= temp;
      for(int shift = distance - 1; shift >= 0; --shift)
      {
        if (!*(iter[largest] + (distance - shift)))
          continue;
        temp >>= (last_shift - shift);
        last_shift = shift;
        *polycoef[largest] ^= temp;
      }
      iter[largest].find1();
      if (iter[largest] == polycoef[largest]->rend())
        break;                          // Found the gcd.
      distance = (iter[largest] - iter[smallest]);
    }
    return *polycoef[smallest];
  }
  
} // namespace libecc

#endif // LIBECC_FIELDMATH_H