cln_mirror/include/cl_modinteger.h

// Modular integer operations.

#ifndef _CL_MODINTEGER_H
#define _CL_MODINTEGER_H

#include "cl_object.h"
#include "cl_ring.h"
#include "cl_integer.h"
#include "cl_random.h"
#include "cl_malloc.h"
#include "cl_io.h"
#include "cl_proplist.h"
#include "cl_condition.h"
#include "cl_abort.h"
#undef random // Linux defines random() as a macro!


// Representation of an element of a ring Z/mZ.

// To protect against mixing elements of different modular rings, such as
// (3 mod 4) + (2 mod 5), every modular integer carries its ring in itself.


// Representation of a ring Z/mZ.

class cl_heap_modint_ring;

class cl_modint_ring : public cl_ring {
public:
	// Default constructor.
	cl_modint_ring ();
	// Constructor. Takes a cl_heap_modint_ring*, increments its refcount.
	cl_modint_ring (cl_heap_modint_ring* r);
	// Copy constructor.
	cl_modint_ring (const cl_modint_ring&);
	// Assignment operator.
	cl_modint_ring& operator= (const cl_modint_ring&);
	// Automatic dereferencing.
	cl_heap_modint_ring* operator-> () const
	{ return (cl_heap_modint_ring*)heappointer; }
};

// Z/0Z
extern const cl_modint_ring cl_modint0_ring;
// Default constructor. This avoids dealing with NULL pointers.
inline cl_modint_ring::cl_modint_ring ()
	: cl_ring (as_cl_private_thing(cl_modint0_ring)) {}
CL_REQUIRE(cl_MI)
// Copy constructor and assignment operator.
CL_DEFINE_COPY_CONSTRUCTOR2(cl_modint_ring,cl_ring)
CL_DEFINE_ASSIGNMENT_OPERATOR(cl_modint_ring,cl_modint_ring)

// Normal constructor for `cl_modint_ring'.
inline cl_modint_ring::cl_modint_ring (cl_heap_modint_ring* r)
	: cl_ring ((cl_private_thing) (cl_inc_pointer_refcount((cl_heap*)r), r)) {}

// Operations on modular integer rings.

inline bool operator== (const cl_modint_ring& R1, const cl_modint_ring& R2)
{ return (R1.pointer == R2.pointer); }
inline bool operator!= (const cl_modint_ring& R1, const cl_modint_ring& R2)
{ return (R1.pointer != R2.pointer); }
inline bool operator== (const cl_modint_ring& R1, cl_heap_modint_ring* R2)
{ return (R1.pointer == R2); }
inline bool operator!= (const cl_modint_ring& R1, cl_heap_modint_ring* R2)
{ return (R1.pointer != R2); }


// Condition raised when a probable prime is discovered to be composite.
struct cl_composite_condition : public cl_condition {
	SUBCLASS_cl_condition()
	cl_I p; // the non-prime
	cl_I factor; // a nontrivial factor, or 0
	// Constructors.
	cl_composite_condition (const cl_I& _p)
		: p (_p), factor (0)
		{ print(cl_stderr); }
	cl_composite_condition (const cl_I& _p, const cl_I& _f)
		: p (_p), factor (_f)
		{ print(cl_stderr); }
	// Implement general condition methods.
	const char * name () const;
	void print (cl_ostream) const;
	~cl_composite_condition () {}
};


// Representation of an element of a ring Z/mZ.

class _cl_MI /* cf. _cl_ring_element */ {
public:
	cl_I rep;		// representative, integer >=0, <m
				// (maybe the Montgomery representative!)
	// Default constructor.
	_cl_MI () : rep () {}
public: /* ugh */
	// Constructor.
	_cl_MI (const cl_heap_modint_ring* R, const cl_I& r) : rep (r) { (void)R; }
	_cl_MI (const cl_modint_ring& R, const cl_I& r) : rep (r) { (void)R; }
public:
	// Conversion.
	CL_DEFINE_CONVERTER(_cl_ring_element)
public:	// Ability to place an object at a given address.
	void* operator new (size_t size) { return cl_malloc_hook(size); }
	void* operator new (size_t size, _cl_MI* ptr) { (void)size; return ptr; }
	void operator delete (void* ptr) { cl_free_hook(ptr); }
};

class cl_MI /* cf. cl_ring_element */ : public _cl_MI {
protected:
	cl_modint_ring _ring;	// ring Z/mZ
public:
	const cl_modint_ring& ring () const { return _ring; }
	// Default constructor.
	cl_MI () : _cl_MI (), _ring () {}
public: /* ugh */
	// Constructor.
	cl_MI (const cl_modint_ring& R, const cl_I& r) : _cl_MI (R,r), _ring (R) {}
	cl_MI (const cl_modint_ring& R, const _cl_MI& r) : _cl_MI (r), _ring (R) {}
public:
	// Conversion.
	CL_DEFINE_CONVERTER(cl_ring_element)
	// Debugging output.
	void debug_print () const;
public:	// Ability to place an object at a given address.
	void* operator new (size_t size) { return cl_malloc_hook(size); }
	void* operator new (size_t size, cl_MI* ptr) { (void)size; return ptr; }
	void operator delete (void* ptr) { cl_free_hook(ptr); }
};


// Representation of an element of a ring Z/mZ or an exception.

class cl_MI_x {
private:
	cl_MI value;
public:
	cl_composite_condition* condition;
	// Constructors.
	cl_MI_x (cl_composite_condition* c) : value (), condition (c) {}
	cl_MI_x (const cl_MI& x) : value (x), condition (NULL) {}
	// Cast operators.
	//operator cl_MI& () { if (condition) cl_abort(); return value; }
	//operator const cl_MI& () const { if (condition) cl_abort(); return value; }
	operator cl_MI () const { if (condition) cl_abort(); return value; }
};


// Ring operations.

struct _cl_modint_setops /* cf. _cl_ring_setops */ {
	// print
	void (* fprint) (cl_heap_modint_ring* R, cl_ostream stream, const _cl_MI& x);
	// equality
	cl_boolean (* equal) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
	// random number
	const _cl_MI (* random) (cl_heap_modint_ring* R, cl_random_state& randomstate);
};
struct _cl_modint_addops /* cf. _cl_ring_addops */ {
	// 0
	const _cl_MI (* zero) (cl_heap_modint_ring* R);
	cl_boolean (* zerop) (cl_heap_modint_ring* R, const _cl_MI& x);
	// x+y
	const _cl_MI (* plus) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
	// x-y
	const _cl_MI (* minus) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
	// -x
	const _cl_MI (* uminus) (cl_heap_modint_ring* R, const _cl_MI& x);
};
struct _cl_modint_mulops /* cf. _cl_ring_mulops */ {
	// 1
	const _cl_MI (* one) (cl_heap_modint_ring* R);
	// canonical homomorphism
	const _cl_MI (* canonhom) (cl_heap_modint_ring* R, const cl_I& x);
	// x*y
	const _cl_MI (* mul) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
	// x^2
	const _cl_MI (* square) (cl_heap_modint_ring* R, const _cl_MI& x);
	// x^y, y Integer >0
	const _cl_MI (* expt_pos) (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y);
	// x^-1
	const cl_MI_x (* recip) (cl_heap_modint_ring* R, const _cl_MI& x);
	// x*y^-1
	const cl_MI_x (* div) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
	// x^y, y Integer
	const cl_MI_x (* expt) (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y);
	// x -> x mod m for x>=0
	const cl_I (* reduce_modulo) (cl_heap_modint_ring* R, const cl_I& x);
	// some inverse of canonical homomorphism
	const cl_I (* retract) (cl_heap_modint_ring* R, const _cl_MI& x);
};
#if defined(__GNUC__) && (__GNUC__ == 2) && (__GNUC_MINOR__ < 8) // workaround two g++-2.7.0 bugs
  #define cl_modint_setops  _cl_modint_setops
  #define cl_modint_addops  _cl_modint_addops
  #define cl_modint_mulops  _cl_modint_mulops
#else
  typedef const _cl_modint_setops  cl_modint_setops;
  typedef const _cl_modint_addops  cl_modint_addops;
  typedef const _cl_modint_mulops  cl_modint_mulops;
#endif

// Representation of the ring Z/mZ.

// Currently rings are garbage collected only when they are not referenced
// any more and when the ring table gets full.

// Modular integer rings are kept unique in memory. This way, ring equality
// can be checked very efficiently by a simple pointer comparison.

class cl_heap_modint_ring /* cf. cl_heap_ring */ : public cl_heap {
	SUBCLASS_cl_heap_ring()
private:
	cl_property_list properties;
protected:
	cl_modint_setops* setops;
	cl_modint_addops* addops;
	cl_modint_mulops* mulops;
public:
	cl_I modulus;	// m, normalized to be >= 0
public:
	// Low-level operations.
	void _fprint (cl_ostream stream, const _cl_MI& x)
		{ setops->fprint(this,stream,x); }
	cl_boolean _equal (const _cl_MI& x, const _cl_MI& y)
		{ return setops->equal(this,x,y); }
	const _cl_MI _random (cl_random_state& randomstate)
		{ return setops->random(this,randomstate); }
	const _cl_MI _zero ()
		{ return addops->zero(this); }
	cl_boolean _zerop (const _cl_MI& x)
		{ return addops->zerop(this,x); }
	const _cl_MI _plus (const _cl_MI& x, const _cl_MI& y)
		{ return addops->plus(this,x,y); }
	const _cl_MI _minus (const _cl_MI& x, const _cl_MI& y)
		{ return addops->minus(this,x,y); }
	const _cl_MI _uminus (const _cl_MI& x)
		{ return addops->uminus(this,x); }
	const _cl_MI _one ()
		{ return mulops->one(this); }
	const _cl_MI _canonhom (const cl_I& x)
		{ return mulops->canonhom(this,x); }
	const _cl_MI _mul (const _cl_MI& x, const _cl_MI& y)
		{ return mulops->mul(this,x,y); }
	const _cl_MI _square (const _cl_MI& x)
		{ return mulops->square(this,x); }
	const _cl_MI _expt_pos (const _cl_MI& x, const cl_I& y)
		{ return mulops->expt_pos(this,x,y); }
	const cl_MI_x _recip (const _cl_MI& x)
		{ return mulops->recip(this,x); }
	const cl_MI_x _div (const _cl_MI& x, const _cl_MI& y)
		{ return mulops->div(this,x,y); }
	const cl_MI_x _expt (const _cl_MI& x, const cl_I& y)
		{ return mulops->expt(this,x,y); }
	const cl_I _reduce_modulo (const cl_I& x)
		{ return mulops->reduce_modulo(this,x); }
	const cl_I _retract (const _cl_MI& x)
		{ return mulops->retract(this,x); }
	// High-level operations.
	void fprint (cl_ostream stream, const cl_MI& x)
	{
		if (!(x.ring() == this)) cl_abort();
		_fprint(stream,x);
	}
	cl_boolean equal (const cl_MI& x, const cl_MI& y)
	{
		if (!(x.ring() == this)) cl_abort();
		if (!(y.ring() == this)) cl_abort();
		return _equal(x,y);
	}
	const cl_MI random (cl_random_state& randomstate = cl_default_random_state)
	{
		return cl_MI(this,_random(randomstate));
	}
	const cl_MI zero ()
	{
		return cl_MI(this,_zero());
	}
	cl_boolean zerop (const cl_MI& x)
	{
		if (!(x.ring() == this)) cl_abort();
		return _zerop(x);
	}
	const cl_MI plus (const cl_MI& x, const cl_MI& y)
	{
		if (!(x.ring() == this)) cl_abort();
		if (!(y.ring() == this)) cl_abort();
		return cl_MI(this,_plus(x,y));
	}
	const cl_MI minus (const cl_MI& x, const cl_MI& y)
	{
		if (!(x.ring() == this)) cl_abort();
		if (!(y.ring() == this)) cl_abort();
		return cl_MI(this,_minus(x,y));
	}
	const cl_MI uminus (const cl_MI& x)
	{
		if (!(x.ring() == this)) cl_abort();
		return cl_MI(this,_uminus(x));
	}
	const cl_MI one ()
	{
		return cl_MI(this,_one());
	}
	const cl_MI canonhom (const cl_I& x)
	{
		return cl_MI(this,_canonhom(x));
	}
	const cl_MI mul (const cl_MI& x, const cl_MI& y)
	{
		if (!(x.ring() == this)) cl_abort();
		if (!(y.ring() == this)) cl_abort();
		return cl_MI(this,_mul(x,y));
	}
	const cl_MI square (const cl_MI& x)
	{
		if (!(x.ring() == this)) cl_abort();
		return cl_MI(this,_square(x));
	}
	const cl_MI expt_pos (const cl_MI& x, const cl_I& y)
	{
		if (!(x.ring() == this)) cl_abort();
		return cl_MI(this,_expt_pos(x,y));
	}
	const cl_MI_x recip (const cl_MI& x)
	{
		if (!(x.ring() == this)) cl_abort();
		return _recip(x);
	}
	const cl_MI_x div (const cl_MI& x, const cl_MI& y)
	{
		if (!(x.ring() == this)) cl_abort();
		if (!(y.ring() == this)) cl_abort();
		return _div(x,y);
	}
	const cl_MI_x expt (const cl_MI& x, const cl_I& y)
	{
		if (!(x.ring() == this)) cl_abort();
		return _expt(x,y);
	}
	const cl_I reduce_modulo (const cl_I& x)
	{
		return _reduce_modulo(x);
	}
	const cl_I retract (const cl_MI& x)
	{
		if (!(x.ring() == this)) cl_abort();
		return _retract(x);
	}
	// Miscellaneous.
	sintL bits; // number of bits needed to represent a representative, or -1
	int log2_bits; // log_2(bits), or -1
	// Property operations.
	cl_property* get_property (const cl_symbol& key)
		{ return properties.get_property(key); }
	void add_property (cl_property* new_property)
		{ properties.add_property(new_property); }
// Constructor.
	cl_heap_modint_ring (cl_I m, cl_modint_setops*, cl_modint_addops*, cl_modint_mulops*);
// This class is intented to be subclassable, hence needs a virtual destructor.
	virtual ~cl_heap_modint_ring () {}
private:
	virtual void dummy ();
};
#define SUBCLASS_cl_heap_modint_ring() \
  SUBCLASS_cl_heap_ring()

// Lookup or create a modular integer ring  Z/mZ
extern const cl_modint_ring cl_find_modint_ring (const cl_I& m);
CL_REQUIRE(cl_MI)

// Runtime typing support.
extern cl_class cl_class_modint_ring;


// Operations on modular integers.

// Output.
inline void fprint (cl_ostream stream, const cl_MI& x)
	{ x.ring()->fprint(stream,x); }
CL_DEFINE_PRINT_OPERATOR(cl_MI)

// Add.
inline const cl_MI operator+ (const cl_MI& x, const cl_MI& y)
	{ return x.ring()->plus(x,y); }
inline const cl_MI operator+ (const cl_MI& x, const cl_I& y)
	{ return x.ring()->plus(x,x.ring()->canonhom(y)); }
inline const cl_MI operator+ (const cl_I& x, const cl_MI& y)
	{ return y.ring()->plus(y.ring()->canonhom(x),y); }

// Negate.
inline const cl_MI operator- (const cl_MI& x)
	{ return x.ring()->uminus(x); }

// Subtract.
inline const cl_MI operator- (const cl_MI& x, const cl_MI& y)
	{ return x.ring()->minus(x,y); }
inline const cl_MI operator- (const cl_MI& x, const cl_I& y)
	{ return x.ring()->minus(x,x.ring()->canonhom(y)); }
inline const cl_MI operator- (const cl_I& x, const cl_MI& y)
	{ return y.ring()->minus(y.ring()->canonhom(x),y); }

// Shifts.
extern const cl_MI operator<< (const cl_MI& x, sintL y); // assume 0 <= y < 2^31
extern const cl_MI operator>> (const cl_MI& x, sintL y); // assume m odd, 0 <= y < 2^31

// Equality.
inline bool operator== (const cl_MI& x, const cl_MI& y)
	{ return x.ring()->equal(x,y); }
inline bool operator!= (const cl_MI& x, const cl_MI& y)
	{ return !x.ring()->equal(x,y); }
inline bool operator== (const cl_MI& x, const cl_I& y)
	{ return x.ring()->equal(x,x.ring()->canonhom(y)); }
inline bool operator!= (const cl_MI& x, const cl_I& y)
	{ return !x.ring()->equal(x,x.ring()->canonhom(y)); }
inline bool operator== (const cl_I& x, const cl_MI& y)
	{ return y.ring()->equal(y.ring()->canonhom(x),y); }
inline bool operator!= (const cl_I& x, const cl_MI& y)
	{ return !y.ring()->equal(y.ring()->canonhom(x),y); }

// Compare against 0.
inline cl_boolean zerop (const cl_MI& x)
	{ return x.ring()->zerop(x); }

// Multiply.
inline const cl_MI operator* (const cl_MI& x, const cl_MI& y)
	{ return x.ring()->mul(x,y); }

// Squaring.
inline const cl_MI square (const cl_MI& x)
	{ return x.ring()->square(x); }

// Exponentiation x^y, where y > 0.
inline const cl_MI expt_pos (const cl_MI& x, const cl_I& y)
	{ return x.ring()->expt_pos(x,y); }

// Reciprocal.
inline const cl_MI recip (const cl_MI& x)
	{ return x.ring()->recip(x); }

// Division.
inline const cl_MI div (const cl_MI& x, const cl_MI& y)
	{ return x.ring()->div(x,y); }
inline const cl_MI div (const cl_MI& x, const cl_I& y)
	{ return x.ring()->div(x,x.ring()->canonhom(y)); }
inline const cl_MI div (const cl_I& x, const cl_MI& y)
	{ return y.ring()->div(y.ring()->canonhom(x),y); }

// Exponentiation x^y.
inline const cl_MI expt (const cl_MI& x, const cl_I& y)
	{ return x.ring()->expt(x,y); }

// Scalar multiplication.
inline const cl_MI operator* (const cl_I& x, const cl_MI& y)
	{ return y.ring()->mul(y.ring()->canonhom(x),y); }
inline const cl_MI operator* (const cl_MI& x, const cl_I& y)
	{ return x.ring()->mul(x.ring()->canonhom(y),x); }

// TODO: implement gcd, index (= gcd), unitp, sqrtp


// Debugging support.
#ifdef CL_DEBUG
extern int cl_MI_debug_module;
static void* const cl_MI_debug_dummy[] = { &cl_MI_debug_dummy,
	&cl_MI_debug_module
};
#endif


#endif /* _CL_MODINTEGER_H */