modp.hpp

#ifndef __MODP_HPP_
#define __MODP_HPP_

#include <stdexcept>
using std::domain_error;

#include <vector>
using std::vector;

#include "mod.hpp"

namespace Rings {

template<typename T, T p>
class Modp : virtual public Field_Element, virtual public Mod<T, p> {

protected:

  using Mod<T,p>::value;
  using Mod<T,p>::invertible;

  static vector<T> inverses;

  virtual void check_inverse() noexcept override { }
  void check_inverses();

public:

  using Field_Element::has_inverse;

  Modp<T,p>();
  Modp<T,p>(T);
  Modp<T,p>(const Mod<T,p> & other);
  Modp<T,p>(const Modp<T,p> & other);

  virtual Modp<T,p> & operator * (const T & other) const override;
  virtual Modp<T,p> & operator / (const Field_Element & other) const override;
  virtual Modp<T,p> & inverse() const override;

  virtual bool has_inverse() const override;

};

template <typename T, T p>
void Modp<T,p>::check_inverses() {
  if (inverses.size() == 0) {
    inverses.resize(p);
    inverses[1] = 1;
    inverses[p - 1] = p - 1;
    for (T i = 2; i < p - 1; ++i) {
      T a, b;
      T d = gcd<T>(i, p, a, b);
      if (d != 1)
          throw domain_error("Attempted a non-prime modulus for Modp");
      if (a < 0) a += p;
      inverses[i] = a;
    }
  }
}

template <typename T, T p>
Modp<T,p>::Modp() : Mod<T,p>() { invertible = true; check_inverses(); }

template <typename T, T p>
Modp<T,p>::Modp(T v) : Mod<T,p>(v) {
  invertible = true;
  check_inverses();
}

template <typename T, T p>
Modp<T,p>::Modp(const Modp<T,p> & other) : Mod<T,p>(other) {
  invertible = true;
  check_inverses();
}

template <typename T, T p>
Modp<T,p>::Modp(const Mod<T,p> & other) : Mod<T,p>(other) {
  invertible = true;
  check_inverses();
}

template <typename T, T p>
Modp<T,p> & Modp<T,p>::operator * (const T & other) const {
  static Modp<T,p> result(value * other);
  return result;
}

template <typename T, T p>
Modp<T,p> & Modp<T,p>::operator / (const Field_Element & other)
    const
{
  static Modp<T,p> result;
  result.value
      = value * inverses[dynamic_cast<const Modp<T,p> &>(other).value];
  result.adjust_value();
  return result;
}

template <typename T, T p>
Modp<T,p> & Modp<T,p>::inverse() const {
  static Modp<T,p> result;
  result.value = inverses[value];
  return result;
}

template <typename T, T p>
bool Modp<T,p>::has_inverse() const { return not is_zero(); }

template<typename T, T p> vector<T> Modp<T,p>::inverses;

template<typename T, T p>
constexpr Modp<T,p> initialize_inverses;

}

#endif