_IntegerGMP.py 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799
  1. # ===================================================================
  2. #
  3. # Copyright (c) 2014, Legrandin <helderijs@gmail.com>
  4. # All rights reserved.
  5. #
  6. # Redistribution and use in source and binary forms, with or without
  7. # modification, are permitted provided that the following conditions
  8. # are met:
  9. #
  10. # 1. Redistributions of source code must retain the above copyright
  11. # notice, this list of conditions and the following disclaimer.
  12. # 2. Redistributions in binary form must reproduce the above copyright
  13. # notice, this list of conditions and the following disclaimer in
  14. # the documentation and/or other materials provided with the
  15. # distribution.
  16. #
  17. # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  18. # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  19. # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  20. # FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  21. # COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  22. # INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  23. # BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  24. # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  25. # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  26. # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  27. # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28. # POSSIBILITY OF SUCH DAMAGE.
  29. # ===================================================================
  30. import sys
  31. import struct
  32. from Crypto.Util.py3compat import is_native_int
  33. from Crypto.Util._raw_api import (backend, load_lib,
  34. c_ulong, c_size_t, c_uint8_ptr)
  35. from ._IntegerBase import IntegerBase
  36. gmp_defs = """typedef unsigned long UNIX_ULONG;
  37. typedef struct { int a; int b; void *c; } MPZ;
  38. typedef MPZ mpz_t[1];
  39. typedef UNIX_ULONG mp_bitcnt_t;
  40. void __gmpz_init (mpz_t x);
  41. void __gmpz_init_set (mpz_t rop, const mpz_t op);
  42. void __gmpz_init_set_ui (mpz_t rop, UNIX_ULONG op);
  43. UNIX_ULONG __gmpz_get_ui (const mpz_t op);
  44. void __gmpz_set (mpz_t rop, const mpz_t op);
  45. void __gmpz_set_ui (mpz_t rop, UNIX_ULONG op);
  46. void __gmpz_add (mpz_t rop, const mpz_t op1, const mpz_t op2);
  47. void __gmpz_add_ui (mpz_t rop, const mpz_t op1, UNIX_ULONG op2);
  48. void __gmpz_sub_ui (mpz_t rop, const mpz_t op1, UNIX_ULONG op2);
  49. void __gmpz_addmul (mpz_t rop, const mpz_t op1, const mpz_t op2);
  50. void __gmpz_addmul_ui (mpz_t rop, const mpz_t op1, UNIX_ULONG op2);
  51. void __gmpz_submul_ui (mpz_t rop, const mpz_t op1, UNIX_ULONG op2);
  52. void __gmpz_import (mpz_t rop, size_t count, int order, size_t size,
  53. int endian, size_t nails, const void *op);
  54. void * __gmpz_export (void *rop, size_t *countp, int order,
  55. size_t size,
  56. int endian, size_t nails, const mpz_t op);
  57. size_t __gmpz_sizeinbase (const mpz_t op, int base);
  58. void __gmpz_sub (mpz_t rop, const mpz_t op1, const mpz_t op2);
  59. void __gmpz_mul (mpz_t rop, const mpz_t op1, const mpz_t op2);
  60. void __gmpz_mul_ui (mpz_t rop, const mpz_t op1, UNIX_ULONG op2);
  61. int __gmpz_cmp (const mpz_t op1, const mpz_t op2);
  62. void __gmpz_powm (mpz_t rop, const mpz_t base, const mpz_t exp, const
  63. mpz_t mod);
  64. void __gmpz_powm_ui (mpz_t rop, const mpz_t base, UNIX_ULONG exp,
  65. const mpz_t mod);
  66. void __gmpz_pow_ui (mpz_t rop, const mpz_t base, UNIX_ULONG exp);
  67. void __gmpz_sqrt(mpz_t rop, const mpz_t op);
  68. void __gmpz_mod (mpz_t r, const mpz_t n, const mpz_t d);
  69. void __gmpz_neg (mpz_t rop, const mpz_t op);
  70. void __gmpz_abs (mpz_t rop, const mpz_t op);
  71. void __gmpz_and (mpz_t rop, const mpz_t op1, const mpz_t op2);
  72. void __gmpz_ior (mpz_t rop, const mpz_t op1, const mpz_t op2);
  73. void __gmpz_clear (mpz_t x);
  74. void __gmpz_tdiv_q_2exp (mpz_t q, const mpz_t n, mp_bitcnt_t b);
  75. void __gmpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d);
  76. void __gmpz_mul_2exp (mpz_t rop, const mpz_t op1, mp_bitcnt_t op2);
  77. int __gmpz_tstbit (const mpz_t op, mp_bitcnt_t bit_index);
  78. int __gmpz_perfect_square_p (const mpz_t op);
  79. int __gmpz_jacobi (const mpz_t a, const mpz_t b);
  80. void __gmpz_gcd (mpz_t rop, const mpz_t op1, const mpz_t op2);
  81. UNIX_ULONG __gmpz_gcd_ui (mpz_t rop, const mpz_t op1,
  82. UNIX_ULONG op2);
  83. void __gmpz_lcm (mpz_t rop, const mpz_t op1, const mpz_t op2);
  84. int __gmpz_invert (mpz_t rop, const mpz_t op1, const mpz_t op2);
  85. int __gmpz_divisible_p (const mpz_t n, const mpz_t d);
  86. int __gmpz_divisible_ui_p (const mpz_t n, UNIX_ULONG d);
  87. size_t __gmpz_size (const mpz_t op);
  88. UNIX_ULONG __gmpz_getlimbn (const mpz_t op, size_t n);
  89. """
  90. if sys.platform == "win32":
  91. raise ImportError("Not using GMP on Windows")
  92. lib = load_lib("gmp", gmp_defs)
  93. implementation = {"library": "gmp", "api": backend}
  94. if hasattr(lib, "__mpir_version"):
  95. raise ImportError("MPIR library detected")
  96. # Lazy creation of GMP methods
  97. class _GMP(object):
  98. def __getattr__(self, name):
  99. if name.startswith("mpz_"):
  100. func_name = "__gmpz_" + name[4:]
  101. elif name.startswith("gmp_"):
  102. func_name = "__gmp_" + name[4:]
  103. else:
  104. raise AttributeError("Attribute %s is invalid" % name)
  105. func = getattr(lib, func_name)
  106. setattr(self, name, func)
  107. return func
  108. _gmp = _GMP()
  109. # In order to create a function that returns a pointer to
  110. # a new MPZ structure, we need to break the abstraction
  111. # and know exactly what ffi backend we have
  112. if implementation["api"] == "ctypes":
  113. from ctypes import Structure, c_int, c_void_p, byref
  114. class _MPZ(Structure):
  115. _fields_ = [('_mp_alloc', c_int),
  116. ('_mp_size', c_int),
  117. ('_mp_d', c_void_p)]
  118. def new_mpz():
  119. return byref(_MPZ())
  120. _gmp.mpz_getlimbn.restype = c_ulong
  121. else:
  122. # We are using CFFI
  123. from Crypto.Util._raw_api import ffi
  124. def new_mpz():
  125. return ffi.new("MPZ*")
  126. # Size of a native word
  127. _sys_bits = 8 * struct.calcsize("P")
  128. class IntegerGMP(IntegerBase):
  129. """A fast, arbitrary precision integer"""
  130. _zero_mpz_p = new_mpz()
  131. _gmp.mpz_init_set_ui(_zero_mpz_p, c_ulong(0))
  132. def __init__(self, value):
  133. """Initialize the integer to the given value."""
  134. self._mpz_p = new_mpz()
  135. self._initialized = False
  136. if isinstance(value, float):
  137. raise ValueError("A floating point type is not a natural number")
  138. if is_native_int(value):
  139. _gmp.mpz_init(self._mpz_p)
  140. self._initialized = True
  141. if value == 0:
  142. return
  143. tmp = new_mpz()
  144. _gmp.mpz_init(tmp)
  145. try:
  146. positive = value >= 0
  147. reduce = abs(value)
  148. slots = (reduce.bit_length() - 1) // 32 + 1
  149. while slots > 0:
  150. slots = slots - 1
  151. _gmp.mpz_set_ui(tmp,
  152. c_ulong(0xFFFFFFFF & (reduce >> (slots * 32))))
  153. _gmp.mpz_mul_2exp(tmp, tmp, c_ulong(slots * 32))
  154. _gmp.mpz_add(self._mpz_p, self._mpz_p, tmp)
  155. finally:
  156. _gmp.mpz_clear(tmp)
  157. if not positive:
  158. _gmp.mpz_neg(self._mpz_p, self._mpz_p)
  159. elif isinstance(value, IntegerGMP):
  160. _gmp.mpz_init_set(self._mpz_p, value._mpz_p)
  161. self._initialized = True
  162. else:
  163. raise NotImplementedError
  164. # Conversions
  165. def __int__(self):
  166. tmp = new_mpz()
  167. _gmp.mpz_init_set(tmp, self._mpz_p)
  168. try:
  169. value = 0
  170. slot = 0
  171. while _gmp.mpz_cmp(tmp, self._zero_mpz_p) != 0:
  172. lsb = _gmp.mpz_get_ui(tmp) & 0xFFFFFFFF
  173. value |= lsb << (slot * 32)
  174. _gmp.mpz_tdiv_q_2exp(tmp, tmp, c_ulong(32))
  175. slot = slot + 1
  176. finally:
  177. _gmp.mpz_clear(tmp)
  178. if self < 0:
  179. value = -value
  180. return int(value)
  181. def __str__(self):
  182. return str(int(self))
  183. def __repr__(self):
  184. return "Integer(%s)" % str(self)
  185. # Only Python 2.x
  186. def __hex__(self):
  187. return hex(int(self))
  188. # Only Python 3.x
  189. def __index__(self):
  190. return int(self)
  191. def to_bytes(self, block_size=0, byteorder='big'):
  192. """Convert the number into a byte string.
  193. This method encodes the number in network order and prepends
  194. as many zero bytes as required. It only works for non-negative
  195. values.
  196. :Parameters:
  197. block_size : integer
  198. The exact size the output byte string must have.
  199. If zero, the string has the minimal length.
  200. byteorder : string
  201. 'big' for big-endian integers (default), 'little' for litte-endian.
  202. :Returns:
  203. A byte string.
  204. :Raise ValueError:
  205. If the value is negative or if ``block_size`` is
  206. provided and the length of the byte string would exceed it.
  207. """
  208. if self < 0:
  209. raise ValueError("Conversion only valid for non-negative numbers")
  210. num_limbs = _gmp.mpz_size(self._mpz_p)
  211. if _sys_bits == 32:
  212. spchar = "L"
  213. num_limbs = max(1, num_limbs, (block_size + 3) // 4)
  214. elif _sys_bits == 64:
  215. spchar = "Q"
  216. num_limbs = max(1, num_limbs, (block_size + 7) // 8)
  217. else:
  218. raise ValueError("Unknown limb size")
  219. # mpz_getlimbn returns 0 if i is larger than the number of actual limbs
  220. limbs = [_gmp.mpz_getlimbn(self._mpz_p, num_limbs - i - 1) for i in range(num_limbs)]
  221. result = struct.pack(">" + spchar * num_limbs, *limbs)
  222. cutoff_len = len(result) - block_size
  223. if block_size == 0:
  224. result = result.lstrip(b'\x00')
  225. elif cutoff_len > 0:
  226. if result[:cutoff_len] != b'\x00' * (cutoff_len):
  227. raise ValueError("Number is too big to convert to "
  228. "byte string of prescribed length")
  229. result = result[cutoff_len:]
  230. elif cutoff_len < 0:
  231. result = b'\x00' * (-cutoff_len) + result
  232. if byteorder == 'little':
  233. result = result[::-1]
  234. elif byteorder == 'big':
  235. pass
  236. else:
  237. raise ValueError("Incorrect byteorder")
  238. if len(result) == 0:
  239. result = b'\x00'
  240. return result
  241. @staticmethod
  242. def from_bytes(byte_string, byteorder='big'):
  243. """Convert a byte string into a number.
  244. :Parameters:
  245. byte_string : byte string
  246. The input number, encoded in network order.
  247. It can only be non-negative.
  248. byteorder : string
  249. 'big' for big-endian integers (default), 'little' for litte-endian.
  250. :Return:
  251. The ``Integer`` object carrying the same value as the input.
  252. """
  253. result = IntegerGMP(0)
  254. if byteorder == 'big':
  255. pass
  256. elif byteorder == 'little':
  257. byte_string = bytearray(byte_string)
  258. byte_string.reverse()
  259. else:
  260. raise ValueError("Incorrect byteorder")
  261. _gmp.mpz_import(
  262. result._mpz_p,
  263. c_size_t(len(byte_string)), # Amount of words to read
  264. 1, # Big endian
  265. c_size_t(1), # Each word is 1 byte long
  266. 0, # Endianess within a word - not relevant
  267. c_size_t(0), # No nails
  268. c_uint8_ptr(byte_string))
  269. return result
  270. # Relations
  271. def _apply_and_return(self, func, term):
  272. if not isinstance(term, IntegerGMP):
  273. term = IntegerGMP(term)
  274. return func(self._mpz_p, term._mpz_p)
  275. def __eq__(self, term):
  276. if not (isinstance(term, IntegerGMP) or is_native_int(term)):
  277. return False
  278. return self._apply_and_return(_gmp.mpz_cmp, term) == 0
  279. def __ne__(self, term):
  280. if not (isinstance(term, IntegerGMP) or is_native_int(term)):
  281. return True
  282. return self._apply_and_return(_gmp.mpz_cmp, term) != 0
  283. def __lt__(self, term):
  284. return self._apply_and_return(_gmp.mpz_cmp, term) < 0
  285. def __le__(self, term):
  286. return self._apply_and_return(_gmp.mpz_cmp, term) <= 0
  287. def __gt__(self, term):
  288. return self._apply_and_return(_gmp.mpz_cmp, term) > 0
  289. def __ge__(self, term):
  290. return self._apply_and_return(_gmp.mpz_cmp, term) >= 0
  291. def __nonzero__(self):
  292. return _gmp.mpz_cmp(self._mpz_p, self._zero_mpz_p) != 0
  293. __bool__ = __nonzero__
  294. def is_negative(self):
  295. return _gmp.mpz_cmp(self._mpz_p, self._zero_mpz_p) < 0
  296. # Arithmetic operations
  297. def __add__(self, term):
  298. result = IntegerGMP(0)
  299. if not isinstance(term, IntegerGMP):
  300. try:
  301. term = IntegerGMP(term)
  302. except NotImplementedError:
  303. return NotImplemented
  304. _gmp.mpz_add(result._mpz_p,
  305. self._mpz_p,
  306. term._mpz_p)
  307. return result
  308. def __sub__(self, term):
  309. result = IntegerGMP(0)
  310. if not isinstance(term, IntegerGMP):
  311. try:
  312. term = IntegerGMP(term)
  313. except NotImplementedError:
  314. return NotImplemented
  315. _gmp.mpz_sub(result._mpz_p,
  316. self._mpz_p,
  317. term._mpz_p)
  318. return result
  319. def __mul__(self, term):
  320. result = IntegerGMP(0)
  321. if not isinstance(term, IntegerGMP):
  322. try:
  323. term = IntegerGMP(term)
  324. except NotImplementedError:
  325. return NotImplemented
  326. _gmp.mpz_mul(result._mpz_p,
  327. self._mpz_p,
  328. term._mpz_p)
  329. return result
  330. def __floordiv__(self, divisor):
  331. if not isinstance(divisor, IntegerGMP):
  332. divisor = IntegerGMP(divisor)
  333. if _gmp.mpz_cmp(divisor._mpz_p,
  334. self._zero_mpz_p) == 0:
  335. raise ZeroDivisionError("Division by zero")
  336. result = IntegerGMP(0)
  337. _gmp.mpz_fdiv_q(result._mpz_p,
  338. self._mpz_p,
  339. divisor._mpz_p)
  340. return result
  341. def __mod__(self, divisor):
  342. if not isinstance(divisor, IntegerGMP):
  343. divisor = IntegerGMP(divisor)
  344. comp = _gmp.mpz_cmp(divisor._mpz_p,
  345. self._zero_mpz_p)
  346. if comp == 0:
  347. raise ZeroDivisionError("Division by zero")
  348. if comp < 0:
  349. raise ValueError("Modulus must be positive")
  350. result = IntegerGMP(0)
  351. _gmp.mpz_mod(result._mpz_p,
  352. self._mpz_p,
  353. divisor._mpz_p)
  354. return result
  355. def inplace_pow(self, exponent, modulus=None):
  356. if modulus is None:
  357. if exponent < 0:
  358. raise ValueError("Exponent must not be negative")
  359. # Normal exponentiation
  360. if exponent > 256:
  361. raise ValueError("Exponent is too big")
  362. _gmp.mpz_pow_ui(self._mpz_p,
  363. self._mpz_p, # Base
  364. c_ulong(int(exponent))
  365. )
  366. else:
  367. # Modular exponentiation
  368. if not isinstance(modulus, IntegerGMP):
  369. modulus = IntegerGMP(modulus)
  370. if not modulus:
  371. raise ZeroDivisionError("Division by zero")
  372. if modulus.is_negative():
  373. raise ValueError("Modulus must be positive")
  374. if is_native_int(exponent):
  375. if exponent < 0:
  376. raise ValueError("Exponent must not be negative")
  377. if exponent < 65536:
  378. _gmp.mpz_powm_ui(self._mpz_p,
  379. self._mpz_p,
  380. c_ulong(exponent),
  381. modulus._mpz_p)
  382. return self
  383. exponent = IntegerGMP(exponent)
  384. elif exponent.is_negative():
  385. raise ValueError("Exponent must not be negative")
  386. _gmp.mpz_powm(self._mpz_p,
  387. self._mpz_p,
  388. exponent._mpz_p,
  389. modulus._mpz_p)
  390. return self
  391. def __pow__(self, exponent, modulus=None):
  392. result = IntegerGMP(self)
  393. return result.inplace_pow(exponent, modulus)
  394. def __abs__(self):
  395. result = IntegerGMP(0)
  396. _gmp.mpz_abs(result._mpz_p, self._mpz_p)
  397. return result
  398. def sqrt(self, modulus=None):
  399. """Return the largest Integer that does not
  400. exceed the square root"""
  401. if modulus is None:
  402. if self < 0:
  403. raise ValueError("Square root of negative value")
  404. result = IntegerGMP(0)
  405. _gmp.mpz_sqrt(result._mpz_p,
  406. self._mpz_p)
  407. else:
  408. if modulus <= 0:
  409. raise ValueError("Modulus must be positive")
  410. modulus = int(modulus)
  411. result = IntegerGMP(self._tonelli_shanks(int(self) % modulus, modulus))
  412. return result
  413. def __iadd__(self, term):
  414. if is_native_int(term):
  415. if 0 <= term < 65536:
  416. _gmp.mpz_add_ui(self._mpz_p,
  417. self._mpz_p,
  418. c_ulong(term))
  419. return self
  420. if -65535 < term < 0:
  421. _gmp.mpz_sub_ui(self._mpz_p,
  422. self._mpz_p,
  423. c_ulong(-term))
  424. return self
  425. term = IntegerGMP(term)
  426. _gmp.mpz_add(self._mpz_p,
  427. self._mpz_p,
  428. term._mpz_p)
  429. return self
  430. def __isub__(self, term):
  431. if is_native_int(term):
  432. if 0 <= term < 65536:
  433. _gmp.mpz_sub_ui(self._mpz_p,
  434. self._mpz_p,
  435. c_ulong(term))
  436. return self
  437. if -65535 < term < 0:
  438. _gmp.mpz_add_ui(self._mpz_p,
  439. self._mpz_p,
  440. c_ulong(-term))
  441. return self
  442. term = IntegerGMP(term)
  443. _gmp.mpz_sub(self._mpz_p,
  444. self._mpz_p,
  445. term._mpz_p)
  446. return self
  447. def __imul__(self, term):
  448. if is_native_int(term):
  449. if 0 <= term < 65536:
  450. _gmp.mpz_mul_ui(self._mpz_p,
  451. self._mpz_p,
  452. c_ulong(term))
  453. return self
  454. if -65535 < term < 0:
  455. _gmp.mpz_mul_ui(self._mpz_p,
  456. self._mpz_p,
  457. c_ulong(-term))
  458. _gmp.mpz_neg(self._mpz_p, self._mpz_p)
  459. return self
  460. term = IntegerGMP(term)
  461. _gmp.mpz_mul(self._mpz_p,
  462. self._mpz_p,
  463. term._mpz_p)
  464. return self
  465. def __imod__(self, divisor):
  466. if not isinstance(divisor, IntegerGMP):
  467. divisor = IntegerGMP(divisor)
  468. comp = _gmp.mpz_cmp(divisor._mpz_p,
  469. divisor._zero_mpz_p)
  470. if comp == 0:
  471. raise ZeroDivisionError("Division by zero")
  472. if comp < 0:
  473. raise ValueError("Modulus must be positive")
  474. _gmp.mpz_mod(self._mpz_p,
  475. self._mpz_p,
  476. divisor._mpz_p)
  477. return self
  478. # Boolean/bit operations
  479. def __and__(self, term):
  480. result = IntegerGMP(0)
  481. if not isinstance(term, IntegerGMP):
  482. term = IntegerGMP(term)
  483. _gmp.mpz_and(result._mpz_p,
  484. self._mpz_p,
  485. term._mpz_p)
  486. return result
  487. def __or__(self, term):
  488. result = IntegerGMP(0)
  489. if not isinstance(term, IntegerGMP):
  490. term = IntegerGMP(term)
  491. _gmp.mpz_ior(result._mpz_p,
  492. self._mpz_p,
  493. term._mpz_p)
  494. return result
  495. def __rshift__(self, pos):
  496. result = IntegerGMP(0)
  497. if pos < 0:
  498. raise ValueError("negative shift count")
  499. if pos > 65536:
  500. if self < 0:
  501. return -1
  502. else:
  503. return 0
  504. _gmp.mpz_tdiv_q_2exp(result._mpz_p,
  505. self._mpz_p,
  506. c_ulong(int(pos)))
  507. return result
  508. def __irshift__(self, pos):
  509. if pos < 0:
  510. raise ValueError("negative shift count")
  511. if pos > 65536:
  512. if self < 0:
  513. return -1
  514. else:
  515. return 0
  516. _gmp.mpz_tdiv_q_2exp(self._mpz_p,
  517. self._mpz_p,
  518. c_ulong(int(pos)))
  519. return self
  520. def __lshift__(self, pos):
  521. result = IntegerGMP(0)
  522. if not 0 <= pos < 65536:
  523. raise ValueError("Incorrect shift count")
  524. _gmp.mpz_mul_2exp(result._mpz_p,
  525. self._mpz_p,
  526. c_ulong(int(pos)))
  527. return result
  528. def __ilshift__(self, pos):
  529. if not 0 <= pos < 65536:
  530. raise ValueError("Incorrect shift count")
  531. _gmp.mpz_mul_2exp(self._mpz_p,
  532. self._mpz_p,
  533. c_ulong(int(pos)))
  534. return self
  535. def get_bit(self, n):
  536. """Return True if the n-th bit is set to 1.
  537. Bit 0 is the least significant."""
  538. if self < 0:
  539. raise ValueError("no bit representation for negative values")
  540. if n < 0:
  541. raise ValueError("negative bit count")
  542. if n > 65536:
  543. return 0
  544. return bool(_gmp.mpz_tstbit(self._mpz_p,
  545. c_ulong(int(n))))
  546. # Extra
  547. def is_odd(self):
  548. return _gmp.mpz_tstbit(self._mpz_p, 0) == 1
  549. def is_even(self):
  550. return _gmp.mpz_tstbit(self._mpz_p, 0) == 0
  551. def size_in_bits(self):
  552. """Return the minimum number of bits that can encode the number."""
  553. if self < 0:
  554. raise ValueError("Conversion only valid for non-negative numbers")
  555. return _gmp.mpz_sizeinbase(self._mpz_p, 2)
  556. def size_in_bytes(self):
  557. """Return the minimum number of bytes that can encode the number."""
  558. return (self.size_in_bits() - 1) // 8 + 1
  559. def is_perfect_square(self):
  560. return _gmp.mpz_perfect_square_p(self._mpz_p) != 0
  561. def fail_if_divisible_by(self, small_prime):
  562. """Raise an exception if the small prime is a divisor."""
  563. if is_native_int(small_prime):
  564. if 0 < small_prime < 65536:
  565. if _gmp.mpz_divisible_ui_p(self._mpz_p,
  566. c_ulong(small_prime)):
  567. raise ValueError("The value is composite")
  568. return
  569. small_prime = IntegerGMP(small_prime)
  570. if _gmp.mpz_divisible_p(self._mpz_p,
  571. small_prime._mpz_p):
  572. raise ValueError("The value is composite")
  573. def multiply_accumulate(self, a, b):
  574. """Increment the number by the product of a and b."""
  575. if not isinstance(a, IntegerGMP):
  576. a = IntegerGMP(a)
  577. if is_native_int(b):
  578. if 0 < b < 65536:
  579. _gmp.mpz_addmul_ui(self._mpz_p,
  580. a._mpz_p,
  581. c_ulong(b))
  582. return self
  583. if -65535 < b < 0:
  584. _gmp.mpz_submul_ui(self._mpz_p,
  585. a._mpz_p,
  586. c_ulong(-b))
  587. return self
  588. b = IntegerGMP(b)
  589. _gmp.mpz_addmul(self._mpz_p,
  590. a._mpz_p,
  591. b._mpz_p)
  592. return self
  593. def set(self, source):
  594. """Set the Integer to have the given value"""
  595. if not isinstance(source, IntegerGMP):
  596. source = IntegerGMP(source)
  597. _gmp.mpz_set(self._mpz_p,
  598. source._mpz_p)
  599. return self
  600. def inplace_inverse(self, modulus):
  601. """Compute the inverse of this number in the ring of
  602. modulo integers.
  603. Raise an exception if no inverse exists.
  604. """
  605. if not isinstance(modulus, IntegerGMP):
  606. modulus = IntegerGMP(modulus)
  607. comp = _gmp.mpz_cmp(modulus._mpz_p,
  608. self._zero_mpz_p)
  609. if comp == 0:
  610. raise ZeroDivisionError("Modulus cannot be zero")
  611. if comp < 0:
  612. raise ValueError("Modulus must be positive")
  613. result = _gmp.mpz_invert(self._mpz_p,
  614. self._mpz_p,
  615. modulus._mpz_p)
  616. if not result:
  617. raise ValueError("No inverse value can be computed")
  618. return self
  619. def inverse(self, modulus):
  620. result = IntegerGMP(self)
  621. result.inplace_inverse(modulus)
  622. return result
  623. def gcd(self, term):
  624. """Compute the greatest common denominator between this
  625. number and another term."""
  626. result = IntegerGMP(0)
  627. if is_native_int(term):
  628. if 0 < term < 65535:
  629. _gmp.mpz_gcd_ui(result._mpz_p,
  630. self._mpz_p,
  631. c_ulong(term))
  632. return result
  633. term = IntegerGMP(term)
  634. _gmp.mpz_gcd(result._mpz_p, self._mpz_p, term._mpz_p)
  635. return result
  636. def lcm(self, term):
  637. """Compute the least common multiplier between this
  638. number and another term."""
  639. result = IntegerGMP(0)
  640. if not isinstance(term, IntegerGMP):
  641. term = IntegerGMP(term)
  642. _gmp.mpz_lcm(result._mpz_p, self._mpz_p, term._mpz_p)
  643. return result
  644. @staticmethod
  645. def jacobi_symbol(a, n):
  646. """Compute the Jacobi symbol"""
  647. if not isinstance(a, IntegerGMP):
  648. a = IntegerGMP(a)
  649. if not isinstance(n, IntegerGMP):
  650. n = IntegerGMP(n)
  651. if n <= 0 or n.is_even():
  652. raise ValueError("n must be positive odd for the Jacobi symbol")
  653. return _gmp.mpz_jacobi(a._mpz_p, n._mpz_p)
  654. @staticmethod
  655. def _mult_modulo_bytes(term1, term2, modulus):
  656. if not isinstance(term1, IntegerGMP):
  657. term1 = IntegerGMP(term1)
  658. if not isinstance(term2, IntegerGMP):
  659. term2 = IntegerGMP(term2)
  660. if not isinstance(modulus, IntegerGMP):
  661. modulus = IntegerGMP(modulus)
  662. if modulus < 0:
  663. raise ValueError("Modulus must be positive")
  664. if modulus == 0:
  665. raise ZeroDivisionError("Modulus cannot be zero")
  666. if (modulus & 1) == 0:
  667. raise ValueError("Odd modulus is required")
  668. product = (term1 * term2) % modulus
  669. return product.to_bytes(modulus.size_in_bytes())
  670. # Clean-up
  671. def __del__(self):
  672. try:
  673. if self._mpz_p is not None:
  674. if self._initialized:
  675. _gmp.mpz_clear(self._mpz_p)
  676. self._mpz_p = None
  677. except AttributeError:
  678. pass