Непонятная ошибка при использовании факториала с numpy или массивов scipy

Question

Следующий код идентифицирует список простых чисел

import math
import numpy as np
import pandas as pd
import scipy
for x in range(20, 50):
    print(x, math.factorial(x-1) % x)

на выходе получается

20 0
21 0
22 0
23 22
24 0
25 0
26 0
27 0
28 0
29 28
30 0
31 30
32 0
33 0
34 0
35 0
36 0
37 36
38 0
39 0
40 0
41 40
42 0
43 42
44 0
45 0
46 0
47 46
48 0
49 0

Остаток вычисления не равен нулю для каждого простого числа. Когда я пытаюсь сделать то же самое вычисление с массивами, мои результаты отличаются.

arr = np.arange(20,50)
modFactArr = factorial(arr-1) % arr
print(np.column_stack((arr,modFactArr)),'\n\n')

smodFactArr=scipy.special.factorial(arr-1) % arr
print(np.column_stack((arr,smodFactArr)),'\n\n')

дает

[[20.  0.]
 [21.  0.]
 [22.  0.]
 [23. 22.]
 [24.  0.]
 [25. 22.]
 [26. 16.]
 [27. 11.]
 [28. 12.]
 [29. 16.]
 [30. 24.]
 [31.  8.]
 [32.  0.]
 [33. 18.]
 [34. 16.]
 [35. 33.]
 [36. 20.]
 [37. 12.]
 [38. 20.]
 [39. 10.]
 [40.  0.]
 [41. 25.]
 [42. 26.]
 [43.  6.]
 [44.  4.]
 [45.  3.]
 [46. 36.]
 [47. 40.]
 [48.  0.]
 [49. 12.]] 


[[20.  0.]
 [21.  0.]
 [22.  0.]
 [23. 22.]
 [24.  0.]
 [25. 22.]
 [26. 16.]
 [27. 11.]
 [28. 12.]
 [29. 16.]
 [30. 24.]
 [31.  8.]
 [32.  0.]
 [33. 18.]
 [34. 16.]
 [35. 33.]
 [36. 20.]
 [37. 12.]
 [38. 20.]
 [39. 10.]
 [40.  0.]
 [41. 25.]
 [42. 26.]
 [43.  6.]
 [44.  4.]
 [45.  3.]
 [46. 36.]
 [47. 40.]
 [48.  0.]
 [49. 12.]] 

Обратите внимание, что теперь такие числа, как 26,27,28 и др. c дают остатки сейчас. Это ошибка с моим кодом? или есть причина, по которой Сципи и numpy делают по модулю арифмети c по-разному?

Answer 1

Это происходит из-за np.dtype. У вас просто переполнение np.int64, в то время как нативное python int не может быть переполнено.

Answer 2

Из документов:

При exact=False факториал равен приблизительно с использованием гамма-функции

exact=False по умолчанию.

Вы можете сказать, что это аппроксимируется, потому что результат - это число с плавающей точкой (отсюда .0), и число с плавающей точкой не может точно хранить интегральные результаты за пределами 2 ** 53.