Using GMP (gmpy2) in Python

在Python中使用GMP（gmpy2） 苏剑林 · 2014-10-28 · Mathematics number theory Python

I previously wrote A First Try at Using PARI/GP in Python, giving a brief introduction to calling PARI/GP from Python. PARI/GP is a fairly powerful number theory library, "designed for fast computations in number theory (factorization of large numbers, algebraic number theory, elliptic curves…)." It can be called from programming languages such as C/C++ or Python, and it is also a self-contained scripting language in its own right. However, if all you need is high-precision arithmetic on large numbers, GMP seems to fit our needs even better.

Readers familiar with C/C++ will know GMP (short for GNU Multiple Precision Arithmetic Library), an open-source library for high-precision arithmetic. It provides not only high-precision operations on ordinary integers, reals, and floating-point numbers, but also random number generation, and in particular a very comprehensive set of interfaces for number-theoretic operations — for example, the Miller-Rabin primality test, large prime generation, the Euclidean algorithm, computing the inverse of an element in a field, the Jacobi symbol, the Legendre symbol, and so on [source]. Although calling GMP from C/C++ isn't particularly complicated, being able to use GMP from Python — a language famous for its high development efficiency — would undoubtedly be a delightful thing. This is exactly what this post is about: gmpy2. more

Introduction to gmpy2

gmpy2 is a Python extension library that wraps GMP; its predecessor was gmpy. Through the author's refinements and wrapping, using gmpy2 has become far simpler. Based on my experience so far, using gmpy2 in Python, compared to using GMP directly in C/C++, has at least the following advantages:

1. Simplified function names: the author has greatly simplified many of GMP's function names. For example, for the probabilistic primality test, the command in C/C++ is mpz_probab_prime_p, while in gmpy2 it's simply is_prime.
2. Convenient operator overloading: in C/C++, to add two mpz integers, you'd need mpz_add(z_i, z_i, z_o); in gmpy2, you just use the + operator. More generally, adding an mpz integer to a regular int also just requires +, and the same goes for subtraction, multiplication, division, modulo, exponentiation, and so on.

My own knowledge is still shallow, so this is only a one-sided assessment. Of course, as a third-party language binding, gmpy2's efficiency is always a bit lower than calling GMP directly from C/C++, but the difference is small, since gmpy2 is essentially just a precompiled C library. For a detailed gmpy2 tutorial, see:

http://gmpy2.readthedocs.org/en/latest/

Basic Usage

This post only gives a brief introduction. All the code below was run in Python 3.4. To initialize a large integer, you simply need

import gmpy2
n=gmpy2.mpz(1257787) #初始化
gmpy2.is_prime(n) #概率性素性测试

This is parallel to C/C++; the argument inside the parentheses can be either an integer or a string. gmpy2 integrates not only the large-integer type mpz, but also the high-precision floating-point type mpfr, whose usage is likewise parallel to C/C++. Below are some basic computations

a+2 #求和，结果是一个mpz
a-2 #求差，结果是一个mpz
a*2 #求积，结果是一个mpz
a/2 #求商，结果是一个mpfr
a//2 #求商，结果是一个mpz
a**2 #求平方，结果是一个mpz
a%2 #求模，结果是一个mpz

The overloading of these operators is both intuitive and convenient. You can also replace the number 2 with an mpz-type variable, which of course goes without saying. In fact, the expression a+2 first converts 2 into the mpz value 2, and then performs addition within mpz.

Comparison with Python's Built-In Support

Python itself supports high-precision arithmetic on large integers, which is sufficient for a general, relatively small range of numbers — it's just not "fast." Readers only need to run the following two lines of code separately

gmpy2.mpz(1257787)**123456 #gmpy2计算
1257787**123456 #Python自带计算

to get a feel for it (of course, if your machine happens to be quite powerful, the two might show no noticeable difference — in that case, just scale up the exponent by a factor of ten~ My own laptop isn't very powerful, I'm afraid ^_^).

Trying Out Large-Number Factorization

Here's a piece of code I wrote myself for factoring large numbers:

from gmpy2 import *
import time
start=time.clock()
n = mpz(63281217910257742583918406571)
x = mpz(2)
y = x**2 + 1
for i in range(n):
    p = gcd(y-x,n)
    if p != 1:
        print(p)
        break
    else:
        y=(((y**2+1)%n)**2+1)%n
        x=(x**2+1)%n
end=time.clock()
print(end-start)

This code calls gmpy2, and the algorithm used is Pollard's rho method — run in its most bare-bones form, without any optimization. On my laptop, this algorithm factored

63281217910257742583918406571 = 125778791843321 × 503115167373251

in 84 seconds (stopping as soon as one factor was found). Of course, this is nothing impressive — Mathematica and PARI/GP handle this number instantaneously. This script was purely an exercise to practice and test calling gmpy2 from Python.

The steps of Pollard's rho method:

1. Given a composite number $n$, and $x_0=2,y_0=x_0^2+1$;
2. Compute $p=Gcd(x_0-y_0,n)$; if p is not 1, then the result is a factor of $n$, and we stop;
3. If p = 1, then update according to $x_{n+1}=x_n^2+1,y_{n+1}=(y_n^2+1)^2+1$, and repeat step 2.

Here, to reduce the amount of computation, each step's computation of $x_n ,y_n$ is reduced by taking it modulo $n$. The steps above are only the simplest part of the computation and don't yet handle the various exceptions that may arise, such as the case where $n|(x_n -y_n)$ occurs; nor have certain details been optimized — please cite with caution.

Reference URL:

http://hi.baidu.com/pytvzcnbuolrsue/item/ae714592f1fda8d87b7f010a