advmpz.rst
1 Multiple-precision Integers (Advanced topics) 2 ============================================= 3 4 The xmpz type 5 ------------- 6 7 gmpy2 provides access to an experimental integer type called *xmpz*. The 8 *xmpz* type is a mutable integer type. In-place operations (+=, //=, etc.) 9 modify the original object and do not create a new object. Instances of 10 *xmpz* cannot be used as dictionary keys. 11 12 :: 13 14 >>> import gmpy2 15 >>> from gmpy2 import xmpz 16 >>> a = xmpz(123) 17 >>> b = a 18 >>> a += 1 19 >>> a 20 xmpz(124) 21 >>> b 22 xmpz(124) 23 24 The ability to change an *xmpz* object in-place allows for efficient and rapid 25 bit manipulation. 26 27 Individual bits can be set or cleared:: 28 29 >>> a[10]=1 30 >>> a 31 xmpz(1148) 32 33 Slice notation is supported. The bits referenced by a slice can be either 'read 34 from' or 'written to'. To clear a slice of bits, use a source value of 0. In 35 2s-complement format, 0 is represented by an arbitrary number of 0-bits. To set 36 a slice of bits, use a source value of ~0. The *tilde* operator inverts, or 37 complements the bits in an integer. (~0 is -1 so you can also use -1.) In 38 2s-complement format, -1 is represented by an arbitrary number of 1-bits. 39 40 If a value for *stop* is specified in a slice assignment and the actual 41 bit-length of the *xmpz* is less than *stop*, then the destination *xmpz* is 42 logically padded with 0-bits to length *stop*. 43 44 :: 45 46 >>> a=xmpz(0) 47 >>> a[8:16] = ~0 48 >>> bin(a) 49 '0b1111111100000000' 50 >>> a[4:12] = ~a[4:12] 51 >>> bin(a) 52 '0b1111000011110000' 53 54 Bits can be reversed:: 55 56 >>> bin(a) 57 '0b10001111100' 58 >>> a[::] = a[::-1] 59 >>> bin(a) 60 '0b111110001' 61 62 The *iter_bits()* method returns a generator that returns True or False for each 63 bit position. The methods *iter_clear()*, and *iter_set()* return generators 64 that return the bit positions that are 1 or 0. The methods support arguments 65 *start* and *stop* that define the beginning and ending bit positions that are 66 used. To mimic the behavior of slices. the bit positions checked include *start* 67 but the last position checked is *stop* - 1. 68 69 :: 70 71 >>> a=xmpz(117) 72 >>> bin(a) 73 '0b1110101' 74 >>> list(a.iter_bits()) 75 [True, False, True, False, True, True, True] 76 >>> list(a.iter_clear()) 77 [1, 3] 78 >>> list(a.iter_set()) 79 [0, 2, 4, 5, 6] 80 >>> list(a.iter_bits(stop=12)) 81 [True, False, True, False, True, True, True, False, False, False, False, False] 82 83 The following program uses the Sieve of Eratosthenes to generate a list of 84 prime numbers. 85 86 :: 87 88 from __future__ import print_function 89 import time 90 import gmpy2 91 92 def sieve(limit=1000000): 93 '''Returns a generator that yields the prime numbers up to limit.''' 94 95 # Increment by 1 to account for the fact that slices do not include 96 # the last index value but we do want to include the last value for 97 # calculating a list of primes. 98 sieve_limit = gmpy2.isqrt(limit) + 1 99 limit += 1 100 101 # Mark bit positions 0 and 1 as not prime. 102 bitmap = gmpy2.xmpz(3) 103 104 # Process 2 separately. This allows us to use p+p for the step size 105 # when sieving the remaining primes. 106 bitmap[4 : limit : 2] = -1 107 108 # Sieve the remaining primes. 109 for p in bitmap.iter_clear(3, sieve_limit): 110 bitmap[p*p : limit : p+p] = -1 111 112 return bitmap.iter_clear(2, limit) 113 114 if __name__ == "__main__": 115 start = time.time() 116 result = list(sieve()) 117 print(time.time() - start) 118 print(len(result)) 119 120 121 Advanced Number Theory Functions 122 -------------------------------- 123 124 The following functions are based on mpz_lucas.c and mpz_prp.c by David 125 Cleaver. 126 127 A good reference for probable prime testing is 128 http://www.pseudoprime.com/pseudo.html 129 130 **is_bpsw_prp(...)** 131 is_bpsw_prp(n) will return True if *n* is a Baillie-Pomerance-Selfridge-Wagstaff 132 probable prime. A BPSW probable prime passes the is_strong_prp() test with base 133 2 and the is_selfridge_prp() test. 134 135 **is_euler_prp(...)** 136 is_euler_prp(n,a) will return True if *n* is an Euler (also known as 137 Solovay-Strassen) probable prime to the base *a*. 138 139 | Assuming: 140 | gcd(n, a) == 1 141 | n is odd 142 | 143 | Then an Euler probable prime requires: 144 | a**((n-1)/2) == 1 (mod n) 145 146 **is_extra_strong_lucas_prp(...)** 147 is_extra_strong_lucas_prp(n,p) will return True if *n* is an extra strong 148 Lucas probable prime with parameters (p,1). 149 150 | Assuming: 151 | n is odd 152 | D = p*p - 4, D != 0 153 | gcd(n, 2*D) == 1 154 | n = s*(2**r) + Jacobi(D,n), s odd 155 | 156 | Then an extra strong Lucas probable prime requires: 157 | lucasu(p,1,s) == 0 (mod n) 158 | or 159 | lucasv(p,1,s) == +/-2 (mod n) 160 | or 161 | lucasv(p,1,s*(2**t)) == 0 (mod n) for some t, 0 <= t < r 162 163 **is_fermat_prp(...)** 164 is_fermat_prp(n,a) will return True if *n* is a Fermat probable prime to the 165 base a. 166 167 | Assuming: 168 | gcd(n,a) == 1 169 | 170 | Then a Fermat probable prime requires: 171 | a**(n-1) == 1 (mod n) 172 173 **is_fibonacci_prp(...)** 174 is_fibonacci_prp(n,p,q) will return True if *n* is a Fibonacci 175 probable prime with parameters (p,q). 176 177 | Assuming: 178 | n is odd 179 | p > 0, q = +/-1 180 | p*p - 4*q != 0 181 | 182 | Then a Fibonacci probable prime requires: 183 | lucasv(p,q,n) == p (mod n). 184 185 **is_lucas_prp(...)** 186 is_lucas_prp(n,p,q) will return True if *n* is a Lucas probable prime with 187 parameters (p,q). 188 189 | Assuming: 190 | n is odd 191 | D = p*p - 4*q, D != 0 192 | gcd(n, 2*q*D) == 1 193 | 194 | Then a Lucas probable prime requires: 195 | lucasu(p,q,n - Jacobi(D,n)) == 0 (mod n) 196 197 **is_selfridge_prp(...)** 198 is_selfridge_prp(n) will return True if *n* is a Lucas probable prime with 199 Selfidge parameters (p,q). The Selfridge parameters are chosen by finding 200 the first element D in the sequence {5, -7, 9, -11, 13, ...} such that 201 Jacobi(D,n) == -1. Let p=1 and q = (1-D)/4 and then perform a Lucas 202 probable prime test. 203 204 **is_strong_bpsw_prp(...)** 205 is_strong_bpsw_prp(n) will return True if *n* is a strong 206 Baillie-Pomerance-Selfridge-Wagstaff probable prime. A strong BPSW 207 probable prime passes the is_strong_prp() test with base 2 and the 208 is_strongselfridge_prp() test. 209 210 **is_strong_lucas_prp(...)** 211 is_strong_lucas_prp(n,p,q) will return True if *n* is a strong Lucas 212 probable prime with parameters (p,q). 213 214 | Assuming: 215 | n is odd 216 | D = p*p - 4*q, D != 0 217 | gcd(n, 2*q*D) == 1 218 | n = s*(2**r) + Jacobi(D,n), s odd 219 | 220 | Then a strong Lucas probable prime requires: 221 | lucasu(p,q,s) == 0 (mod n) 222 | or 223 | lucasv(p,q,s*(2**t)) == 0 (mod n) for some t, 0 <= t < r 224 225 **is_strong_prp(...)** 226 is_strong_prp(n,a) will return True if *n* is a strong (also known as 227 Miller-Rabin) probable prime to the base a. 228 229 | Assuming: 230 | gcd(n,a) == 1 231 | n is odd 232 | n = s*(2**r) + 1, with s odd 233 | 234 | Then a strong probable prime requires one of the following is true: 235 | a**s == 1 (mod n) 236 | or 237 | a**(s*(2**t)) == -1 (mod n) for some t, 0 <= t < r. 238 239 **is_strong_selfridge_prp(...)** 240 is_strong_selfridge_prp(n) will return True if *n* is a strong Lucas 241 probable prime with Selfidge parameters (p,q). The Selfridge parameters are 242 chosen by finding the first element D in the sequence 243 {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) == -1. Let p=1 and 244 q = (1-D)/4 and then perform a strong Lucas probable prime test. 245 246 **lucasu(...)** 247 lucasu(p,q,k) will return the k-th element of the Lucas U sequence defined 248 by p,q. p*p - 4*q must not equal 0; k must be greater than or equal to 0. 249 250 **lucasu_mod(...)** 251 lucasu_mod(p,q,k,n) will return the k-th element of the Lucas U sequence 252 defined by p,q (mod n). p*p - 4*q must not equal 0; k must be greater than 253 or equal to 0; n must be greater than 0. 254 255 **lucasv(...)** 256 lucasv(p,q,k) will return the k-th element of the Lucas V sequence defined 257 by parameters (p,q). p*p - 4*q must not equal 0; k must be greater than or 258 equal to 0. 259 260 **lucasv_mod(...)** 261 lucasv_mod(p,q,k,n) will return the k-th element of the Lucas V sequence 262 defined by parameters (p,q) (mod n). p*p - 4*q must not equal 0; k must be 263 greater than or equal to 0; n must be greater than 0. 264