include/cln/modinteger.h

   1 // Modular integer operations.
   2
   3 #ifndef _CL_MODINTEGER_H
   4 #define _CL_MODINTEGER_H
   5
   6 #include "cln/object.h"
   7 #include "cln/ring.h"
   8 #include "cln/integer.h"
   9 #include "cln/random.h"
  10 #include "cln/malloc.h"
  11 #include "cln/io.h"
  12 #include "cln/proplist.h"
  13 #include "cln/condition.h"
  14 #include "cln/exception.h"
  15 #undef random // Linux defines random() as a macro!
  16
  17 namespace cln {
  18
  19 // Representation of an element of a ring Z/mZ.
  20
  21 // To protect against mixing elements of different modular rings, such as
  22 // (3 mod 4) + (2 mod 5), every modular integer carries its ring in itself.
  23
  24
  25 // Representation of a ring Z/mZ.
  26
  27 class cl_heap_modint_ring;
  28
  29 class cl_modint_ring : public cl_ring {
  30 public:
  31         // Default constructor.
  32         cl_modint_ring ();
  33         // Constructor. Takes a cl_heap_modint_ring*, increments its refcount.
  34         cl_modint_ring (cl_heap_modint_ring* r);
  35         // Copy constructor.
  36         cl_modint_ring (const cl_modint_ring&);
  37         // Assignment operator.
  38         cl_modint_ring& operator= (const cl_modint_ring&);
  39         // Automatic dereferencing.
  40         cl_heap_modint_ring* operator-> () const
  41         { return (cl_heap_modint_ring*)heappointer; }
  42 };
  43
  44 // Z/0Z
  45 extern const cl_modint_ring cl_modint0_ring;
  46 // Default constructor. This avoids dealing with NULL pointers.
  47 inline cl_modint_ring::cl_modint_ring ()
  48         : cl_ring (as_cl_private_thing(cl_modint0_ring)) {}
  49
  50 class cl_MI_init_helper
  51 {
  52         static int count;
  53 public:
  54         cl_MI_init_helper();
  55         ~cl_MI_init_helper();
  56 };
  57 static cl_MI_init_helper cl_MI_init_helper_instance;
  58
  59 // Copy constructor and assignment operator.
  60 CL_DEFINE_COPY_CONSTRUCTOR2(cl_modint_ring,cl_ring)
  61 CL_DEFINE_ASSIGNMENT_OPERATOR(cl_modint_ring,cl_modint_ring)
  62
  63 // Normal constructor for `cl_modint_ring'.
  64 inline cl_modint_ring::cl_modint_ring (cl_heap_modint_ring* r)
  65         : cl_ring ((cl_private_thing) (cl_inc_pointer_refcount((cl_heap*)r), r)) {}
  66
  67 // Operations on modular integer rings.
  68
  69 inline bool operator== (const cl_modint_ring& R1, const cl_modint_ring& R2)
  70 { return (R1.pointer == R2.pointer); }
  71 inline bool operator!= (const cl_modint_ring& R1, const cl_modint_ring& R2)
  72 { return (R1.pointer != R2.pointer); }
  73 inline bool operator== (const cl_modint_ring& R1, cl_heap_modint_ring* R2)
  74 { return (R1.pointer == R2); }
  75 inline bool operator!= (const cl_modint_ring& R1, cl_heap_modint_ring* R2)
  76 { return (R1.pointer != R2); }
  77
  78
  79 // Condition raised when a probable prime is discovered to be composite.
  80 struct cl_composite_condition : public cl_condition {
  81         SUBCLASS_cl_condition()
  82         cl_I p; // the non-prime
  83         cl_I factor; // a nontrivial factor, or 0
  84         // Constructors.
  85         cl_composite_condition (const cl_I& _p)
  86                 : p (_p), factor (0)
  87                 { print(std::cerr); }
  88         cl_composite_condition (const cl_I& _p, const cl_I& _f)
  89                 : p (_p), factor (_f)
  90                 { print(std::cerr); }
  91         // Implement general condition methods.
  92         const char * name () const;
  93         void print (std::ostream&) const;
  94         ~cl_composite_condition () {}
  95 };
  96
  97
  98 // Representation of an element of a ring Z/mZ.
  99
 100 class _cl_MI /* cf. _cl_ring_element */ {
 101 public:
 102         cl_I rep;               // representative, integer >=0, <m
 103                                 // (maybe the Montgomery representative!)
 104         // Default constructor.
 105         _cl_MI () : rep () {}
 106 public: /* ugh */
 107         // Constructor.
 108         _cl_MI (const cl_heap_modint_ring* R, const cl_I& r) : rep (r) { (void)R; }
 109         _cl_MI (const cl_modint_ring& R, const cl_I& r) : rep (r) { (void)R; }
 110 public:
 111         // Conversion.
 112         CL_DEFINE_CONVERTER(_cl_ring_element)
 113 public: // Ability to place an object at a given address.
 114         void* operator new (size_t size) { return malloc_hook(size); }
 115         void* operator new (size_t size, void* ptr) { (void)size; return ptr; }
 116         void operator delete (void* ptr) { free_hook(ptr); }
 117 };
 118
 119 class cl_MI /* cf. cl_ring_element */ : public _cl_MI {
 120 protected:
 121         cl_modint_ring _ring;   // ring Z/mZ
 122 public:
 123         const cl_modint_ring& ring () const { return _ring; }
 124         // Default constructor.
 125         cl_MI () : _cl_MI (), _ring () {}
 126 public: /* ugh */
 127         // Constructor.
 128         cl_MI (const cl_modint_ring& R, const cl_I& r) : _cl_MI (R,r), _ring (R) {}
 129         cl_MI (const cl_modint_ring& R, const _cl_MI& r) : _cl_MI (r), _ring (R) {}
 130 public:
 131         // Conversion.
 132         CL_DEFINE_CONVERTER(cl_ring_element)
 133         // Debugging output.
 134         void debug_print () const;
 135 public: // Ability to place an object at a given address.
 136         void* operator new (size_t size) { return malloc_hook(size); }
 137         void* operator new (size_t size, void* ptr) { (void)size; return ptr; }
 138         void operator delete (void* ptr) { free_hook(ptr); }
 139 };
 140
 141
 142 // Representation of an element of a ring Z/mZ or an exception.
 143
 144 class cl_MI_x {
 145 private:
 146         cl_MI value;
 147 public:
 148         cl_composite_condition* condition;
 149         // Constructors.
 150         cl_MI_x (cl_composite_condition* c) : value (), condition (c) {}
 151         cl_MI_x (const cl_MI& x) : value (x), condition (NULL) {}
 152         // Cast operators.
 153         //operator cl_MI& () { if (condition) throw runtime_exception(); return value; }
 154         //operator const cl_MI& () const { if (condition) throw runtime_exception(); return value; }
 155         operator cl_MI () const { if (condition) throw runtime_exception(); return value; }
 156 };
 157
 158
 159 // Ring operations.
 160
 161 struct _cl_modint_setops /* cf. _cl_ring_setops */ {
 162         // print
 163         void (* fprint) (cl_heap_modint_ring* R, std::ostream& stream, const _cl_MI& x);
 164         // equality
 165         bool (* equal) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
 166         // random number
 167         const _cl_MI (* random) (cl_heap_modint_ring* R, random_state& randomstate);
 168 };
 169 struct _cl_modint_addops /* cf. _cl_ring_addops */ {
 170         // 0
 171         const _cl_MI (* zero) (cl_heap_modint_ring* R);
 172         bool (* zerop) (cl_heap_modint_ring* R, const _cl_MI& x);
 173         // x+y
 174         const _cl_MI (* plus) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
 175         // x-y
 176         const _cl_MI (* minus) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
 177         // -x
 178         const _cl_MI (* uminus) (cl_heap_modint_ring* R, const _cl_MI& x);
 179 };
 180 struct _cl_modint_mulops /* cf. _cl_ring_mulops */ {
 181         // 1
 182         const _cl_MI (* one) (cl_heap_modint_ring* R);
 183         // canonical homomorphism
 184         const _cl_MI (* canonhom) (cl_heap_modint_ring* R, const cl_I& x);
 185         // x*y
 186         const _cl_MI (* mul) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
 187         // x^2
 188         const _cl_MI (* square) (cl_heap_modint_ring* R, const _cl_MI& x);
 189         // x^y, y Integer >0
 190         const _cl_MI (* expt_pos) (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y);
 191         // x^-1
 192         const cl_MI_x (* recip) (cl_heap_modint_ring* R, const _cl_MI& x);
 193         // x*y^-1
 194         const cl_MI_x (* div) (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y);
 195         // x^y, y Integer
 196         const cl_MI_x (* expt) (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y);
 197         // x -> x mod m for x>=0
 198         const cl_I (* reduce_modulo) (cl_heap_modint_ring* R, const cl_I& x);
 199         // some inverse of canonical homomorphism
 200         const cl_I (* retract) (cl_heap_modint_ring* R, const _cl_MI& x);
 201 };
 202   typedef const _cl_modint_setops  cl_modint_setops;
 203   typedef const _cl_modint_addops  cl_modint_addops;
 204   typedef const _cl_modint_mulops  cl_modint_mulops;
 205
 206 // Representation of the ring Z/mZ.
 207
 208 // Currently rings are garbage collected only when they are not referenced
 209 // any more and when the ring table gets full.
 210
 211 // Modular integer rings are kept unique in memory. This way, ring equality
 212 // can be checked very efficiently by a simple pointer comparison.
 213
 214 class cl_heap_modint_ring /* cf. cl_heap_ring */ : public cl_heap {
 215         SUBCLASS_cl_heap_ring()
 216 private:
 217         cl_property_list properties;
 218 protected:
 219         cl_modint_setops* setops;
 220         cl_modint_addops* addops;
 221         cl_modint_mulops* mulops;
 222 public:
 223         cl_I modulus;   // m, normalized to be >= 0
 224 public:
 225         // Low-level operations.
 226         void _fprint (std::ostream& stream, const _cl_MI& x)
 227                 { setops->fprint(this,stream,x); }
 228         bool _equal (const _cl_MI& x, const _cl_MI& y)
 229                 { return setops->equal(this,x,y); }
 230         const _cl_MI _random (random_state& randomstate)
 231                 { return setops->random(this,randomstate); }
 232         const _cl_MI _zero ()
 233                 { return addops->zero(this); }
 234         bool _zerop (const _cl_MI& x)
 235                 { return addops->zerop(this,x); }
 236         const _cl_MI _plus (const _cl_MI& x, const _cl_MI& y)
 237                 { return addops->plus(this,x,y); }
 238         const _cl_MI _minus (const _cl_MI& x, const _cl_MI& y)
 239                 { return addops->minus(this,x,y); }
 240         const _cl_MI _uminus (const _cl_MI& x)
 241                 { return addops->uminus(this,x); }
 242         const _cl_MI _one ()
 243                 { return mulops->one(this); }
 244         const _cl_MI _canonhom (const cl_I& x)
 245                 { return mulops->canonhom(this,x); }
 246         const _cl_MI _mul (const _cl_MI& x, const _cl_MI& y)
 247                 { return mulops->mul(this,x,y); }
 248         const _cl_MI _square (const _cl_MI& x)
 249                 { return mulops->square(this,x); }
 250         const _cl_MI _expt_pos (const _cl_MI& x, const cl_I& y)
 251                 { return mulops->expt_pos(this,x,y); }
 252         const cl_MI_x _recip (const _cl_MI& x)
 253                 { return mulops->recip(this,x); }
 254         const cl_MI_x _div (const _cl_MI& x, const _cl_MI& y)
 255                 { return mulops->div(this,x,y); }
 256         const cl_MI_x _expt (const _cl_MI& x, const cl_I& y)
 257                 { return mulops->expt(this,x,y); }
 258         const cl_I _reduce_modulo (const cl_I& x)
 259                 { return mulops->reduce_modulo(this,x); }
 260         const cl_I _retract (const _cl_MI& x)
 261                 { return mulops->retract(this,x); }
 262         // High-level operations.
 263         void fprint (std::ostream& stream, const cl_MI& x)
 264         {
 265                 if (!(x.ring() == this)) throw runtime_exception();
 266                 _fprint(stream,x);
 267         }
 268         bool equal (const cl_MI& x, const cl_MI& y)
 269         {
 270                 if (!(x.ring() == this)) throw runtime_exception();
 271                 if (!(y.ring() == this)) throw runtime_exception();
 272                 return _equal(x,y);
 273         }
 274         const cl_MI random (random_state& randomstate = default_random_state)
 275         {
 276                 return cl_MI(this,_random(randomstate));
 277         }
 278         const cl_MI zero ()
 279         {
 280                 return cl_MI(this,_zero());
 281         }
 282         bool zerop (const cl_MI& x)
 283         {
 284                 if (!(x.ring() == this)) throw runtime_exception();
 285                 return _zerop(x);
 286         }
 287         const cl_MI plus (const cl_MI& x, const cl_MI& y)
 288         {
 289                 if (!(x.ring() == this)) throw runtime_exception();
 290                 if (!(y.ring() == this)) throw runtime_exception();
 291                 return cl_MI(this,_plus(x,y));
 292         }
 293         const cl_MI minus (const cl_MI& x, const cl_MI& y)
 294         {
 295                 if (!(x.ring() == this)) throw runtime_exception();
 296                 if (!(y.ring() == this)) throw runtime_exception();
 297                 return cl_MI(this,_minus(x,y));
 298         }
 299         const cl_MI uminus (const cl_MI& x)
 300         {
 301                 if (!(x.ring() == this)) throw runtime_exception();
 302                 return cl_MI(this,_uminus(x));
 303         }
 304         const cl_MI one ()
 305         {
 306                 return cl_MI(this,_one());
 307         }
 308         const cl_MI canonhom (const cl_I& x)
 309         {
 310                 return cl_MI(this,_canonhom(x));
 311         }
 312         const cl_MI mul (const cl_MI& x, const cl_MI& y)
 313         {
 314                 if (!(x.ring() == this)) throw runtime_exception();
 315                 if (!(y.ring() == this)) throw runtime_exception();
 316                 return cl_MI(this,_mul(x,y));
 317         }
 318         const cl_MI square (const cl_MI& x)
 319         {
 320                 if (!(x.ring() == this)) throw runtime_exception();
 321                 return cl_MI(this,_square(x));
 322         }
 323         const cl_MI expt_pos (const cl_MI& x, const cl_I& y)
 324         {
 325                 if (!(x.ring() == this)) throw runtime_exception();
 326                 return cl_MI(this,_expt_pos(x,y));
 327         }
 328         const cl_MI_x recip (const cl_MI& x)
 329         {
 330                 if (!(x.ring() == this)) throw runtime_exception();
 331                 return _recip(x);
 332         }
 333         const cl_MI_x div (const cl_MI& x, const cl_MI& y)
 334         {
 335                 if (!(x.ring() == this)) throw runtime_exception();
 336                 if (!(y.ring() == this)) throw runtime_exception();
 337                 return _div(x,y);
 338         }
 339         const cl_MI_x expt (const cl_MI& x, const cl_I& y)
 340         {
 341                 if (!(x.ring() == this)) throw runtime_exception();
 342                 return _expt(x,y);
 343         }
 344         const cl_I reduce_modulo (const cl_I& x)
 345         {
 346                 return _reduce_modulo(x);
 347         }
 348         const cl_I retract (const cl_MI& x)
 349         {
 350                 if (!(x.ring() == this)) throw runtime_exception();
 351                 return _retract(x);
 352         }
 353         // Miscellaneous.
 354         sintC bits; // number of bits needed to represent a representative, or -1
 355         int log2_bits; // log_2(bits), or -1
 356         // Property operations.
 357         cl_property* get_property (const cl_symbol& key)
 358                 { return properties.get_property(key); }
 359         void add_property (cl_property* new_property)
 360                 { properties.add_property(new_property); }
 361 // Constructor / destructor.
 362         cl_heap_modint_ring (cl_I m, cl_modint_setops*, cl_modint_addops*, cl_modint_mulops*);
 363         ~cl_heap_modint_ring () {}
 364 };
 365 #define SUBCLASS_cl_heap_modint_ring() \
 366   SUBCLASS_cl_heap_ring()
 367
 368 // Lookup or create a modular integer ring  Z/mZ
 369 extern const cl_modint_ring find_modint_ring (const cl_I& m);
 370 static cl_MI_init_helper cl_MI_init_helper_instance2;
 371
 372 // Operations on modular integers.
 373
 374 // Output.
 375 inline void fprint (std::ostream& stream, const cl_MI& x)
 376         { x.ring()->fprint(stream,x); }
 377 CL_DEFINE_PRINT_OPERATOR(cl_MI)
 378
 379 // Add.
 380 inline const cl_MI operator+ (const cl_MI& x, const cl_MI& y)
 381         { return x.ring()->plus(x,y); }
 382 inline const cl_MI operator+ (const cl_MI& x, const cl_I& y)
 383         { return x.ring()->plus(x,x.ring()->canonhom(y)); }
 384 inline const cl_MI operator+ (const cl_I& x, const cl_MI& y)
 385         { return y.ring()->plus(y.ring()->canonhom(x),y); }
 386
 387 // Negate.
 388 inline const cl_MI operator- (const cl_MI& x)
 389         { return x.ring()->uminus(x); }
 390
 391 // Subtract.
 392 inline const cl_MI operator- (const cl_MI& x, const cl_MI& y)
 393         { return x.ring()->minus(x,y); }
 394 inline const cl_MI operator- (const cl_MI& x, const cl_I& y)
 395         { return x.ring()->minus(x,x.ring()->canonhom(y)); }
 396 inline const cl_MI operator- (const cl_I& x, const cl_MI& y)
 397         { return y.ring()->minus(y.ring()->canonhom(x),y); }
 398
 399 // Shifts.
 400 extern const cl_MI operator<< (const cl_MI& x, sintC y); // assume 0 <= y < 2^(intCsize-1)
 401 extern const cl_MI operator>> (const cl_MI& x, sintC y); // assume m odd, 0 <= y < 2^(intCsize-1)
 402
 403 // Equality.
 404 inline bool operator== (const cl_MI& x, const cl_MI& y)
 405         { return x.ring()->equal(x,y); }
 406 inline bool operator!= (const cl_MI& x, const cl_MI& y)
 407         { return !x.ring()->equal(x,y); }
 408 inline bool operator== (const cl_MI& x, const cl_I& y)
 409         { return x.ring()->equal(x,x.ring()->canonhom(y)); }
 410 inline bool operator!= (const cl_MI& x, const cl_I& y)
 411         { return !x.ring()->equal(x,x.ring()->canonhom(y)); }
 412 inline bool operator== (const cl_I& x, const cl_MI& y)
 413         { return y.ring()->equal(y.ring()->canonhom(x),y); }
 414 inline bool operator!= (const cl_I& x, const cl_MI& y)
 415         { return !y.ring()->equal(y.ring()->canonhom(x),y); }
 416
 417 // Compare against 0.
 418 inline bool zerop (const cl_MI& x)
 419         { return x.ring()->zerop(x); }
 420
 421 // Multiply.
 422 inline const cl_MI operator* (const cl_MI& x, const cl_MI& y)
 423         { return x.ring()->mul(x,y); }
 424
 425 // Squaring.
 426 inline const cl_MI square (const cl_MI& x)
 427         { return x.ring()->square(x); }
 428
 429 // Exponentiation x^y, where y > 0.
 430 inline const cl_MI expt_pos (const cl_MI& x, const cl_I& y)
 431         { return x.ring()->expt_pos(x,y); }
 432
 433 // Reciprocal.
 434 inline const cl_MI recip (const cl_MI& x)
 435         { return x.ring()->recip(x); }
 436
 437 // Division.
 438 inline const cl_MI div (const cl_MI& x, const cl_MI& y)
 439         { return x.ring()->div(x,y); }
 440 inline const cl_MI div (const cl_MI& x, const cl_I& y)
 441         { return x.ring()->div(x,x.ring()->canonhom(y)); }
 442 inline const cl_MI div (const cl_I& x, const cl_MI& y)
 443         { return y.ring()->div(y.ring()->canonhom(x),y); }
 444
 445 // Exponentiation x^y.
 446 inline const cl_MI expt (const cl_MI& x, const cl_I& y)
 447         { return x.ring()->expt(x,y); }
 448
 449 // Scalar multiplication.
 450 inline const cl_MI operator* (const cl_I& x, const cl_MI& y)
 451         { return y.ring()->mul(y.ring()->canonhom(x),y); }
 452 inline const cl_MI operator* (const cl_MI& x, const cl_I& y)
 453         { return x.ring()->mul(x.ring()->canonhom(y),x); }
 454
 455 // TODO: implement gcd, index (= gcd), unitp, sqrtp
 456
 457
 458 // Debugging support.
 459 #ifdef CL_DEBUG
 460 extern int cl_MI_debug_module;
 461 CL_FORCE_LINK(cl_MI_debug_dummy, cl_MI_debug_module)
 462 #endif
 463
 464 }  // namespace cln
 465
 466 #endif /* _CL_MODINTEGER_H */