Package sage :: Package groups :: Module generic
[hide private]
[frames] | no frames]

Source Code for Module sage.groups.generic

   1  r""" 
   2  Miscellaneous generic functions 
   3   
   4  A collection of functions implementing generic algorithms in arbitrary 
   5  groups, including additive and multiplicative groups. 
   6   
   7  In all cases the group operation is specified by a parameter 
   8  'operation', which is a string either one of the set of 
   9  multiplication_names or addition_names specified below, or 'other'. 
  10  In the latter case, the caller must provide an identity, inverse() and 
  11  op() functions. 
  12   
  13  Also included are a generic function for computing multiples (or 
  14  powers), and an iterator for general multiples and powers. 
  15   
  16  EXAMPLES: 
  17   
  18  Some examples in the multiplicative group of a finite field: 
  19   
  20  Discrete logs: 
  21      sage: K = GF(3^6,'b') 
  22      sage: b = K.gen() 
  23      sage: a = b^210 
  24      sage: discrete_log(a, b, K.order()-1) 
  25      210 
  26   
  27  Linear relation finder: 
  28      sage: F.<a>=GF(3^6,'a') 
  29      sage: a.multiplicative_order().factor() 
  30      2^3 * 7 * 13 
  31      sage: b=a^7 
  32      sage: c=a^13 
  33      sage: linear_relation(b,c,'*') 
  34      (13, 7) 
  35      sage: b^13==c^7 
  36      True 
  37   
  38  Orders of elements: 
  39      sage: k.<a> = GF(5^5) 
  40      sage: b = a^4 
  41      sage: order_from_multiple(b,5^5-1,operation='*') 
  42      781 
  43      sage: order_from_bounds(b,(5^4,5^5),operation='*') 
  44      781 
  45   
  46  Some examples in the group of points of an elliptic curve over a finite field: 
  47   
  48  Discrete logs: 
  49      sage: F=GF(37^2,'a') 
  50      sage: E=EllipticCurve(F,[1,1]) 
  51      sage: F.<a>=GF(37^2,'a') 
  52      sage: E=EllipticCurve(F,[1,1]) 
  53      sage: P=E(25*a + 16 , 15*a + 7 ) 
  54      sage: P.order() 
  55      672 
  56      sage: Q=39*P; Q 
  57      (36*a + 32 : 5*a + 12 : 1) 
  58      sage: discrete_log(Q,P,P.order(),'+') 
  59      39 
  60   
  61  Linear relation finder: 
  62      sage: F.<a>=GF(3^6,'a') 
  63      sage: E=EllipticCurve([a^5 + 2*a^3 + 2*a^2 + 2*a, a^4 + a^3 + 2*a + 1]) 
  64      sage: P=E(a^5 + a^4 + a^3 + a^2 + a + 2 , 0) 
  65      sage: Q=E(2*a^3 + 2*a^2 + 2*a , a^3 + 2*a^2 + 1) 
  66      sage: linear_relation(P,Q,'+') 
  67      (1, 2) 
  68      sage: P == 2*Q 
  69      True 
  70   
  71  Orders of elements: 
  72      sage: k.<a> = GF(5^5) 
  73      sage: E = EllipticCurve(k,[2,4]) 
  74      sage: P = E(3*a^4 + 3*a , 2*a + 1 ) 
  75      sage: M = E.cardinality(); M 
  76      3227 
  77      sage: plist = M.prime_factors() 
  78      sage: order_from_multiple(P, M, plist, operation='+') 
  79      3227 
  80      sage: Q = E(0,2) 
  81      sage: order_from_multiple(Q, M, plist, operation='+') 
  82      7 
  83      sage: order_from_bounds(Q, Hasse_bounds(5^5), operation='+') 
  84      7 
  85   
  86  """ 
  87   
  88  ########################################################################### 
  89  #       Copyright (C) 2008 William Stein <wstein@gmail.com> 
  90  #                          John Cremona  <john.cremona@gmail.com>    
  91  # 
  92  #  Distributed under the terms of the GNU General Public License (GPL) 
  93  #                  http://www.gnu.org/licenses/ 
  94  ########################################################################### 
  95   
  96  from copy import copy 
  97  import math 
  98   
  99  import sage.misc.misc as misc 
 100  import sage.misc.search  
 101  import sage.rings.integer_ring as integer_ring 
 102  import sage.rings.integer 
 103   
 104  # 
 105  # Lists of names (as strings) which the user may use to identify one 
 106  # of the standard operations: 
 107  # 
 108  multiplication_names = ( 'multiplication', 'times', 'product', '*') 
 109  addition_names       = ( 'addition', 'plus', 'sum', '+') 
 110   
 111   
 112  # 
 113  # This function was moved from sage/rings/arith.py 
 114  # 
115 -def power(a, m, one=1):
116 """ 117 The m-th power of a, where m is a non-negative 118 integer and a is a Python object on which 119 multiplication is defined. The exponentiation 120 is done using the standard binary powering algorithm. 121 122 EXAMPLES: 123 sage: power(2,5) 124 32 125 sage: power(RealField()('2.5'),4) 126 39.0625000000000 127 sage: power(0,0) 128 Traceback (most recent call last): 129 ... 130 ArithmeticError: 0^0 is undefined. 131 sage: power(2,-3) 132 1/8 133 """ 134 return sage.structure.element.generic_power(a,m,one)
135
136 -def multiple(a, n, operation='*', identity=None, inverse=None, op=None):
137 """ 138 Returns n*a [or a**n], where n is any integer and a is a Python 139 object on which a group operaton such as addition or 140 multiplication is defined. Uses the standard binary algorithm. 141 142 EXAMPLES: 143 sage: multiple(2,5) 144 32 145 sage: multiple(RealField()('2.5'),4) 146 39.0625000000000 147 sage: multiple(2,-3) 148 1/8 149 sage: multiple(2,100,'+') == 100*2 150 True 151 sage: multiple(2,100) == 2**100 152 True 153 sage: multiple(2,-100,) == 2**-100 154 True 155 sage: R.<x>=ZZ[] 156 sage: multiple(x,100) 157 x^100 158 sage: multiple(x,100,'+') 159 100*x 160 sage: multiple(x,-10) 161 1/x^10 162 163 Idempotence is detected making this fast: 164 sage: multiple(1,10^1000) 165 1 166 167 sage: E=EllipticCurve('389a1') 168 sage: P=E(-1,1) 169 sage: multiple(P,10,'+') 170 (645656132358737542773209599489/22817025904944891235367494656 : 525532176124281192881231818644174845702936831/3446581505217248068297884384990762467229696 : 1) 171 sage: multiple(P,-10,'+') 172 (645656132358737542773209599489/22817025904944891235367494656 : -528978757629498440949529703029165608170166527/3446581505217248068297884384990762467229696 : 1) 173 174 175 """ 176 from operator import inv, mul, neg, add 177 178 if operation in multiplication_names: 179 identity = a.parent()(1) 180 inverse = inv 181 op = mul 182 elif operation in addition_names: 183 identity = a.parent()(0) 184 inverse = neg 185 op = add 186 else: 187 if identity==None or inverse==None or op==None: 188 raise ValueError, "identity, inverse and operation must all be specified" 189 190 if n == 0: 191 return identity 192 193 if n < 0: 194 n = -n 195 a = inverse(a) 196 197 if n == 1: 198 return a 199 200 # check for idempotence, and store the result otherwise 201 aa = op(a,a) 202 if aa == a: 203 return a 204 205 if n == 2: 206 return aa 207 208 if n == 3: 209 return op(aa,a) 210 211 if n == 4: 212 return op(aa,aa) 213 214 # since we've computed a^2, let's start squaring there 215 # so, let's keep the least-significant bit around, just 216 # in case. 217 m = n & 1 218 n = n >> 1 219 220 # One multiplication can be saved by starting with 221 # the second-smallest power needed rather than with 1 222 # we've already squared a, so let's start there. 223 apow = aa 224 while n&1 == 0: 225 apow = op(apow,apow) 226 n = n >> 1 227 power = apow 228 n = n >> 1 229 230 # now multiply that least-significant bit in... 231 if m: 232 power = op(power,a) 233 234 # and this is straight from the book. 235 while n != 0: 236 apow = op(apow,apow) 237 if n&1 != 0: 238 power = op(power,apow) 239 n = n >> 1 240 241 return power
242 243 244 # 245 # Generic iterator for looping through multiples or powers 246 # 247
248 -class multiples:
249 """ 250 Return an iterator which runs through P0+i*P for i in range(n) 251 252 P and P0 must be Sage objects in some group; if the operation in 253 multiplcation then the returned values are P0*P**i. 254 255 EXAMPLES 256 sage: list(multiples(1,10)) 257 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] 258 sage: list(multiples(1,10,100)) 259 [100, 101, 102, 103, 104, 105, 106, 107, 108, 109] 260 261 sage: E=EllipticCurve('389a1') 262 sage: P=E(-1,1) 263 sage: for Q in multiples(P,5): print Q, Q.height()/P.height() 264 (0 : 1 : 0) 0.000000000000000 265 (-1 : 1 : 1) 1.00000000000000 266 (10/9 : -35/27 : 1) 4.00000000000000 267 (26/361 : -5720/6859 : 1) 9.00000000000000 268 (47503/16641 : 9862190/2146689 : 1) 16.0000000000000 269 270 sage: R.<x>=ZZ[] 271 sage: list(multiples(x,5)) 272 [0, x, 2*x, 3*x, 4*x] 273 sage: list(multiples(x,5,operation='*')) 274 [1, x, x^2, x^3, x^4] 275 sage: list(multiples(x,5,indexed=True)) 276 [(0, 0), (1, x), (2, 2*x), (3, 3*x), (4, 4*x)] 277 sage: list(multiples(x,5,indexed=True,operation='*')) 278 [(0, 1), (1, x), (2, x^2), (3, x^3), (4, x^4)] 279 sage: for i,y in multiples(x,5,indexed=True): print "%s times %s = %s"%(i,x,y) 280 0 times x = 0 281 1 times x = x 282 2 times x = 2*x 283 3 times x = 3*x 284 4 times x = 4*x 285 286 sage: for i,n in multiples(3,5,indexed=True,operation='*'): print "3 to the power %s = %s"%(i,n) 287 3 to the power 0 = 1 288 3 to the power 1 = 3 289 3 to the power 2 = 9 290 3 to the power 3 = 27 291 3 to the power 4 = 81 292 """
293 - def __init__(self,P,n,P0=None,indexed=False, operation='+', op=None):
294 """ 295 Create a multiples iterator 296 297 INPUT: 298 299 P -- step value: any Sage object on which a binary 300 operation is defined 301 302 n -- number of multiples: non-negative integer 303 P0 -- offset (default 0): Sage object which can be 'added' to P 304 indexed -- boolean (default False) 305 306 operation -- string: '+' (default ) or 307 '*' or other. 308 309 If other, a function op() must be supplied (a 310 function of 2 arguments) defining the 311 group binary operation; also either P0 must be supplied. 312 313 If indexed==False the iterator delivers P0+i*P (if 314 operation=='+') or P0*P**i (if 315 operation=='*'), for i in range(n) . 316 317 If indexed==True the iterator delivers tuples (i,P0+i*P) 318 or (i,P0*P**i) 319 """ 320 if n<0: 321 raise ValueError, 'n cannot be negative in multiples' 322 323 from operator import mul, add 324 325 if operation in multiplication_names: 326 if P0==None: P0 = P.parent()(1) 327 self.op = mul 328 elif operation in addition_names: 329 if P0==None: P0 = P.parent()(0) 330 self.op = add 331 else: 332 self.op = op 333 if P0==None: 334 raise ValueError, "P0 must be supplied when operation is neither addition nor multiplication" 335 if op==None: 336 raise ValueError, "op() must both be supplied when operation is neither addition nor multiplication" 337 338 self.P=copy(P) 339 self.Q=copy(P0) 340 assert self.P != None and self.Q!=None 341 self.i = 0 342 self.bound = n 343 self.indexed = indexed
344
345 - def next(self):
346 if self.i >= self.bound: 347 raise StopIteration 348 i = self.i 349 val = self.Q 350 self.i +=1 351 self.Q=self.op(self.Q,self.P) 352 if self.indexed: 353 return (i,val) 354 else: 355 return val
356 - def __iter__(self):
357 return self
358 359
360 -def bsgs(a, b, bounds, operation='*', identity=None, inverse=None, op=None):
361 r""" 362 Totally generic discrete baby-step giant-step function. 363 364 Solves n*a=b (or a^n=b) with lb<=n<=ub where bounds==(lb,ub), 365 raising an error if no such n exists. 366 367 a and b must be elements of some group with given identity, 368 inverse of x given by inverse(x), and group operation on x,y by 369 op(x,y). 370 371 If operation is '*' or '+' then the other 372 arguments are provided automatically; otherwise they must be 373 provided by the caller. 374 375 INPUT: 376 a -- group element 377 b -- group element 378 bounds -- a 2-tuple of integers (lower,upper) with 0<=lower<=upper 379 operation -- string: '*', '+', 'other' 380 identity -- the identity element of the group 381 inverse() -- function of 1 argument x returning inverse of x 382 op() -- function of 2 arguments x,y returning x*y in group 383 384 OUTPUT: 385 Returns an integer $n$ such that $a^n = b$ (or $n*a = b$). 386 If no such $n$ exists, this function raises a ValueError exception. 387 388 NOTE: This is a generalization of discrete logarithm. One 389 situation where this version is useful is to find the order of 390 an element in a group where we only have bounds on the group 391 order (see the elliptic curve example below). 392 393 ALGORITHM: Baby step giant step. Time and space are soft 394 O(sqrt(n)) where n is the difference between upper and lower 395 bounds. 396 397 EXAMPLES: 398 sage: b = Mod(2,37); a = b^20 399 sage: bsgs(b, a, (0,36)) 400 20 401 402 sage: p=next_prime(10^20) 403 sage: a=Mod(2,p); b=a^(10^25) 404 sage: bsgs(a, b, (10^25-10^6,10^25+10^6)) == 10^25 405 True 406 407 sage: K = GF(3^6,'b') 408 sage: a = K.gen() 409 sage: b = a^210 410 sage: bsgs(a, b, (0,K.order()-1)) 411 210 412 413 sage: K.<z>=CyclotomicField(230) 414 sage: w=z^500 415 sage: bsgs(z,w,(0,229)) 416 40 417 418 An additive example in an elliptic curve group: 419 sage: F.<a>=GF(37^5,'a') 420 sage: E=EllipticCurve(F,[1,1]) 421 sage: P=E.lift_x(a); P 422 (a : 9*a^4 + 22*a^3 + 23*a^2 + 30 : 1) 423 424 This will return a multiple of the order of P: 425 sage: bsgs(P,P.parent()(0),Hasse_bounds(F.order()),operation='+') 426 69327408 427 428 AUTHOR: 429 -- John Cremona (2008-03-15) 430 """ 431 Z = integer_ring.ZZ 432 433 from operator import inv, mul, neg, add 434 435 if operation in multiplication_names: 436 identity = a.parent()(1) 437 inverse = inv 438 op = mul 439 elif operation in addition_names: 440 identity = a.parent()(0) 441 inverse = neg 442 op = add 443 else: 444 if identity==None or inverse==None or op==None: 445 raise ValueError, "identity, inverse and operation must be given" 446 447 lb, ub = bounds 448 if lb<0 or ub<lb: 449 raise ValueError, "bsgs() requires 0<=lb<=ub" 450 451 if a.is_zero() and not b.is_zero(): 452 raise ValueError, "No solution in bsgs()" 453 454 ran = 1 + ub - lb # the length of the interval 455 456 c = op(inverse(b),multiple(a,lb,operation=operation)) 457 458 if ran < 30: # use simple search for small ranges 459 i = lb 460 d = c 461 # for i,d in multiples(a,ran,c,indexed=True,operation=operation): 462 for i0 in range(ran): 463 i = lb + i0 464 if identity == d: # identity == b^(-1)*a^i, so return i 465 return Z(i) 466 d = op(a,d) 467 raise ValueError, "No solution in bsgs()" 468 469 m = ran.isqrt()+1 # we need sqrt(ran) rounded up 470 table = dict() # will hold pairs (a^(lb+i),lb+i) for i in range(m) 471 472 d=c 473 for i0 in misc.srange(m): 474 i = lb + i0 475 if identity==d: # identity == b^(-1)*a^i, so return i 476 return Z(i) 477 table[d] = i 478 d=op(d,a) 479 480 c = op(c,inverse(d)) # this is now a**(-m) 481 d=identity 482 for i in misc.srange(m): 483 j = table.get(d) 484 if not j==None: # then d == b*a**(-i*m) == a**j 485 return Z(i*m + j) 486 d=op(c,d) 487 488 raise ValueError, "Log of %s to the base %s does not exist in %s."%(b,a,bounds)
489
490 -def discrete_log(a, base, ord=None, operation='*', identity=None, inverse=None, op=None):
491 """ 492 Totally generic discrete log function. 493 494 a and base must be elements of some group with identity given by 495 identity, inverse of x by inverse(x), and group operation on x,y 496 by op(x,y). 497 498 If operation is '*' or '+' then the other 499 arguments are provided automatically; otherwise they must be 500 provided by the caller. 501 502 WARNING: If x has a log method, it is likely to be vastly faster 503 than using this function. E.g., if x is an integer modulo n, use 504 its log method instead! 505 506 INPUT: 507 a -- group element 508 base -- group element (the base) 509 ord -- integer (multiple of order of base, or None) 510 operation -- string: '*', '+', 'other' 511 identity -- the group's identity 512 inverse() -- function of 1 argument x returning inverse of x 513 op() -- function of 2 arguments x,y returning x*y in group 514 515 OUTPUT: 516 Returns an integer $n$ such that $b^n = a$ (or $n*b = a$), 517 assuming that ord is a multiple of the order of the base $b$. 518 If ord is not specified an attempt is made to compute it. 519 520 If no such $n$ exists, this function raises a ValueError exception. 521 522 ALGORITHM: Baby step giant step. 523 524 EXAMPLES: 525 sage: b = Mod(2,37); a = b^20 526 sage: discrete_log(a, b) 527 20 528 sage: b = Mod(2,997); a = b^20 529 sage: discrete_log(a, b) 530 20 531 532 sage: K = GF(3^6,'b') 533 sage: b = K.gen() 534 sage: a = b^210 535 sage: discrete_log(a, b, K.order()-1) 536 210 537 538 sage: b = Mod(1,37); x = Mod(2,37) 539 sage: discrete_log(x, b) 540 Traceback (most recent call last): 541 ... 542 ValueError: No discrete log of 2 found to base 1 543 sage: b = Mod(1,997); x = Mod(2,997) 544 sage: discrete_log(x, b) 545 Traceback (most recent call last): 546 ... 547 ValueError: No discrete log of 2 found to base 1 548 549 See trac\#2356: 550 sage: F.<w> = GF(121) 551 sage: v = w^120 552 sage: v.log(w) 553 0 554 555 sage: K.<z>=CyclotomicField(230) 556 sage: w=z^50 557 sage: discrete_log(w,z) 558 50 559 560 An example where the order is infinite: note that we must give 561 an upper bound here: 562 sage: K.<a> = QuadraticField(23) 563 sage: eps = 5*a-24 # a fundamental unit 564 sage: eps.multiplicative_order() 565 +Infinity 566 sage: eta = eps^100 567 sage: discrete_log(eta,eps,1000) 568 100 569 570 In this case we cannot detect negative powers: 571 sage: eta = eps^(-3) 572 sage: discrete_log(eta,eps,100) 573 Traceback (most recent call last): 574 ... 575 ValueError: No discrete log of -11515*a - 55224 found to base 5*a - 24 576 577 But we can invert the base (and negate the result) instead: 578 sage: - discrete_log(eta^-1,eps,100) 579 -3 580 581 An additive example: elliptic curve DLOG: 582 sage: F=GF(37^2,'a') 583 sage: E=EllipticCurve(F,[1,1]) 584 sage: F.<a>=GF(37^2,'a') 585 sage: E=EllipticCurve(F,[1,1]) 586 sage: P=E(25*a + 16 , 15*a + 7 ) 587 sage: P.order() 588 672 589 sage: Q=39*P; Q 590 (36*a + 32 : 5*a + 12 : 1) 591 sage: discrete_log(Q,P,P.order(),'+') 592 39 593 594 AUTHOR: 595 -- William Stein and David Joyner (2005-01-05) 596 -- John Cremona (2008-02-29) rewrite using dict() and make generic 597 598 """ 599 if ord == None: 600 if operation in multiplication_names: 601 try: 602 ord = base.multiplicative_order() 603 except: 604 ord = base.order() 605 elif operation in addition_names: 606 try: 607 ord = base.additive_order() 608 except: 609 ord = base.order() 610 else: 611 try: 612 ord = base.order() 613 except: 614 raise ValueError, "ord must be specified" 615 try: 616 return bsgs(base,a,(0,ord),operation=operation) 617 except ValueError: 618 raise ValueError, "No discrete log of %s found to base %s"%(a,base)
619 620 ################################################################ 621 # 622 # Older version of discrete_log. Works fine but has been 623 # superceded by the version which simply calls the more general 624 # bsgs() function. 625 # 626 ################################################################ 627
628 -def old_discrete_log(a, base, ord=None, operation='*', 629 identity=None, inverse=None, op=None):
630 r""" 631 Totally generic discrete log function. 632 633 a and base must be elements of some group with identity given by 634 identity, inverse of x by inverse(x), and group operation on x,y 635 by op(x,y). 636 637 If operation is '*' or '+' then the other 638 arguments are provided automatically; otherwise they must be 639 provided by the caller. 640 641 WARNING: If x has a log method, it is likely to be vastly faster 642 than using this function. E.g., if x is an integer modulo n, use 643 its log method instead! 644 645 INPUT: 646 a -- group element 647 base -- group element (the base) 648 ord -- integer (multiple of order of base, or None) 649 operation -- string: '*', '+', 'other' 650 identity -- the group's identity 651 inverse() -- function of 1 argument x returning inverse of x 652 op() -- function of 2 arguments x,y returning x*y in group 653 654 OUTPUT: 655 Returns an integer $n$ such that $b^n = a$ (or $n*b = a$), 656 assuming that ord is a multiple of the order of the base $b$. 657 If ord is not specified an attempt is made to compute it. 658 659 If no such $n$ exists, this function raises a ValueError exception. 660 661 ALGORITHM: Baby step giant step. 662 663 EXAMPLES: 664 sage: b = Mod(2,37); a = b^20 665 sage: old_discrete_log(a, b) 666 20 667 sage: b = Mod(2,997); a = b^20 668 sage: old_discrete_log(a, b) 669 20 670 671 sage: K = GF(3^6,'b') 672 sage: b = K.gen() 673 sage: a = b^210 674 sage: old_discrete_log(a, b, K.order()-1) 675 210 676 677 sage: b = Mod(1,37); x = Mod(2,37) 678 sage: old_discrete_log(x, b) 679 Traceback (most recent call last): 680 ... 681 ValueError: Log of 2 to the base 1 does not exist. 682 sage: b = Mod(1,997); x = Mod(2,997) 683 sage: old_discrete_log(x, b) 684 Traceback (most recent call last): 685 ... 686 ValueError: Log of 2 to the base 1 does not exist. 687 688 See trac\#2356: 689 sage: F.<w> = GF(121) 690 sage: v = w^120 691 sage: v.log(w) 692 0 693 694 sage: K.<z>=CyclotomicField(230) 695 sage: w=z^50 696 sage: old_discrete_log(w,z) 697 50 698 699 An example where the order is infinite: note that we must give 700 an upper bound here: 701 sage: K.<a> = QuadraticField(23) 702 sage: eps = 5*a-24 # a fundamental unit 703 sage: eps.multiplicative_order() 704 +Infinity 705 sage: eta = eps^100 706 sage: old_discrete_log(eta,eps,1000) 707 100 708 709 In this case we cannot detect negative powers: 710 sage: eta = eps^(-3) 711 sage: old_discrete_log(eta,eps,100) 712 Traceback (most recent call last): 713 ... 714 ValueError: Log of -11515*a - 55224 to the base 5*a - 24 does not exist. 715 716 But we can invert the base (and negate the result) instead: 717 sage: - old_discrete_log(eta^-1,eps,100) 718 -3 719 720 An additive example: elliptic curve DLOG: 721 sage: F=GF(37^2,'a') 722 sage: E=EllipticCurve(F,[1,1]) 723 sage: F.<a>=GF(37^2,'a') 724 sage: E=EllipticCurve(F,[1,1]) 725 sage: P=E(25*a + 16 , 15*a + 7 ) 726 sage: P.order() 727 672 728 sage: Q=39*P; Q 729 (36*a + 32 : 5*a + 12 : 1) 730 sage: old_discrete_log(Q,P,P.order(),'+') 731 39 732 733 AUTHOR: 734 -- William Stein and David Joyner (2005-01-05) 735 -- John Cremona (2008-02-29) rewrite using dict() and make generic 736 """ 737 Z = integer_ring.ZZ 738 b = base 739 740 from operator import inv, mul, neg, add 741 742 if operation in multiplication_names: 743 identity = b.parent()(1) 744 inverse = inv 745 op = mul 746 if ord==None: 747 ord = b.multiplicative_order() 748 elif operation in addition_names: 749 identity = b.parent()(0) 750 inverse = neg 751 op = add 752 if ord==None: 753 ord = b.order() 754 else: 755 if ord==None or identity==None or inverse==None or op==None: 756 raise ValueError, "order, identity, inverse and operation must all be specified" 757 758 ord = Z(ord) 759 if ord < 100: 760 c = identity 761 for i in range(ord): 762 if c == a: # is b^i 763 return Z(i) 764 c = op(c,b) 765 raise ValueError, "Log of %s to the base %s does not exist."%(a,b) 766 767 m = ord.isqrt()+1 # we need sqrt(ord) rounded up 768 table = dict() # will hold pairs (b^j,j) for j in range(m) 769 g = identity # will run through b**j 770 for j in range(m): 771 if a==g: 772 return Z(j) 773 table[g] = j 774 g = op(g,b) 775 776 g = inverse(g) # this is now b**(-m) 777 h = op(a,g) # will run through a*g**i = a*b**(-i*m) 778 for i in range(1,m): 779 j = table.get(h) 780 if not j==None: # then a*b**(-i*m) == b**j 781 return Z(i*m + j) 782 if i < m-1: 783 h = op(h,g) 784 785 raise ValueError, "Log of %s to the base %s does not exist."%(a,b)
786 787 discrete_log_generic = discrete_log 788 789 790 ################################################################ 791 # 792 # Generic linear relation finder 793 # 794 ################################################################ 795
796 -def linear_relation(P, Q, operation='+', identity=None, inverse=None, op=None):
797 """ 798 Function which solves the equation a*P=m*Q or P^a=Q^m. 799 800 Additive version: returns (a,m) with minimal m>0 such that 801 a*P==m*Q. Special case: if <P> and <Q> intersect only in {0} then 802 (a,m)=(0,Q.additive_order()). 803 804 Multiplicative version: returns (a,m) with minimal m>0 such that 805 P^a==Q^m. Special case: if <P> and <Q> intersect only in {1} then 806 (a,m)=(0,Q.multiplicative_order()). 807 808 Works in general finite abelian groups: uses bsgs() 809 810 EXAMPLE: 811 812 An additive example (in an elliptic curve group): 813 sage: F.<a>=GF(3^6,'a') 814 sage: E=EllipticCurve([a^5 + 2*a^3 + 2*a^2 + 2*a, a^4 + a^3 + 2*a + 1]) 815 sage: P=E(a^5 + a^4 + a^3 + a^2 + a + 2 , 0) 816 sage: Q=E(2*a^3 + 2*a^2 + 2*a , a^3 + 2*a^2 + 1) 817 sage: linear_relation(P,Q,'+') 818 (1, 2) 819 sage: P == 2*Q 820 True 821 822 An multiplcative example (in a finite field's multiplicative group): 823 sage: F.<a>=GF(3^6,'a') 824 sage: a.multiplicative_order().factor() 825 2^3 * 7 * 13 826 sage: b=a^7 827 sage: c=a^13 828 sage: linear_relation(b,c,'*') 829 (13, 7) 830 sage: b^13==c^7 831 True 832 """ 833 834 from operator import mul, add 835 Z = integer_ring.ZZ 836 837 if operation in multiplication_names: 838 op = mul 839 try: 840 n = P.multiplicative_order() 841 m = Q.multiplicative_order() 842 except: 843 n = P.order() 844 m = Q.order() 845 elif operation in addition_names: 846 op = add 847 try: 848 n = P.additive_order() 849 m = Q.additive_order() 850 except: 851 n = P.order() 852 m = Q.order() 853 else: 854 if op==None: 855 raise ValueError, "operation must be specified" 856 n = P.order() 857 m = Q.order() 858 859 g = sage.rings.arith.gcd(n,m) 860 if g==1: return (m,Z(0)) 861 n1 = n//g 862 m1 = m//g 863 P1 = multiple(P,n1,operation=operation) # has exact order g 864 Q1 = multiple(Q,m1,operation=operation) # has exact order g 865 866 # now see if Q1 is a multiple of P1; the only multiples we 867 # need check are h*Q1 where h divides g 868 for h in g.divisors(): # positive divisors! 869 try: 870 Q2 = multiple(Q1,h,operation=operation) 871 return (n1 * bsgs(P1,Q2,(0,g-1),operation=operation), 872 m1 * h) 873 except ValueError: 874 pass # to next h 875 raise ValueError, "No solution found in linear_relation!"
876 877 ################################################################ 878 # 879 # Generic functions to find orders of elements 880 # 881 # 1. order_from_multiple: finds the order given a multiple of the order 882 # 883 # 2. order_from_bounds: finds the order given an interval containing a 884 # multiple of the order 885 # 886 ################################################################ 887
888 -def order_from_multiple(P, m, plist=None, operation='+', 889 identity=None, inverse=None, op=None):
890 """ 891 Generic function to find order of P given a multiple of the order 892 893 INPUT: 894 895 P -- a Sage object which is a group element 896 897 m -- a Sage integer which is a multiple of the order of P, 898 i.e. we require that m*P=0 (or P^m=1). 899 900 plist -- a list of the prime factors of m, or None in 901 which case this function will need to factor m. 902 903 operation -- string: '+' (default ) or 904 '*' or other. 905 906 If other, the following must be supplied: 907 908 identity: the identity element for the group; 909 910 inverse(): a function of one argument giving 911 the inverse of a group element; 912 913 op(): a function of 2 arguments defining 914 the group binary operation. 915 916 NOTE: It is more efficient for the caller to factor m and cache the 917 factors for subsequent calls. 918 919 EXAMPLES: 920 sage: k.<a> = GF(5^5) 921 sage: b = a^4 922 sage: order_from_multiple(b,5^5-1,operation='*') 923 781 924 sage: E = EllipticCurve(k,[2,4]) 925 sage: P = E(3*a^4 + 3*a , 2*a + 1 ) 926 sage: M = E.cardinality(); M 927 3227 928 sage: plist = M.prime_factors() 929 sage: order_from_multiple(P, M, plist, operation='+') 930 3227 931 sage: Q = E(0,2) 932 sage: order_from_multiple(Q, M, plist, operation='+') 933 7 934 935 sage: K.<z>=CyclotomicField(230) 936 sage: w=z^50 937 sage: order_from_multiple(w,230,operation='*') 938 23 939 940 """ 941 from operator import mul, add 942 Z = integer_ring.ZZ 943 944 if operation in multiplication_names: 945 op = mul 946 identity = P.parent()(1) 947 elif operation in addition_names: 948 op = add 949 identity = P.parent()(0) 950 else: 951 if op==None: 952 raise ValueError, "operation and identity must be specified" 953 954 M=Z(m) 955 assert multiple(P,M,operation=operation) == identity 956 957 if plist==None: 958 plist = M.prime_factors() 959 960 # For each p in plist we determine the power of p dividing 961 # the order, accumulating the order in N 962 963 N=Z(1) 964 for p in plist: 965 Q = multiple(P,M.prime_to_m_part(p),operation=operation) 966 # so Q has p-power order 967 while Q != identity: 968 Q = multiple(Q,p,operation=operation) 969 N *= p 970 971 # now N is the exact order of self 972 return N
973
974 -def order_from_bounds(P, bounds, d=None, operation='+', 975 identity=None, inverse=None, op=None):
976 """ 977 Generic function to find order of P given bounds on group order. 978 979 upper and lower bounds for a multiple of the order (e.g. bounds on the 980 order of the group of which P is an element) 981 982 INPUT: 983 984 P -- a Sage object which is a group element 985 986 bounds -- a 2-tuple (lb,ub) such that m*P=0 (or P^m=1) for 987 some m with lb<=m<=ub 988 989 d -- (optional) a positive integer; only m which are 990 multiples of this will be considered. 991 992 operation -- string: '+' (default ) or 993 '*' or other. 994 995 If other, the following must be supplied: 996 997 identity: the identity element for the group; 998 999 inverse(): a function of one argument giving 1000 the inverse of a group element; 1001 1002 op(): a function of 2 arguments defining 1003 the group binary operation. 1004 1005 1006 NOTE: Typically lb and ub will be bounds on the group order, 1007 and from previous calculation we know that the group order is 1008 divisible by d. 1009 1010 EXAMPLES: 1011 sage: k.<a> = GF(5^5) 1012 sage: b = a^4 1013 sage: order_from_bounds(b,(5^4,5^5),operation='*') 1014 781 1015 sage: E = EllipticCurve(k,[2,4]) 1016 sage: P = E(3*a^4 + 3*a , 2*a + 1 ) 1017 sage: bounds = Hasse_bounds(5^5) 1018 sage: Q = E(0,2) 1019 sage: order_from_bounds(Q, bounds, operation='+') 1020 7 1021 sage: order_from_bounds(P, bounds, 7, operation='+') 1022 3227 1023 1024 sage: K.<z>=CyclotomicField(230) 1025 sage: w=z^50 1026 sage: order_from_bounds(w,(200,250),operation='*') 1027 23 1028 1029 """ 1030 from operator import mul, add 1031 Z = integer_ring.ZZ 1032 1033 if operation in multiplication_names: 1034 op = mul 1035 identity = P.parent()(1) 1036 elif operation in addition_names: 1037 op = add 1038 identity = P.parent()(0) 1039 else: 1040 if op==None: 1041 raise ValueError, "operation and identity must be specified" 1042 1043 Q = P 1044 if d == None: d = 1 1045 if d > 1: 1046 Q = multiple(P,d,operation=operation) 1047 lb, ub = bounds 1048 bounds = ( sage.rings.arith.integer_ceil(lb/d), 1049 sage.rings.arith.integer_floor(ub/d) ) 1050 1051 # Use generic bsgs to find n=d*m with lb<=n<=ub and n*P=0 1052 1053 m = d * bsgs(Q, identity, bounds, operation=operation) 1054 1055 # Now use the order_from_multiple() function to finish the job: 1056 1057 return order_from_multiple(P, m, operation=operation)
1058
1059 -def merge_points(P1,P2, operation='+', 1060 identity=None, inverse=None, op=None):
1061 """ 1062 Returns a group element whose order is the lcm of the given elements 1063 1064 INPUT: 1065 P1 -- a pair (g1,n1) where g1 is a group element of order n1 1066 P2 -- a pair (g2,n2) where g2 is a group element of order n2 1067 operation -- string: '+' (default ) or 1068 '*' or other. 1069 1070 If other, the following must be supplied: 1071 1072 identity: the identity element for the group; 1073 1074 inverse(): a function of one argument giving 1075 the inverse of a group element; 1076 1077 op(): a function of 2 arguments defining 1078 the group binary operation. 1079 1080 1081 OUTPUT: 1082 A pair (g3,n3) where g3 has order n3=lcm(n1,n2) 1083 1084 EXAMPLES: 1085 sage: F.<a>=GF(3^6,'a') 1086 sage: b = a^7 1087 sage: c = a^13 1088 sage: ob = (3^6-1)//7 1089 sage: oc = (3^6-1)//13 1090 sage: merge_points((b,ob),(c,oc),operation='*') 1091 (a^4 + 2*a^3 + 2*a^2, 728) 1092 sage: d,od = merge_points((b,ob),(c,oc),operation='*') 1093 sage: od == d.multiplicative_order() 1094 True 1095 sage: od == lcm(ob,oc) 1096 True 1097 1098 sage: E=EllipticCurve([a^5 + 2*a^3 + 2*a^2 + 2*a, a^4 + a^3 + 2*a + 1]) 1099 sage: P=E(2*a^5 + 2*a^4 + a^3 + 2 , a^4 + a^3 + a^2 + 2*a + 2) 1100 sage: P.order() 1101 7 1102 sage: Q=E(2*a^5 + 2*a^4 + 1 , a^5 + 2*a^3 + 2*a + 2 ) 1103 sage: Q.order() 1104 4 1105 sage: R,m = merge_points((P,7),(Q,4), operation='+') 1106 sage: R.order() == m 1107 True 1108 sage: m == lcm(7,4) 1109 True 1110 """ 1111 from operator import mul, add 1112 Z = integer_ring.ZZ 1113 1114 g1, n1 = P1 1115 g2, n2 = P2 1116 1117 if operation in multiplication_names: 1118 op = mul 1119 identity = g1.parent()(1) 1120 elif operation in addition_names: 1121 op = add 1122 identity = g1.parent()(0) 1123 else: 1124 if op==None: 1125 raise ValueError, "operation and identity must be specified" 1126 1127 assert multiple(g1,n1,operation=operation) == identity 1128 assert multiple(g2,n2,operation=operation) == identity 1129 1130 # trivial cases 1131 if n1.divides(n2): 1132 return (g2,n2) 1133 if n2.divides(n1): 1134 return (g1,n1) 1135 1136 m,k1,k2 = sage.rings.arith.xlcm(n1,n2); 1137 m1 = n1//k1 1138 m2 = n2//k2 1139 g1 = multiple(g1,m1,operation=operation) 1140 g2 = multiple(g2,m2,operation=operation) 1141 return (op(g1,g2), m)
1142