python int VS gmpy2.mpz

python_int_VS_gmpy2mpz.py

from gmpy2 import mpz
from timeit import timeit
from random import randint
from time import perf_counter

# Здесь мы создаём логарифмическую шкалу длин
from math import log10
import numpy as np
int_lens = list(np.arange(1,100,1)) + list(np.logspace(2, log10(500000), num=150).astype(np.int32))
int_lens = [int(x) for x in int_lens]

# Длины чисел в 2**30 системе счисления, забиты константами, на случай, если нет numpy
# int_lens = \
# [1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49,
#  51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 71, 73, 75, 77, 79, 81, 83, 85, 87, 89, 91, 93, 95, 97,
#  99, 100, 103, 107, 111, 115, 119, 123, 128, 132, 137, 142, 147, 153, 158, 164, 170, 176, 183,
#  189, 196, 203, 211, 218, 226, 234, 243, 252, 261, 270, 280, 290, 301, 312, 323, 335, 347, 359,
#  372, 386, 400, 414, 429, 445, 461, 478, 495, 513, 531, 551, 571, 591, 613, 635, 658, 682, 706,
#  732, 759, 786, 815, 844, 875, 906, 939, 973, 1008, 1045, 1083, 1122, 1163, 1205, 1248, 1293,
#  1340, 1389, 1439, 1491, 1545, 1601, 1659, 1719, 1782, 1846, 1913, 1982, 2054, 2128, 2205,
#  2285, 2368, 2454, 2543, 2635, 2730, 2829, 2931, 3037, 3147, 3261, 3379, 3502, 3629, 3760,
#  3896, 4037, 4183, 4335, 4492, 4654, 4823, 4997, 5178, 5366, 5560, 5761, 5970, 6186, 6410,
#  6642, 6882, 7132, 7390, 7657, 7934, 8222, 8519, 8828, 9147, 9478, 9822, 10177, 10546, 10927,
#  11323, 11733, 12158, 12598, 13054, 13526, 14016, 14523, 15049, 15594, 16158, 16743, 17350,
#  17978, 18628, 19303, 20000]

TESTS = 2  # Сколько раз повторять самый тяжёлый тест
# В списке будет столько чисел, чтобы время всех эксперимернов было примерно одинаковым
reps = [round((max(int_lens) / i)**1.01 * TESTS) for i in int_lens]

# Генерим случайные числа
randsp1 = [[randint(1, 1<<(30*int_len)) for _ in range(rep)] for int_len, rep in zip(int_lens, reps)]
randsp2 = [[randint(1, 1<<(30*int_len)) for _ in range(rep)] for int_len, rep in zip(int_lens, reps)]

# Сюда мы будем записывать тайминги
py_int_times = []
mpz_int_times = []

# Тестируем умножение
for int_len, x1s, x2s in zip(int_lens, randsp1, randsp2):
    cur_dur = 0.0
    for a, b in zip(x1s, x2s):
        st = perf_counter()
        x = a * b
        en = perf_counter()
        cur_dur += en - st
    dur_pm = cur_dur / len(x1s)
    py_int_times.append(dur_pm)
    print('py:  {} bits, {:.3} spm, {:.3} tot'.format(30 * int_len, dur_pm, cur_dur))

# Тест тяжёлый, поэтому экономим память
randsm1 = [[mpz(x) for x in row] for row in randsp1]
del(randsp1)
randsm2 = [[mpz(x) for x in row] for row in randsp2]
del(randsp2)

# то же самое для mpz
for int_len, x1s, x2s in zip(int_lens, randsm1, randsm2):
    cur_dur = 0.0
    for a, b in zip(x1s, x2s):
        st = perf_counter()
        x = a * b
        en = perf_counter()
        cur_dur += en - st
    dur_pm = cur_dur / len(x1s)
    mpz_int_times.append(dur_pm)
    print('mpz: {} bits, {:.3} spm, {:.3} tot'.format(30 * int_len, dur_pm, cur_dur))

# Экономим память
del(randsm1)
del(randsm2)




# Теперь начинаем рисовать графики
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rc('font', family='Verdana')  # Для русского языка

# Просто сами графики
plt.plot(int_lens, py_int_times, linewidth=3.0, label="python int", color='g')
plt.plot(int_lens, mpz_int_times, linewidth=3.0, label="gmpy2.mpz int", color='r')

# Теперь ищем асимптоты на лог-графике методом наименьших квадратов
from scipy.stats import linregress
import numpy as np
ST_F = 100
slope, intercept, _, _, _ = linregress(np.log(np.array(int_lens[ST_F:])), np.log(np.array(py_int_times[ST_F:])))
py_mul, py_pow = np.exp(intercept), slope
plt.plot(int_lens, np.exp(intercept) * np.power(int_lens, slope), linewidth=2.0, linestyle='--', label="python asmp", color='m')

slope, intercept, _, _, _ = linregress(np.log(np.array(int_lens[ST_F:])), np.log(np.array(mpz_int_times[ST_F:])))
mpz_mul, mpz_pow = np.exp(intercept), slope
plt.plot(int_lens, np.exp(intercept) * np.power(int_lens, slope), linewidth=2.0, linestyle='--', label="gmpy2.mpz asmp", color='b')


# На самом деле у mpz асимптотика n log(n) log(log(n)). Ищем подходящую константу С
xx = int_lens[ST_F:]
yy = mpz_int_times[ST_F:]
FFT = np.array(int_lens[ST_F:]) * np.log(int_lens[ST_F:]) * np.log(np.log(int_lens[ST_F:]))
lft, rht = np.min(yy/FFT), np.max(yy/FFT)
dv1000 = np.linspace(lft, rht, 1000)
pen = [np.sum((yy - FFT*x)**2) for x in dv1000]
best_C = dv1000[np.argmin(pen)]
# Нашли константу, рисуем асимптотику на графике
FFT_plt = best_C * np.array(int_lens[2:]) * np.log(int_lens[2:]) * np.log(np.log(int_lens[2:]))
plt.plot(int_lens[2:], FFT_plt, linewidth=2.0, linestyle='--', label="FFT asmp", color='y')


# Подписываем функции
plt.annotate(r'$({:.03f})\cdot10^{{-9}}\cdot n^{{{:.05}}}$'.format(py_mul*10**9, py_pow).replace('.', '{.}'),
             xy=(int_lens[-30], py_int_times[-30]), xytext=(int_lens[50], py_int_times[-25]),
             arrowprops=dict(facecolor='m', shrink=0.05))
plt.annotate(r'$({:.03f})\cdot10^{{-9}}\cdot n^{{{:.05}}}$'.format(mpz_mul*10**9, mpz_pow).replace('.', '{.}'),
             xy=(int_lens[-30], mpz_int_times[-30]), xytext=(int_lens[150], mpz_int_times[50]),
             arrowprops=dict(facecolor='b', shrink=0.05))

plt.annotate(r'$({:.03f})\cdot10^{{-9}}\cdot n \cdot \log n \cdot \log(\log n)$'.format(best_C*10**9).replace('.', '{.}'),
             xy=(int_lens[-20], mpz_int_times[-20]), xytext=(int_lens[-50], mpz_int_times[5]),
             arrowprops=dict(facecolor='y', shrink=0.05))


# Легенды и красивости
plt.legend(loc=4)
plt.xscale('log')
plt.yscale('log')
plt.title('python int VS gmpy2.mpz')
plt.xlabel('Количество $2^{30}$-ричных цифр')
plt.ylabel('Время одного умножения, c')
plt.grid(True)

# Поехали! :)
plt.show()