Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-19 08:08:31

0001 // Modular integer operations.
0002 
0003 #ifndef _CL_MODINTEGER_H
0004 #define _CL_MODINTEGER_H
0005 
0006 #include "cln/object.h"
0007 #include "cln/ring.h"
0008 #include "cln/integer.h"
0009 #include "cln/random.h"
0010 #include "cln/malloc.h"
0011 #include "cln/io.h"
0012 #include "cln/proplist.h"
0013 #include "cln/condition.h"
0014 #include "cln/exception.h"
0015 #undef random // Linux defines random() as a macro!
0016 
0017 namespace cln {
0018 
0019 // Representation of an element of a ring Z/mZ.
0020 
0021 // To protect against mixing elements of different modular rings, such as
0022 // (3 mod 4) + (2 mod 5), every modular integer carries its ring in itself.
0023 
0024 
0025 // Representation of a ring Z/mZ.
0026 
0027 class cl_heap_modint_ring;
0028 
0029 class cl_modint_ring : public cl_ring {
0030 public:
0031     // Default constructor.
0032     cl_modint_ring ();
0033     // Constructor. Takes a cl_heap_modint_ring*, increments its refcount.
0034     cl_modint_ring (cl_heap_modint_ring* r);
0035     // Copy constructor.
0036     cl_modint_ring (const cl_modint_ring&);
0037     // Assignment operator.
0038     cl_modint_ring& operator= (const cl_modint_ring&);
0039     // Automatic dereferencing.
0040     cl_heap_modint_ring* operator-> () const
0041     { return (cl_heap_modint_ring*)heappointer; }
0042 };
0043 
0044 // Z/0Z
0045 extern const cl_modint_ring cl_modint0_ring;
0046 // Default constructor. This avoids dealing with NULL pointers.
0047 inline cl_modint_ring::cl_modint_ring ()
0048     : cl_ring (as_cl_private_thing(cl_modint0_ring)) {}
0049 
0050 class cl_MI_init_helper
0051 {
0052     static int count;
0053 public:
0054     cl_MI_init_helper();
0055     ~cl_MI_init_helper();
0056 };
0057 static cl_MI_init_helper cl_MI_init_helper_instance;
0058 
0059 // Copy constructor and assignment operator.
0060 CL_DEFINE_COPY_CONSTRUCTOR2(cl_modint_ring,cl_ring)
0061 CL_DEFINE_ASSIGNMENT_OPERATOR(cl_modint_ring,cl_modint_ring)
0062 
0063 // Normal constructor for `cl_modint_ring'.
0064 inline cl_modint_ring::cl_modint_ring (cl_heap_modint_ring* r)
0065     : cl_ring ((cl_private_thing) (cl_inc_pointer_refcount((cl_heap*)r), r)) {}
0066 
0067 // Operations on modular integer rings.
0068 
0069 inline bool operator== (const cl_modint_ring& R1, const cl_modint_ring& R2)
0070 { return (R1.pointer == R2.pointer); }
0071 inline bool operator!= (const cl_modint_ring& R1, const cl_modint_ring& R2)
0072 { return (R1.pointer != R2.pointer); }
0073 inline bool operator== (const cl_modint_ring& R1, cl_heap_modint_ring* R2)
0074 { return (R1.pointer == R2); }
0075 inline bool operator!= (const cl_modint_ring& R1, cl_heap_modint_ring* R2)
0076 { return (R1.pointer != R2); }
0077 
0078 
0079 // Condition raised when a probable prime is discovered to be composite.
0080 struct cl_composite_condition : public cl_condition {
0081     SUBCLASS_cl_condition()
0082     cl_I p; // the non-prime
0083     cl_I factor; // a nontrivial factor, or 0
0084     // Constructors.
0085     cl_composite_condition (const cl_I& _p)
0086         : p (_p), factor (0)
0087         { print(std::cerr); }
0088     cl_composite_condition (const cl_I& _p, const cl_I& _f)
0089         : p (_p), factor (_f)
0090         { print(std::cerr); }
0091     // Implement general condition methods.
0092     const char * name () const;
0093     void print (std::ostream&) const;
0094     ~cl_composite_condition () {}
0095 };
0096 
0097 
0098 // Representation of an element of a ring Z/mZ.
0099 
0100 class _cl_MI /* cf. _cl_ring_element */ {
0101 public:
0102     cl_I rep;       // representative, integer >=0, <m
0103                 // (maybe the Montgomery representative!)
0104     // Default constructor.
0105     _cl_MI () : rep () {}
0106 public: /* ugh */
0107     // Constructor.
0108     _cl_MI (const cl_heap_modint_ring* R, const cl_I& r) : rep (r) { (void)R; }
0109     _cl_MI (const cl_modint_ring& R, const cl_I& r) : rep (r) { (void)R; }
0110 public:
0111     // Conversion.
0112     CL_DEFINE_CONVERTER(_cl_ring_element)
0113 public: // Ability to place an object at a given address.
0114     void* operator new (size_t size) { return malloc_hook(size); }
0115     void* operator new (size_t size, void* ptr) { (void)size; return ptr; }
0116     void operator delete (void* ptr) { free_hook(ptr); }
0117 };
0118 
0119 class cl_MI /* cf. cl_ring_element */ : public _cl_MI {
0120 protected:
0121     cl_modint_ring _ring;   // ring Z/mZ
0122 public:
0123     const cl_modint_ring& ring () const { return _ring; }
0124     // Default constructor.
0125     cl_MI () : _cl_MI (), _ring () {}
0126 public: /* ugh */
0127     // Constructor.
0128     cl_MI (const cl_modint_ring& R, const cl_I& r) : _cl_MI (R,r), _ring (R) {}
0129     cl_MI (const cl_modint_ring& R, const _cl_MI& r) : _cl_MI (r), _ring (R) {}
0130 public:
0131     // Conversion.
0132     CL_DEFINE_CONVERTER(cl_ring_element)
0133     // Debugging output.
0134     void debug_print () const;
0135 public: // Ability to place an object at a given address.
0136     void* operator new (size_t size) { return malloc_hook(size); }
0137     void* operator new (size_t size, void* ptr) { (void)size; return ptr; }
0138     void operator delete (void* ptr) { free_hook(ptr); }
0139 };
0140 
0141 
0142 // Representation of an element of a ring Z/mZ or an exception.
0143 
0144 class cl_MI_x {
0145 private:
0146     cl_MI value;
0147 public:
0148     cl_composite_condition* condition;
0149     // Constructors.
0150     cl_MI_x (cl_composite_condition* c) : value (), condition (c) {}
0151     cl_MI_x (const cl_MI& x) : value (x), condition (NULL) {}
0152     // Cast operators.
0153     //operator cl_MI& () { if (condition) throw runtime_exception(); return value; }
0154     //operator const cl_MI& () const { if (condition) throw runtime_exception(); return value; }
0155     operator cl_MI () const { if (condition) throw runtime_exception(); return value; }
0156 };
0157 
0158 
0159 // Ring operations.
0160 
0161 struct _cl_modint_setops /* cf. _cl_ring_setops */ {
0162     // print
0163     void (* fprint) (cl_heap_modint_ring* R, std::ostream& stream, const _cl_MI& x);
0164     // equality
0165     bool (* equal) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
0166     // random number
0167     const _cl_MI (* random) (cl_heap_modint_ring* R, random_state& randomstate);
0168 };
0169 struct _cl_modint_addops /* cf. _cl_ring_addops */ {
0170     // 0
0171     const _cl_MI (* zero) (cl_heap_modint_ring* R);
0172     bool (* zerop) (cl_heap_modint_ring* R, const _cl_MI& x);
0173     // x+y
0174     const _cl_MI (* plus) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
0175     // x-y
0176     const _cl_MI (* minus) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
0177     // -x
0178     const _cl_MI (* uminus) (cl_heap_modint_ring* R, const _cl_MI& x);
0179 };
0180 struct _cl_modint_mulops /* cf. _cl_ring_mulops */ {
0181     // 1
0182     const _cl_MI (* one) (cl_heap_modint_ring* R);
0183     // canonical homomorphism
0184     const _cl_MI (* canonhom) (cl_heap_modint_ring* R, const cl_I& x);
0185     // x*y
0186     const _cl_MI (* mul) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
0187     // x^2
0188     const _cl_MI (* square) (cl_heap_modint_ring* R, const _cl_MI& x);
0189     // x^y, y Integer >0
0190     const _cl_MI (* expt_pos) (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y);
0191     // x^-1
0192     const cl_MI_x (* recip) (cl_heap_modint_ring* R, const _cl_MI& x);
0193     // x*y^-1
0194     const cl_MI_x (* div) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
0195     // x^y, y Integer
0196     const cl_MI_x (* expt) (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y);
0197     // x -> x mod m for x>=0
0198     const cl_I (* reduce_modulo) (cl_heap_modint_ring* R, const cl_I& x);
0199     // some inverse of canonical homomorphism
0200     const cl_I (* retract) (cl_heap_modint_ring* R, const _cl_MI& x);
0201 };
0202   typedef const _cl_modint_setops  cl_modint_setops;
0203   typedef const _cl_modint_addops  cl_modint_addops;
0204   typedef const _cl_modint_mulops  cl_modint_mulops;
0205 
0206 // Representation of the ring Z/mZ.
0207 
0208 // Currently rings are garbage collected only when they are not referenced
0209 // any more and when the ring table gets full.
0210 
0211 // Modular integer rings are kept unique in memory. This way, ring equality
0212 // can be checked very efficiently by a simple pointer comparison.
0213 
0214 class cl_heap_modint_ring /* cf. cl_heap_ring */ : public cl_heap {
0215     SUBCLASS_cl_heap_ring()
0216 private:
0217     cl_property_list properties;
0218 protected:
0219     cl_modint_setops* setops;
0220     cl_modint_addops* addops;
0221     cl_modint_mulops* mulops;
0222 public:
0223     cl_I modulus;   // m, normalized to be >= 0
0224 public:
0225     // Low-level operations.
0226     void _fprint (std::ostream& stream, const _cl_MI& x)
0227         { setops->fprint(this,stream,x); }
0228     bool _equal (const _cl_MI& x, const _cl_MI& y)
0229         { return setops->equal(this,x,y); }
0230     const _cl_MI _random (random_state& randomstate)
0231         { return setops->random(this,randomstate); }
0232     const _cl_MI _zero ()
0233         { return addops->zero(this); }
0234     bool _zerop (const _cl_MI& x)
0235         { return addops->zerop(this,x); }
0236     const _cl_MI _plus (const _cl_MI& x, const _cl_MI& y)
0237         { return addops->plus(this,x,y); }
0238     const _cl_MI _minus (const _cl_MI& x, const _cl_MI& y)
0239         { return addops->minus(this,x,y); }
0240     const _cl_MI _uminus (const _cl_MI& x)
0241         { return addops->uminus(this,x); }
0242     const _cl_MI _one ()
0243         { return mulops->one(this); }
0244     const _cl_MI _canonhom (const cl_I& x)
0245         { return mulops->canonhom(this,x); }
0246     const _cl_MI _mul (const _cl_MI& x, const _cl_MI& y)
0247         { return mulops->mul(this,x,y); }
0248     const _cl_MI _square (const _cl_MI& x)
0249         { return mulops->square(this,x); }
0250     const _cl_MI _expt_pos (const _cl_MI& x, const cl_I& y)
0251         { return mulops->expt_pos(this,x,y); }
0252     const cl_MI_x _recip (const _cl_MI& x)
0253         { return mulops->recip(this,x); }
0254     const cl_MI_x _div (const _cl_MI& x, const _cl_MI& y)
0255         { return mulops->div(this,x,y); }
0256     const cl_MI_x _expt (const _cl_MI& x, const cl_I& y)
0257         { return mulops->expt(this,x,y); }
0258     const cl_I _reduce_modulo (const cl_I& x)
0259         { return mulops->reduce_modulo(this,x); }
0260     const cl_I _retract (const _cl_MI& x)
0261         { return mulops->retract(this,x); }
0262     // High-level operations.
0263     void fprint (std::ostream& stream, const cl_MI& x)
0264     {
0265         if (!(x.ring() == this)) throw runtime_exception();
0266         _fprint(stream,x);
0267     }
0268     bool equal (const cl_MI& x, const cl_MI& y)
0269     {
0270         if (!(x.ring() == this)) throw runtime_exception();
0271         if (!(y.ring() == this)) throw runtime_exception();
0272         return _equal(x,y);
0273     }
0274     const cl_MI random (random_state& randomstate = default_random_state)
0275     {
0276         return cl_MI(this,_random(randomstate));
0277     }
0278     const cl_MI zero ()
0279     {
0280         return cl_MI(this,_zero());
0281     }
0282     bool zerop (const cl_MI& x)
0283     {
0284         if (!(x.ring() == this)) throw runtime_exception();
0285         return _zerop(x);
0286     }
0287     const cl_MI plus (const cl_MI& x, const cl_MI& y)
0288     {
0289         if (!(x.ring() == this)) throw runtime_exception();
0290         if (!(y.ring() == this)) throw runtime_exception();
0291         return cl_MI(this,_plus(x,y));
0292     }
0293     const cl_MI minus (const cl_MI& x, const cl_MI& y)
0294     {
0295         if (!(x.ring() == this)) throw runtime_exception();
0296         if (!(y.ring() == this)) throw runtime_exception();
0297         return cl_MI(this,_minus(x,y));
0298     }
0299     const cl_MI uminus (const cl_MI& x)
0300     {
0301         if (!(x.ring() == this)) throw runtime_exception();
0302         return cl_MI(this,_uminus(x));
0303     }
0304     const cl_MI one ()
0305     {
0306         return cl_MI(this,_one());
0307     }
0308     const cl_MI canonhom (const cl_I& x)
0309     {
0310         return cl_MI(this,_canonhom(x));
0311     }
0312     const cl_MI mul (const cl_MI& x, const cl_MI& y)
0313     {
0314         if (!(x.ring() == this)) throw runtime_exception();
0315         if (!(y.ring() == this)) throw runtime_exception();
0316         return cl_MI(this,_mul(x,y));
0317     }
0318     const cl_MI square (const cl_MI& x)
0319     {
0320         if (!(x.ring() == this)) throw runtime_exception();
0321         return cl_MI(this,_square(x));
0322     }
0323     const cl_MI expt_pos (const cl_MI& x, const cl_I& y)
0324     {
0325         if (!(x.ring() == this)) throw runtime_exception();
0326         return cl_MI(this,_expt_pos(x,y));
0327     }
0328     const cl_MI_x recip (const cl_MI& x)
0329     {
0330         if (!(x.ring() == this)) throw runtime_exception();
0331         return _recip(x);
0332     }
0333     const cl_MI_x div (const cl_MI& x, const cl_MI& y)
0334     {
0335         if (!(x.ring() == this)) throw runtime_exception();
0336         if (!(y.ring() == this)) throw runtime_exception();
0337         return _div(x,y);
0338     }
0339     const cl_MI_x expt (const cl_MI& x, const cl_I& y)
0340     {
0341         if (!(x.ring() == this)) throw runtime_exception();
0342         return _expt(x,y);
0343     }
0344     const cl_I reduce_modulo (const cl_I& x)
0345     {
0346         return _reduce_modulo(x);
0347     }
0348     const cl_I retract (const cl_MI& x)
0349     {
0350         if (!(x.ring() == this)) throw runtime_exception();
0351         return _retract(x);
0352     }
0353     // Miscellaneous.
0354     sintC bits; // number of bits needed to represent a representative, or -1
0355     int log2_bits; // log_2(bits), or -1
0356     // Property operations.
0357     cl_property* get_property (const cl_symbol& key)
0358         { return properties.get_property(key); }
0359     void add_property (cl_property* new_property)
0360         { properties.add_property(new_property); }
0361 // Constructor / destructor.
0362     cl_heap_modint_ring (cl_I m, cl_modint_setops*, cl_modint_addops*, cl_modint_mulops*);
0363     ~cl_heap_modint_ring () {}
0364 };
0365 #define SUBCLASS_cl_heap_modint_ring() \
0366   SUBCLASS_cl_heap_ring()
0367 
0368 // Lookup or create a modular integer ring  Z/mZ
0369 extern const cl_modint_ring find_modint_ring (const cl_I& m);
0370 static cl_MI_init_helper cl_MI_init_helper_instance2;
0371 
0372 // Operations on modular integers.
0373 
0374 // Output.
0375 inline void fprint (std::ostream& stream, const cl_MI& x)
0376     { x.ring()->fprint(stream,x); }
0377 CL_DEFINE_PRINT_OPERATOR(cl_MI)
0378 
0379 // Add.
0380 inline const cl_MI operator+ (const cl_MI& x, const cl_MI& y)
0381     { return x.ring()->plus(x,y); }
0382 inline const cl_MI operator+ (const cl_MI& x, const cl_I& y)
0383     { return x.ring()->plus(x,x.ring()->canonhom(y)); }
0384 inline const cl_MI operator+ (const cl_I& x, const cl_MI& y)
0385     { return y.ring()->plus(y.ring()->canonhom(x),y); }
0386 
0387 // Negate.
0388 inline const cl_MI operator- (const cl_MI& x)
0389     { return x.ring()->uminus(x); }
0390 
0391 // Subtract.
0392 inline const cl_MI operator- (const cl_MI& x, const cl_MI& y)
0393     { return x.ring()->minus(x,y); }
0394 inline const cl_MI operator- (const cl_MI& x, const cl_I& y)
0395     { return x.ring()->minus(x,x.ring()->canonhom(y)); }
0396 inline const cl_MI operator- (const cl_I& x, const cl_MI& y)
0397     { return y.ring()->minus(y.ring()->canonhom(x),y); }
0398 
0399 // Shifts.
0400 extern const cl_MI operator<< (const cl_MI& x, sintC y); // assume 0 <= y < 2^(intCsize-1)
0401 extern const cl_MI operator>> (const cl_MI& x, sintC y); // assume m odd, 0 <= y < 2^(intCsize-1)
0402 
0403 // Equality.
0404 inline bool operator== (const cl_MI& x, const cl_MI& y)
0405     { return x.ring()->equal(x,y); }
0406 inline bool operator!= (const cl_MI& x, const cl_MI& y)
0407     { return !x.ring()->equal(x,y); }
0408 inline bool operator== (const cl_MI& x, const cl_I& y)
0409     { return x.ring()->equal(x,x.ring()->canonhom(y)); }
0410 inline bool operator!= (const cl_MI& x, const cl_I& y)
0411     { return !x.ring()->equal(x,x.ring()->canonhom(y)); }
0412 inline bool operator== (const cl_I& x, const cl_MI& y)
0413     { return y.ring()->equal(y.ring()->canonhom(x),y); }
0414 inline bool operator!= (const cl_I& x, const cl_MI& y)
0415     { return !y.ring()->equal(y.ring()->canonhom(x),y); }
0416 
0417 // Compare against 0.
0418 inline bool zerop (const cl_MI& x)
0419     { return x.ring()->zerop(x); }
0420 
0421 // Multiply.
0422 inline const cl_MI operator* (const cl_MI& x, const cl_MI& y)
0423     { return x.ring()->mul(x,y); }
0424 
0425 // Squaring.
0426 inline const cl_MI square (const cl_MI& x)
0427     { return x.ring()->square(x); }
0428 
0429 // Exponentiation x^y, where y > 0.
0430 inline const cl_MI expt_pos (const cl_MI& x, const cl_I& y)
0431     { return x.ring()->expt_pos(x,y); }
0432 
0433 // Reciprocal.
0434 inline const cl_MI recip (const cl_MI& x)
0435     { return x.ring()->recip(x); }
0436 
0437 // Division.
0438 inline const cl_MI div (const cl_MI& x, const cl_MI& y)
0439     { return x.ring()->div(x,y); }
0440 inline const cl_MI div (const cl_MI& x, const cl_I& y)
0441     { return x.ring()->div(x,x.ring()->canonhom(y)); }
0442 inline const cl_MI div (const cl_I& x, const cl_MI& y)
0443     { return y.ring()->div(y.ring()->canonhom(x),y); }
0444 
0445 // Exponentiation x^y.
0446 inline const cl_MI expt (const cl_MI& x, const cl_I& y)
0447     { return x.ring()->expt(x,y); }
0448 
0449 // Scalar multiplication.
0450 inline const cl_MI operator* (const cl_I& x, const cl_MI& y)
0451     { return y.ring()->mul(y.ring()->canonhom(x),y); }
0452 inline const cl_MI operator* (const cl_MI& x, const cl_I& y)
0453     { return x.ring()->mul(x.ring()->canonhom(y),x); }
0454 
0455 // TODO: implement gcd, index (= gcd), unitp, sqrtp
0456 
0457 
0458 // Debugging support.
0459 #ifdef CL_DEBUG
0460 extern int cl_MI_debug_module;
0461 CL_FORCE_LINK(cl_MI_debug_dummy, cl_MI_debug_module)
0462 #endif
0463 
0464 }  // namespace cln
0465 
0466 #endif /* _CL_MODINTEGER_H */