libecc/point.h

Go to the documentation of this file.

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

#include <gmpxx.h>
#include <libecc/polynomial.h>
#include <libecc/rds.h>

namespacelibecc {

enum order_algorithm {
  order_brute_force,
  order_default = order_brute_force
};

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  classpoint {
  public:
    typedef Polynomial polynomial_type;
  private:
    polynomial_type M_x;        // x coordinate of the point.
    polynomial_type M_y;        // y coordinate of the point.
    bool M_zero;                // Set when this is the 'point at infinity', in which case M_x and M_y are meaningless.
  public:
    point(void) : M_zero(true) { }
    point(std::string x, std::string y) : M_x(x), M_y(y), M_zero(false) { }
    point(polynomial_type const& x, polynomial_type const& y) : M_x(x), M_y(y), M_zero(false) { }
    point& operator=(point const& p1);
    point& operator+=(point const& p1);
    void MULTIPLY_and_assign(point const& pnt, mpz_class const& scalar);
    mpz_class order(order_algorithm algorithm = order_default) const;
    bool operator==(point const& p1) const;
    bool operator!=(point const& p1) const;
    bool check(void) const;
    polynomial_type const& get_x(void) const{ return M_x; }
    polynomial_type const& get_y(void) const{ return M_y; }
    bool is_zero(void) const{ return M_zero; }
    void print_on(std::ostream& os) const;
    void randomize(rds& random_source);
  };

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  void
  point<Polynomial, a, b>::randomize(rds& random_source)
  {
    M_zero = false;
    for(bool notdone = true; notdone;)
    {
      random_source.read(M_x.get_bitset());
      bitset_digit_t x2buf[polynomial_type::square_digits];
      polynomial_type& x2 = M_x.square(x2buf);                    // x^2
      polynomial_type r((M_x + a) * x2 + b);                      // x^3 + a x^2 + b
      try
      {
        M_y = polynomial_type(M_x, r);
        notdone = false;
      }
      catch(std::domain_error const&)
      {
      }
    }
    bitset_digit_t digit;
    random_source.read(&digit, 1);
    if ((digit & 1))
    {
      if (!M_x.get_bitset().any())
        M_zero = true;
      else
        M_y += M_x;
    }
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  bool
  point<Polynomial, a, b>::check(void) const
 {
    if (M_zero)
      return true;
    bitset_digit_t x2buf[polynomial_type::square_digits];
    polynomial_type& x2 = M_x.square(x2buf);                    // x^2
    polynomial_type r((M_x + a) * x2 + b);                      // x^3 + a x^2 + b
    bitset_digit_t y2buf[polynomial_type::square_digits];
    polynomial_type& y2 = M_y.square(y2buf);
    polynomial_type l(y2 + M_x * M_y);                          // y^2 + x y
    return (r == l);
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  inline bool
  point<Polynomial, a, b>::operator!=(point const& p1) const
 {
    return (M_zero != p1.M_zero || (!M_zero && (M_x != p1.M_x || M_y != p1.M_y)));
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  inline bool
  point<Polynomial, a, b>::operator==(point const& p1) const
 {
    return !operator!=(p1);
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  std::ostream&
  operator<<(std::ostream& os, point<Polynomial, a, b> const& p1)
  {
    if (p1.is_zero())
      os << "O";
    else
      os << '(' << p1.get_x() << ", " << p1.get_y() << ')';
    return os;
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  void
  point<Polynomial, a, b>::print_on(std::ostream& os) const
 {
    if (M_zero)
      os << "O";
    else
      os << '(' << M_x << ", " << M_y << ')';
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  point<Polynomial, a, b>&
  point<Polynomial, a, b>::operator=(point<Polynomial, a, b> const& p1)
  {
    M_zero = p1.M_zero;
    if (!M_zero)
    {
      M_x = p1.M_x;
      M_y = p1.M_y;
    }
    return *this;
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  void
  point<Polynomial, a, b>::MULTIPLY_and_assign(point<Polynomial, a, b> const& pnt, mpz_class const& scalar)
  {
    point tmp(pnt);
    M_zero = true;
    unsigned long int p = 0;
    size_t end = mpz_sizeinbase(scalar.get_mpz_t(), 2);
    for (unsigned long int bit = mpz_scan1(scalar.get_mpz_t(), 0); bit < end; bit = mpz_scan1(scalar.get_mpz_t(), bit + 1))
    {
      for (; p < bit; ++p)
        tmp += tmp;
      *this += tmp;
    }
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  point<Polynomial, a, b>&
  point<Polynomial, a, b>::operator+=(point<Polynomial, a, b> const& p1)
  {
    if (!p1.M_zero)
    {
      if (M_zero)
      {
        M_x = p1.M_x;
        M_y = p1.M_y;
        M_zero = false;
      }
      else if (M_x != p1.M_x)
      {
        polynomial_type s = M_y + p1.M_y;
        s /= M_x + p1.M_x;
        bitset_digit_t s2buf[polynomial_type::square_digits];
        polynomial_type& new_x = s.square(s2buf);
        new_x += s + M_x;
        new_x += p1.M_x + a;
        M_y += s * (M_x + new_x) + new_x;
        M_x = new_x;
      }
      else if (M_y == p1.M_y && M_x.get_bitset().any())
      {
        polynomial_type s = M_x + M_y / M_x;
        bitset_digit_t s2buf[polynomial_type::square_digits];
        polynomial_type& new_x = s.square(s2buf);
        new_x += s + a;
        bitset_digit_t x2buf[polynomial_type::square_digits];
        polynomial_type& x2 = M_x.square(x2buf);
        M_y = x2 + new_x * (s + polynomial_type::unity());
        M_x = new_x;
      }
      else
        M_zero = true;
    }
    return *this;
  }

template<typename Polynomial, Polynomial const& a, Polynomial const& b>
  mpz_class
  point<Polynomial, a, b>::order(order_algorithm algorithm) const
 {
    mpz_class n = 1;
    switch(algorithm)
    {
      case order_brute_force:
      {
        point nP(*this);
        while(!nP.M_zero && nP.M_x.get_bitset().any())
        {
          nP += *this;
          ++n;
        }
        if (!nP.M_zero)
          n += n;
        break;
      }
    }
    return n;
  }

} // namespace libecc

#endif // LIBECC_POINT_H