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