include/cl_modinteger.h

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