應數一 整數論 以 Python 實驗 之講義草稿1

對象: 高雄師範大學數學系大一應用數學組同學
其預備知識之假設: 1. Python 沒學過, 2. 懂基礎之高中數學

本課程內容:

1. Pyhton 求因數: 1.1 用最直接的方法 試除法 求因數 1.2 加速 試除法 求因數,
2. Pyhton 求質數: 2.1 Pyhton 求質數最基本(慢)的設計
3. Python 實作 輾轉相除法，及理論
4. 丟方圖方程、百雞問題: 張丘建算經6世紀
5. 展現 Python 呈現超長整數的能力: 孿生質數猜想: 用 Python 呈現 2009 為止最大的孿生質數 Twin Prime conjecture (Using Python to show the biggest twin prime pairs up to 2009)
6. Fermat’s Little Theorem 費馬小定理
7. Euler’s Theorem - the generalization of Fermat’s Little Theorem 歐拉定理 是費馬小定理的推廣
8. 凱薩密碼
9. 凱薩密碼破解
10. RSA密碼

Ref:

• Burton, Elementary Number Theory, McGraw Hill 2011.
• AI Sweigart, Cracking Codes With Python
• 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 人民郵電, 2020.

---

---

0. Python 入門

Python 的安裝使用, 請參考作者之系列文章:
從turtle海龜動畫 學習 Python - 高中彈性課程系列 2 安裝 Python, 線上執行 Python link

1. Python 求因數

註: 比較程式執行時間請看後面 附錄 Appendix

1.1 用最直接的方法 試除法 求因數 Using Python to explore Number Theory — output factors of a given number

Using the most basic method (try by dividing) to output the factors of an input integer number
with Pyhton.
By Prof. P-J Lai MATH NKNU 201403

以下是最簡單直接的寫法, 用試除法, 輸入n, 印出 n 的所有因數
註: codes 存為 FactorList.py, 按F5執行後,函數名例如factorList(10)才會有定義,
執行動作時, 是呼叫函數名, 例如 factorList(20)

def Factor(n):
List=[]
   for i in range(1,n+1):
        if n % i==0:
            List.append(i)
   for i in List:
        print(i)

print List 內容, 以下亦可, 但是較冗長

for i in range(len(List)):
    print(List[i])

執行畫面:

>>> Factor(100)
1
2
4
5
10
20
25
50
100

Ex: 將以上factorList(n) 函數改成, 用 return 關鍵字, 傳回 factors之list
Ans: 用 return(List)

def Factor(n):
    '''這裡是本子程式的說明

           此程式Factor(n) 會回應一個List含有所有n的因數'''
    
    List=[]
    for i in range(1,n+1):
        if n % i==0:
            List.append(i)

    return(List)

Ex: 將以上factorList(n) 函數改成, print list 時, 不跳行, 每個因數用 空白 間隔
Ans: 用 print(i, end=’ ‘)

from math import *
#import math
def factorList(n):
    ‘’’這裡是本子程式的說明
    此程式Factor(n) 會列印出所有n的因數在螢幕上’’’
    List=[]
    for i in range(1, n+1):
        if n % i==0:
            List.append(i)
    for i in List:
        print(i, end=’ ‘)

1.2 加速 試除法 求因數 Using Python to explore Number Theory 2 — output factors of a given number```python

將以上factorList(n) 函數改成, 因數只檢查到 n/2, print list 時, 不跳行, 每個因數用 空白 間隔,
當 n 很大時, 只檢測到 n//2 +1, 會大大加速!
(更快速的加速可以看下一節: 1.3 最簡短的 Python codes 算因數)

注意此時列出之因數不會有 n 這個數, 例如以下執行之圖片列出 12000之 因數, 只會列出 到 6000 以內的所有因數. 可以在最後加上
fooList.append(n)
手動將 n 加入因數 fooList 中

以下實做時
將試除法之 for loop 改成只檢查到 n/2 附近, 要小心 n 有可能是奇數:
for i in range(1, n//2+1)

# 20230316_因數分解_加速版本試除到一半_fun_return_half_n_4.py.py
# math.floor(n/2) 等價 n//2

def fooFactorHalfN(n):
    fooList = []
    for i in range(1, n//2+1):
        if n % i == 0:
            fooList.append(i)
    fooList.append(n)
    return fooList
    # return(fooList) also OK.


##>>> 
##===== RESTART: C:/Users/user/Desktop/20230316_因數分解_加速版本試除到一半_fun_return_half_n_4.py =====
##>>> fooFactorHalfN(12)
##[1, 2, 3, 4, 6, 12]
##>>> fooFactorHalfN(12001)
##[1, 11, 1091, 12001]
##>>> fooFactorHalfN(1091)
##[1, 1091]

執行結果

>>> 
=== RESTART: C:/Users/user/Desktop/20230316_因數分解_加速版本試除到一半_fun_return_half_n_4.py ==
>>> fooFactorHalfN(120)
[1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 24, 30, 40, 60, 120]

…

1.3 最簡短的列出因數的 Python codes (Matlab, R) The very short Python codes to compute factors of a given integer number N (also Matlab, R)

用最精簡的程式碼列出因數
請參考本人 CSDN 另一篇:
https://blog.csdn.net/m0_47985483/article/details/106241917)
link

…

1.4 測量以上各種不同求因數之寫法, 使用 time.time(), 比較它們的效能

time.time(), 早期是 time.clock(), 目前已棄用 deprecated!

Ex: 使用 time.time() 測量以上各種不同求因數之寫法, 比較它們的效能.

# 20230316_因數分解_最簡版本3與HalfN4版本時間比較_time.time_format_5.1.py
# By Prof. P-J Lai MATH NKNU 20230317
# time.time(), 早期是 time.clock(), 目前已棄用 deprecated!
# 使用 f.'字串'之輸出方式


# 20230316_因數分解_最簡版本_List.append_fun_return_3.py.py
def fooFactor(n):
    fooList = []
    for i in range(1,n+1):
        if n % i == 0:
            fooList.append(i)
    return fooList


# 20230316_因數分解_加速版本試除到一半_fun_return_half_n_4.py.py
# math.floor(n/2) 等價 n//2
def fooFactorHalfN(n):
    fooList = []
    for i in range(1, n//2+1):
        if n % i == 0:
            fooList.append(i)
    fooList.append(n)
    return fooList


import time

n=12000000

startTime = time.time()
fooFactor(n)
endTime = time.time()
print(f'The excution time of fooFactor({n})',  endTime - startTime )

startTime = time.time()
fooFactorHalfN(n)
endTime = time.time()
print(f'The excution time of fooFactorHalfN({n})', endTime - startTime )



# 執行結果
##>>> 
##= RESTART: E:/NEW_筆電的/網路免費軟體資料/Python教學/Python數論/用Python列出因數_List要放在函數內/20230316_因數分解_最簡版本3與HalfN4版本時間比較_time.time_format_5.1.py
##The excution time of fooFactor(12000000) 0.7441365718841553
##The excution time of fooFactor(12000000) 0.37409281730651855

1.5 測量以上各種不同求因數之寫法, 使用 timeit.timeit, 比較它們的效能

Ex: 使用 timeit.timeit 測量以上各種不同求因數之寫法, 比較它們的效能.
Ans:

# 20230316_因數分解_最簡版本3與HalfN4版本時間比較_5.py
# By Prof. P-J Lai MATH NKNU 20230317
##另一種測量時間之語法:
# import timeit
# timeit.timeit()
##針對時間很短暫的程式之測量用
##預設是執行1000000，一百萬次
# 使用 f.'字串'之輸出方式


# 20230316_因數分解_最簡版本_List.append_fun_return_3.py.py
def fooFactor(n):
    fooList = []
    for i in range(1,n+1):
        if n % i == 0:
            fooList.append(i)
    return fooList


# 20230316_因數分解_加速版本試除到一半_fun_return_half_n_4.py.py
# math.floor(n/2) 等價 n//2
def fooFactorHalfN(n):
    fooList = []
    for i in range(1, n//2+1):
        if n % i == 0:
            fooList.append(i)
    fooList.append(n)
    return fooList


import timeit
n=120
print(f'fooFactor({n})', timeit.timeit(f'fooFactor({n})', setup='from __main__ import fooFactor'))
print(f'fooFactorHalfN({n})', timeit.timeit(f'fooFactorHalfN({n})', setup='from __main__ import fooFactorHalfN'))
# 執行結果
##>>> 
##= RESTART: E:\NEW_筆電的\網路免費軟體資料\Python教學\Python數論\用Python列出因數_List要放在函數內\20230316_因數分解_最簡版本3與HalfN4版本時間比較_timeit_6.py
##fooFactor(120) 6.5922105
##fooFactorHalfN(120) 3.7491107



### 如果少了 setup=',,,設定, 會出問題!
##from timeit import timeit
##print(timeit('fooFactor(120)'))
##print(timeit('fooFactorHalfN(120)'))
##>>> 
##======= RESTART: C:/Users/user/Desktop/20230316_因數分解_最簡版本3與HalfN4版本時間比較_5.py =======
##Traceback (most recent call last):
##  File "C:/Users/user/Desktop/20230316_因數分解_最簡版本3與HalfN4版本時間比較_5.py", line 34, in <module>
##    print(timeit('fooFactor(120)'))
##  File "C:\Users\user\AppData\Local\Programs\Python\Python39\lib\timeit.py", line 233, in timeit
##    return Timer(stmt, setup, timer, globals).timeit(number)
##  File "C:\Users\user\AppData\Local\Programs\Python\Python39\lib\timeit.py", line 177, in timeit
##    timing = self.inner(it, self.timer)
##  File "<timeit-src>", line 6, in inner
##NameError: name 'fooFactor' is not defined




# 使用__name__=='__main__' 
##if __name__=='__main__':
##    from timeit import timeit
##    print(timeit('fooFactor(120)',  setup='from __main__ import fooFactor'))
##    print(timeit('fooFactorHalfN(120)',  setup='from __main__ import fooFactorHalfN'))

#執行的次數，預設是1000000，一百萬
## n =12 時 
##======= RESTART: C:/Users/user/Desktop/20230316_因數分解_最簡版本3與HalfN4版本時間比較_5.py =======
##1.7224516999999997
##1.1633394
##
##
##
## n =120 時 
##>>> 
##======= RESTART: C:/Users/user/Desktop/20230316_因數分解_最簡版本3與HalfN4版本時間比較_5.py =======
##6.7437
##4.3714336

2. Python 求質數

2.1 Python 求質數最基本(慢)的設計

輸入 n, 列出 n 以內的 所有質數, 分三個函數執行:

       程式 Factor(n) 會回應 (return) 一個 List 含有所有 n 的因數
       程式 IsPrime(n), 對輸入的 n,  會回應一個 boole value, 表示 n 是否為質數
       程式 Primelist(n)  會列出所有小於等於 n 的質數'''

這樣的設計, 速度會較慢, 但是程式碼較清楚, 之後會討論加速的設計

# By Prof. P-J Lai MATH NKNU 201407, 20230323 revised
# 20230323_質數_print方式_最簡版本_複雜度分析.py
# collect all function in a single file

'''這裡是本 script file 的說明

           程式 Factor(n) 會回應 (return) 一個 List 含有所有 n 的因數
           程式 IsPrime(n), 對輸入的n,  會回應一個 boole value, 表示 n 是否為質數
           程式 Primelist(n)  會列出所有小於等於 n 的質數'''


def fooFactor(n):
    '''這裡是本子程式的說明

           此程式Factor(n) 會回應一個List含有所有n的因數'''
    
    List=[]
    for i in range(1,n+1):
        if n % i==0:
            List.append(i)

    return(List)



def fooIsPrime(n):
    if len(fooFactor(n))==2:
        return(True)
    else:
        return(False)


# 不同的輸出呈現, 來觀察 fooIsPrime(n) 到底在做甚麼事:
##>>> for i in range(1,6):
##	fooIsPrime(i)
##
##	
##False
##True
##True
##False
##True


##>>> for i in range(1,6):
##	print("Is i prime?",fooIsPrime(i))
##
##	
##Is i prime? False
##Is i prime? True
##Is i prime? True
##Is i prime? False
##Is i prime? True


##>>> for i in range(1,6):
##	print(f"Is {i} prime?",fooIsPrime(i))
##
##	
##Is 1 prime? False
##Is 2 prime? True
##Is 3 prime? True
##Is 4 prime? False
##Is 5 prime? True



##n=6
##for i in range(1, n):
##    if fooIsPrime(i):
##        print(f"Is {i} prime?", " Yes!")
##    else:
##        print(f"Is {i} prime?", " No!")
##
##>>> 
##==== RESTART: C:/Users/user/Desktop/test.py ===
##Is 1 prime?  No!
##Is 2 prime?  Yes!
##Is 3 prime?  Yes!
##Is 4 prime?  No!
##Is 5 prime?  Yes!


def fooPrimeList(n):
    List=[]
    for i in range(1,n+1):
        if fooIsPrime(i):
           List.append(i)

    return(List)


# Complexity analysis 計算複雜度分析
##O(n)
##O(n^3)
##O(e^n)
##O(2^n)
##O(ln(n))
##
##
##O(fooFactor(n))= O(n)=K*n
##O(fooIsPrime(n))= O(fooFactor(n))= O(n)
##
##O(fooPrimeList(n))=range(1,n+1)*O(fooIsPrime(i))
##=n*O(fooIsPrime(i))=n*O(n) =O(n^2)

2.2 Python 求質數最簡加速 試除法只檢測到 $\frac{n}{2}$ 使用 break 或 return 返回

如果 n 是合成數, n 的因數一定會有至少一個是介於 1 到 $\frac{n}{2}$,
也就是:
Lemma:
假設 n 是合成數, 至少存在一個 因數 k 使得
$$1<k\le \frac{n}{2}.$$
基於這樣的觀察, 我們就可以改進之前的求質數最基本(慢)的設計, 對他做加速,
程式碼試除法只檢測到 $\frac{n}{2}$ :
也就是 n 依序除以(檢測) 1~ n//2,
觀察
$$n\%2,n\%3,n\%4,\cdots ,n\%(n//2)$$

一旦偵測到有真因數, 就使用 break 跳出迴圈,
或是用 return 直接結束函數之執行.

也就是檢測 1 到 $\frac{n}{2}$,一偵測到有真因數,就 return 結束.
因為 $\frac{n}{2}$ 有可能有小數位數, 所以 可以用Python 的 整除 // 算子.

注意底下迴圈需使用: for i in range(1,n//2+1)
而非 for i in range(1,n//2)

# By Prof. P-J Lai MATH NKNU 201407
## 202303023 revised
##如果 n 是合成數,  n 的因數一定會有至少一個是介於 1 到 $\frac{n}{2}$, 
##也就是:
##**Lemma:**
##假設  n 是合成數, 至少存在一個 因數 k 使得
## $$1 < k \;\le \; \frac{n}{2}$$.
## 
## 基於這樣的觀察, 我們就可以改進之前的求質數最基本(慢)的設計, 對他做加速, 
##程式碼只檢測到 $\frac{n}{2}$ : n 除以(檢測) `1~ n//2+1`,
##   也就是檢測 1 到 $\frac{n}{2}$,一偵測到有真因數,就 `return` 結束.
##   因為 $\frac{n}{2}$ 有可能有小數位數, 所以 可以用Python 的 整除 `//` 算子

'''這裡是本 script file 的說明

        程式 FactorHalf(n) 會回應一個 List 含有所有n的因數
        程式 IsPrimeHalf(n)會回應一個 boole value 表示 n 是否為質數
        程式 PrimeListHalf  會列出所有小於等於n的質數'''


import time

# Be careful not to use: for i in range(1,n//2),  should be range(1,n//2+1):
def FactorHalf(n):
    '''這裡是本子程式的說明
    此程式Factor(n) 會回應一個 List 含有所有 n 的因數'''
    
    List=[]
    for i in range(1,n//2+1):
        if n % i==0:
            List.append(i)
    List.append(n)
    
    return(List)


def IsPrimeHalf(n):
    if len(FactorHalf(n))==2:
        return(True)
    else:
        return(False)


def PrimeListHalf(n):
    List=[]
    for i in range(1,n+1):
        if IsPrimeHalf(i):
           List.append(i)

    return(List)


##>>> PrimeListHalf(12000)
##[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
## 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103,
## 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
## ,,,,,,,,,,,
## 11939, 11941, 11953, 11959, 11969, 11971, 11981, 11987]

2.3 Python 求質數簡單加速 試除法只檢測到 $\sqrt{n}$ 使用 break 或 return 返回

進一步觀察,
如果 n 是合成數, n 的因數一定會有至少一個是介於 1 到 $\sqrt{n}$,
也就是:
Lemma:
假設 n 是合成數, 至少存在一個 因數 k 使得
$$1<k\le \sqrt{n}.$$

故, 我們可以進一步改進效能,
因數只檢測到 $\sqrt{n}$,
一旦偵測到有真因數, 就使用 break 跳出迴圈,
或是用 return 直接結束函數之執行.

# optimized codes 2
# PrimeList_sqrt.py

def IsPrime_sqrt(n):
    '''這裡是本子程式的說明
    列出質數加速版本之一
    改進效能, 因數只檢測到 $\sqrt{n}$, 
    一旦偵測到有真因數, 使用 `break` 跳出迴圈, 
    或是用 `return` 直接結束函數之執行.
    另一個是 檢測到 n/2: 除以(檢測)1~ n//2+1,
    也就是檢測1~ n/2,一偵測到有真因數,就return 結束.
    By Prof. P-J Lai MATH NKNU 20220302
    '''
    
    List=[]
    for i in range(2,math.ceil(math.sqrt(n+1))):
        if n % i==0:
           return False
    return True
    

def PrimeList_sqrt(n):
    List=[]
    for i in range(2,n+1):
        if IsPrime_sqrt(i):
           List.append(i)

    return List, len(List)

# return(List,count) OK

Ex: 檢測 試除到 $\frac{n}{2}$ 與 試除到 $\sqrt[2]{n}$ 之時間加速比較.

並比較 試除到 $\frac{n}{2}$ 與 試除到 $\sqrt[2]{n}$ 之時間
程式碼如下

# 試除到 n/2, sqrt(n), 比較時間
import time
import math
#########################################
# optimized codes 1

def IsPrime_half(n):
    '''這裡是本子程式的說明
       此程式Factor(n) 會回應一個List含有所有n的因數'''
    List=[]
    ##flag=True
    for i in range(2,math.ceil((n+1)/2)):
        if n % i==0:
##           flag=False
##           break
            return( False)
    return True

def PrimeList_half(n):
    List=[]
    for i in range(2,n+1):
        if IsPrime_half(i):
           List.append(i)

    return List,len(List)

# return(List,count) OK

#########################################
# optimized codes 2

def IsPrime_sqrt(n):
    '''這裡是本子程式的說明
       此程式Factor(n) 會回應一個List含有所有n的因數'''
    List=[]
    for i in range(2,math.ceil(math.sqrt(n+1))):
        if n % i==0:
           return False
    return True
    

def PrimeList_sqrt(n):
    List=[]
    for i in range(2,n+1):
        if IsPrime_sqrt(i):
           List.append(i)

    return List,len(List)

# return(List,count) OK


# 試除到 n/2, sqrt(n), 比較時間

start = time.time()
PrimeList_half(100)
print('PrimeList_half(100)', time.time() - start)

start = time.time()
PrimeList_sqrt(100)
print('PrimeList_sqrt(100)', time.time() - start)

# 練習用timeit 模組量時間

from timeit import timeit
print('PrimeList_half(100)', timeit('PrimeList_half(100)',  setup='from __main__ import PrimeList_half', number=3))
print('PrimeList_sqrt(100)',timeit('PrimeList_sqrt(100)', setup='from __main__ import PrimeList_sqrt', number=3))


##>>> 
##= RESTART: D:\NEW_筆電的\網路免費軟體資料\Python教學\Python數論\用Python列出質數_加速版本\PrimeListAllCodes_optimize1.py
##PrimeList_half(100) 0.0
##PrimeList_sqrt(100) 0.0
##PrimeList_half(100) 0.0011965000000000447
##PrimeList_sqrt(100) 0.0002674999999999761

2.4 測量以上不同求求質數之寫法, 使用 time.time, 及 timeit.timeit, 比較它們的效能

2.5 Sieve of Eratosthenes 質數篩去法

理論介紹可以參考: 李永樂老師, 如何快速篩選質數？費馬素性檢驗和米勒-拉賓測試, https://youtu.be/E5vtn_OIe-c link

程式碼可以觀摩
SciPy 官網說明手冊的版本
1_0_prime_sieve.py,
此版本主要是用 NumPy 的 array,
如果不使用 NumPy 的寫法, 可以參考後面作者改寫河西朝雄的 C codes 為Python.

2.5.0 最直觀但是效能較低的版本: 使用 list.remove() 函數

最直觀但是效能較低的版本
直接以Python list 的 .remove()函數, 將偵測到的合成數直接 remove 掉, 會較直觀, 但是是效能較低的寫法,

Remark: 一般C語言的寫法是在下標為 k 之處, 存放0, 1 來表示 k 是質數或合成數.

# Eratosthenes_remove_0.py
# 待除錯,  會出現 4
#primeList=[ k for k in range(2,N+1)] 是質數會出現 4的原因
# 已找到錯誤處 20230330

import time

def Eratosthenes_remove_0(N):
    '''Find the primes in 2~N'''

    primeList=[ k for k in range(0,N+1)] 

    for i in primeList[2:]:
        for j in primeList[i+1:]:
            if j%i==0:
                primeList.remove(j)

    return primeList[2:]

print(Eratosthenes_remove_0(100))          
##>>> 
##= RESTART: E:\NEW_筆電的\網路免費軟體資料\Python教學\Python數論\用Python列出質數_加速版本\Eratosthenes厄爾多拉氏篩法用python_Burton3.2\Eratosthenes\Eratosthenes_remove_0.py
##[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

2.5.0.1 直觀但是效能稍快的版本: 使用 下標為 k 之處, 存放0, 1 來表示 k 是質數或合成數.

# Eratosthenes_simple_Lai.py
# 此版本沒有做最佳化加速, 以最直觀簡單的方式寫!
# 不使用Python list 的 `.remove()`函數
# 使用一般C語言的寫法是在下標為 k 之處, 存放 0, 1 來表示 k 是質數或合成數.

import math
import time

def Eratosthenes(N):
    Prime=[]
    primeList=[]
    for i in range(0,N+1):
        Prime.append(1)
      
    for i in range(2,N+1):
        if Prime[i]==1:
            for j in range(i+1,N+1):
                if j%i==0:
                    Prime[j]=0

    for i in range(2,N+1):
        if Prime[i]==1:
            primeList.append(i)

    return primeList


print(Eratosthenes(100))

##>>> 
##= RESTART: E:/NEW_筆電的/網路免費軟體資料/Python教學/Python數論/用Python列出質數_加速版本/Eratosthenes厄爾多拉氏篩法用python_Burton3.2/Eratosthenes/Eratosthenes_simple_Lai.py
##[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

2.5.1 SciPy 官網說明手冊的版本是最佳化過的

SciPy 官網說明手冊的版本的寫法是最佳化過的!

Problem: 注意加速的部分是,
外部 loop 只跑到 ceil(sqrt(n))+1,
內部 loop 則自 i*i 開始跑,
跑的時候是間格 i,
(numPy的程式碼是 mask[j*j::j] = False, 河西的版本沒使用間格 i, SciPy的有 )
請同學自己證明一下為何可以如此?

注意原版本是 Python2 有 xrange(), Python3 要改成 range() 即可.

"""
Computing prime numbers with the archimedean sieve.

(Of course, this is not an optimal way for computing prime numbers...)

1_0_prime_sieve.py
"""

# Pyhton 的步長是放在最後, Matlab則放中間
# mask[j*j::j], step 是i
# mask = np.ones([N], dtype=bool)效果跟
# mask = np.ones(N, dtype=bool) 一樣

##SciPy 官網說明手冊的版本
##1_0_prime_sieve.py
##注意原版本是 Python2 有 xrange(), Python3 要改成 range() 即可.
##By Prof. P-J Lai MATH NKNU 20220302

import numpy as np

eratosthenes = True

# maximum number
N = 1000

# mask for prime numbers
mask = np.ones([N], dtype=bool)

if not eratosthenes:
    # simple prime sieve
    mask[:2] = False
    #for j in xrange(2, int(np.sqrt(N)) + 1):
    for j in range(2, int(np.sqrt(N)) + 1):
        mask[j*j::j] = False

else:
    # Eratosthenes sieve
    mask[:2] = False
    #for j in xrange(2, int(np.sqrt(N)) + 1):
    for j in range(2, int(np.sqrt(N)) + 1):
        if mask[j]:
            mask[j*j::j] = False

# print indices where mask is True
print(np.nonzero(mask)[0])

##>>> 
##= RESTART: D:\NEW_筆電的\網路免費軟體資料\Python教學\Python數論\用Python列出質數_加速版本\Eratosthenes厄爾多拉氏篩法用python_Burton3.2\scipy-lecture-notes-master-codes_intro_numpy_solutions\1_0_prime_sieve.py
##[   2    3    5 ... 9949 9967 9973]

2.5.2 作者改寫河西朝雄 C codes 的版本

河西朝雄的寫法是跟 SciPy 官網說明手冊的版本的最佳化的寫法一樣的設計!

Problem: 注意加速的部分是,
外部 loop 只跑到 ceil(sqrt(n))+1,
內部 loop 則自 i*i 開始跑,
(內部 loop跑的時候可以更加速為間格 i )
請同學自己證明一下為何可以如此?

Ex: 讀者可以想一下, 為何以下的寫法, 只需篩去 $2\sim \sqrt{N}$的倍數, 而無須求 $2\sim N$ 的倍數?
Ans: 假設 m=p*q, 則 p 與 q 至少一個小於等於 $\sqrt{N}$.

# primelist_Eratosthenes
# revised from codes in 河西朝雄 20160325
# 河西朝雄的寫法是跟 SciPy 官網說明手冊的版本的最佳化的寫法一樣的設計!
#(內部 loop跑的時候可以更加速為間格 i, 河西的版本沒使用, SciPy的有 )

##**Problem:** 注意加速的部分是,
##外部  loop 只跑到 ceil(sqrt(n))+1,
##內部 loop 則自 i*i 開始跑, 
## (內部 loop跑的時候可以更加速為間格  i, 河西的版本沒使用, SciPy的有 )
## 請同學自己證明一下為何可以如此?

import math
import time

def Eratosthenes(N):
    Prime=[]
    primeList=[]
    for i in range(0,N+1):
        Prime.append(1)
      
    # Prime=[ 1 for x in range(1,N+1)]
        
    Limit=math.ceil(math.sqrt(N))

    for i in range(2,Limit+1):
        if Prime[i]==1:
            for j in range(i*i,N+1):
                if j%i==0:
                    Prime[j]=0

    for i in range(2,N+1):
        if Prime[i]==1:
            primeList.append(i)

    return primeList


print(Eratosthenes(10000))

2.5.3 Python turtle 跑希臘數學家 Eratosthenes 發明的 Sieve of Eratosthenes 質數篩去法

• 對於 Python 是初學者且對海龜繪圖完全陌生的讀者可以參考本文作者的網誌系列,
  專門設計以海龜繪圖帶領讀者學會Python:
  從turtle海龜動畫學習Python-高中彈性課程1 http://t.csdn.cn/O7NVW link

• Python turtle 海龜繪圖的指令 , 尤其 海龜在畫布寫字的指令 turtle.write() 初學者可以參考官網手冊海龜繪圖指令:
  https://docs.python.org/3/library/turtle.html#turtle.write link

以下是本網誌作者自已設計的用 Python 海龜繪圖呈現 Sieve of Eratosthenes 質數篩去法的圖:

以下我們分數個階段帶領讀者了解程式碼的設計

先測試用 turtle 在畫布上寫下數字 turtle.write()

此處用到
海龜在畫布寫字的指令 turtle.write()

# By Prof P-J Lai MATH NKNU 2016/4/8
# 測試在畫布上寫下數字
# draw 1,2,3,4,,,,on the canvas 在畫面依序寫出1,2,3,4,,,,50

from turtle import *
import random
import time

T=Turtle()
T.shape("turtle")
T.color("red")
#T.color("yellow")
T.turtlesize(2)
bgcolor(0,0,0)#RGB   
T.pensize(5)
#T.pencolor("Green")

T.penup()
T.setpos(-650,0)
for i in range(1,51):
    #T.write(i,True)
    T.write(i) # 'Flase' means words will overlape the path
    T.fd(20)

T.setpos(-650,-50)
for i in range(51,80):
    T.write(i,True,font=('Arial', 15, 'italic'))
    T.fd(20)

Ex: 請同學練習, 在畫布上寫下整數 1 到 20, 1 到 10 排列在第一橫排, 11 到 20 排列在第二橫排, 大小與間隙自己調整, 使用何種字體, 與是否斜體、粗體等, 自己嘗試.

在寫下的某個數字上, 令 turtle 在此數字蓋印, 代表刪去

# By Prof P-J Lai MATH NKNU 2016/4/8
# drawInteger_2_eraseByStamp_Lai.py
# 在寫下的某個數字上, 令 turtle 在此數字蓋印, 代表刪去.

from turtle import *
import random
import time

T=Turtle()
T.shape("turtle")
T.color("red")
#T.color("yellow")
T.turtlesize(2)
#bgcolor(0,0,0)#RGB   
T.pensize(5)
#T.pencolor("Green")

T.penup()
T.setpos(-650,0)
for i in range(1,51):
    #T.write(i,True)
    T.write(i) # 'Flase' means words will overlape the path
    T.fd(20)

T.setpos(-650,-50)
for i in range(51,80):
    T.write(i,True,font=('Arial', 15, 'italic'))
    T.fd(20)

#reset 3 seconds
time.sleep(1)

# erase some integers 把一些數字擦除
T.shape('square')
T.turtlesize(2)
T.color('black')
T.setpos(-650,0)
for i in range(1,10):
    T.stamp()
    T.fd(20)

注意:
turtle.write(arg, move=False, align='left', font=('Arial', 8, 'normal'))
預設是 align=‘left’
move = False

如果 move = True
則海龜會讓出字的空間, 產生偏移,
如果要用 黑色蓋印把數字蓋住, 會出現沒法蓋對位置的情況, 因為
數字與海龜原有的位置會有偏移,
只能用大一點的印
但是數字一多,偏移量就會累積.

即使 move=False, 還是會偏移, 偏移量會小一些.

依據篩去法程式, 令 turtle 在篩去的數字蓋印, 代表刪去

以下是博主的程式碼, 並沒有系統化,
以下是N=500, 來設計跟呈現動畫面,
(注意: 當 N 改變時, 例如 N=100, 版面配置(數字的大小與位置)就會不適合. 後面會請同學練習, 當N=50 時 , 如何修改程式碼, 讓版面數字的排放是適當的?)

# By Prof P-J Lai MATH NKNU 2016/4/8
# drawInteger_Eratosthenes_Lai.py
# draw 1,2,3,4,,,,on the canvas 在畫面依序寫出1,2,3,4,,,,50
# And let turtle jump and stamp to show Eratosthenes sieve
# T.write(i,True,font=('Arial', 15, 'italic'))

#此處可以改進成一般最快的寫法:
#for j in range(i*i, m, i):
#               numbers[j] = False
#20220307

from turtle import *
from random import *
import math

T=Turtle()
T.shape("turtle")
#T.color("red")
T.color("yellow")
T.turtlesize(1)
bgcolor(0,0,0)#RGB   
T.pensize(5)
#T.pencolor("Green")
#T.speed(0)
T.speed(0)
#tracer(2,10)
T.ht()

#####################################
# print number 1 ~ 494
T.penup()
T.setpos(-650,300)
for i in range(1,39):
    T.write(i)
    T.fd(35)

T.setpos(-650,250)
for i in range(39,39+38):
    T.write(i)
    T.fd(35)

T.setpos(-650,200)
for i in range(39+38,39+2*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,150)
for i in range(39+2*38,39+3*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,100)
for i in range(39+3*38,39+4*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,50)
for i in range(39+4*38,39+5*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,0)
for i in range(39+5*38,39+6*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,-50)
for i in range(39+6*38,39+7*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,-100)
for i in range(39+7*38,39+8*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,-150)
for i in range(39+8*38,39+9*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,-200)
for i in range(39+9*38,39+10*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,-250)
for i in range(39+10*38,39+11*38):
    T.write(i)
    T.fd(35)

T.setpos(-650,-300)
for i in range(39+11*38,39+12*38):
    T.write(i)
    T.fd(35)

############################################

def coord(n):
    row=(n-1)//38
    column=(n-1)%38
    return(-650+column*35,300-row*50)
    

#Sieve of Eratosthenes
#此處可以改進成一般最快的寫法:
#for j in range(i*i, m, i):
##                numbers[j] = False
#20220307



N=494
Limit=math.ceil(math.sqrt(N))
Prime=[]
primeList=[]

for i in range(0,N+1):
    Prime.append(1)
  
# Prime=[ 1 for x in range(1,N+1)]
    
T.speed(6)
for i in range(2,Limit+1):
    if Prime[i]==1:
        T.color(random(),random(),random())
        for j in range(i*i,N+1):
            if j%i==0:
                T.setpos(coord(j))
                T.stamp()
                Prime[j]=0

for i in range(2,N+1):
    if Prime[i]==1:
        primeList.append(i)

print(primeList)

輸出的畫面:

>>> 
= RESTART: D:\NEW_筆電的\網路免費軟體資料\Python教學\Python數論\turtle跑Eratosthenes\drawInteger_Eratosthenes_Lai.py
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491]

Ex1: 當 N 改變時, 例如 N=100, 原程式碼呈現的版面配置(數字的大小與位置)就會不適合. 後面會請同學練習, 當N=50 時 , 如何修改程式碼, 讓版面數字的排放是適當的?

Ex2: 改進成, 當 N 改變時, 例如 N=100, 版面配置(數字的大小位置)會自動調整.

Ex: 此處可以改進成一般網路上效能最快的寫法(例如SciPy官網篩去法例子的寫法) (revised 20220307):
迴圈增加一個 間格 i 的跳躍!

for j in range(i*i, m, i):
    numbers[j] = False

Ex2: 改進成, 當 N 改變時, 例如 N=100, 版面配置(數字的大小位置)會自動調整.

2.6 Fermat primlity test 費馬質數檢測

理論介紹可以參考:

• 李永樂老師, 如何快速篩選質數？費馬素性檢驗和米勒-拉賓測試, https://youtu.be/E5vtn_OIe-c link
• 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 14.3.1.
• Fermat primality test, Wiki,
  https://en.wikipedia.org/wiki/Fermat_primality_test link.
• Carmichael number, Wiki, https://en.wikipedia.org/wiki/Carmichael_number link

以下常稱為費馬小定理:

Theorem 5.1 Fermat’s theorem Let $p$ be a prime and suppose that $p\nmid a$. Then $a^{p-1}\equiv 1(modp)$.

費馬質數檢測 則是利用 費馬小定理的反向敘述, 反向敘述不一定成立,

反向敘述 發生的狀況可以分兩種:

1. n 是合成數(例如 n=341=11*31), 可以找到 a 使得 $a^{341-1}\equiv 1(mod341)$, 例如 a=2, 4, 8 等等,
   此時我們稱 n=341 是關於 2 的費馬偽質數, 或 341 是關於 8 的費馬偽質數, 等等.
   (Ref: 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 人民郵電, 2020.)
2. 存在一種數, 叫 Carmichael number 卡米歇爾數, 定義 n 為一個卡米歇爾數, 如果使得 所有與 n 互質的數 a, 都有 $a^{n-1}\equiv 1(modn)$, 例如 n= 561, 1105 等等, 都符合 卡米歇爾數的定義.

以下我們用 列表解析式 list comprehension 的語法, 輕鬆檢驗 341 的狀況, 發現會讓 341 通過檢驗的 a 很多, 例如 a=2, 4, 8 等等, 即使存在很多這樣的 2, 4, 8 等等的使得 $2^{341-1}\equiv 1(mod341)$, 還是無法理論上保證, 341是質數,
如果早已知道341是合成數, 就稱此種 a 稱為 Fermat liar費馬騙子數.
如果是使得 $a^{341-1}\not\equiv 1(mod341)$, 此種 a 稱為 Fermat witness費馬證人數, 例如 a=3, 6, 9 等等, 這樣的 Fermat witness 費馬證人數, 只要找到一個, 就可以使我們立刻斷定 341 是合成數!
($\because (P\Rightarrow Q)\equiv (\sim Q\Rightarrow \sim P)$).

註: 列表解析式 list comprehension 也有 iterator 的語法, 就是把中括號改成小括號即可.
例如把
[a**340 % 341 for a in range(2,341)]
改成
(a**340 % 341 for a in range(2,341)])

列表解析式觀察費馬偽質數

列表解析式 list comprehension 的語法:
[a**340 % 341 for a in range(2,341)]
以下我們發現有的 算出來是 1, 有的不是, 但是必須是與341互質的數算出的不是1, 我們才能斷定 341是合成數,
底下列出的結果, 由於只有部分是與 341 互質, 所以難以一目瞭然,

>>> 31*11
341
>>>[ i**340 % 341 for i in range(2,341)]
[1, 56, 1, 67, 56, 56, 1, 67, 67, 253, 56, 67, 56, 1, 1, 56, 67, 56, 67, 67,
 253, 1, 56, 56, 67, 1, 56, 1, 1, 155, 1, 187, 56, 1, 67, 56, 56, 1, 67, 67,
 67, 56, 253, 56, 1, 1, 56, 67, 56, 67, 67, 67, 1, 242, 56, 67, 1, 56, 1,
 1, 155, 1, 1, 56, 187, 67, 56, 56, 1, 67, 67, 67, 56, 67, 56, 187, 1, 56,
 67, 56, 67, 67, 67, 1, 56, 56, 253, 1, 56, 1, 1, 155, 1, 1, 56, 1, 67,
 242, 56, 1, 67, 67, 67, 56, 67, 56, 1, 1, 242, 67, 56, 67, 67, 67, 1,
 56, 56, 67, 1, 242, 1, 1, 155, 1, 1, 56, 1, 67, 56, 56, 187, 67, 67,
67, 56, 67, 56, 1, 1, 56, 67, 242, 67, 67, 67, 1, 56, 56, 67, 1, 56, 1, 187, 155, 1, 1, 56, 1, 67, 56, 56, 1, 67, 253, 67, 56, 67, 56, 1, 1, 56,
67, 56, 67, 253, 67, 1, 56, 56, 67, 1, 56, 1, 1, 155, 187, 1, 56, 1,
67, 56, 56, 1, 67, 67, 67, 242, 67, 56, 1, 1, 56, 67, 56, 67, 67,
67, 187, 56, 56, 67, 1, 56, 1, 1, 155, 1, 1, 242, 1, 67, 56, 56,
1, 67, 67, 67, 56, 67, 242, 1, 1, 56, 67, 56, 67, 67, 67, 1, 56,
  242, 67, 1, 56, 1, 1, 155, 1, 1, 56, 1, 253, 56, 56, 1, 67, 67, 67, 56, 67, 56, 1, 187, 56, 67, 56, 67, 67, 67, 1, 56, 56, 67, 187, 56, 1, 
 1, 155, 1, 1, 56, 1, 67, 56, 242, 1, 67, 67, 67, 56, 67, 56, 1, 1, 56, 253, 56, 67, 67, 67, 1, 56, 56, 67, 1, 56, 187, 1, 155, 1, 1, 56, 1,
 67, 56, 56, 1, 253, 67, 67, 56, 67, 56, 1, 1, 56, 67, 56, 253,
  67, 67, 1, 56, 56, 67, 1, 56, 1, 1]

Ex: 因為只需檢驗與 n=341 互質的 a, 請使用 math.gcd(),
結合 列表解析式 list comprehension 的 if 語法:
將不互質的 a 去掉, 再看結果, 是否更清楚.
Ans:

[a**340 % 341 for a in range(2,341) if math.gcd(a,341)==1]

將不互質的 a 去掉, 則結果相當清楚, 一目了然,
仍然出現很多不是 1 的結果, 顯示 341 是合成數.

import math
[ i**340 % 341 for i in range(2,341) if math.gcd(i,341)==1]
[1, 56, 1, 67, 56, 56, 1, 67, 67, 56, 67, 56, 1, 1, 56, 67, 56, 67, 67, 1,
 56, 56, 67, 1, 56, 1, 1, 1, 56, 1, 67, 56, 56, 1, 67, 67, 67, 56, 56, 1, 1,
56, 67, 56, 67, 67, 67, 1, 56, 67, 1, 56, 1, 1, 1, 1, 56, 67, 56, 56, 
1, 67, 67, 67, 56, 67, 56, 1, 56, 67, 56, 67, 67, 67, 1, 56, 56, 1, 56, 1, 1, 1, 1, 56, 1, 67, 56, 1, 67, 67, 67, 56, 67, 56, 1, 1, 67, 56, 67, 67, 67,
 1, 56, 56, 67, 1, 1, 1, 1, 1, 56, 1, 67, 56, 56, 67, 67, 67, 56, 67, 56, 1, 1, 56, 67, 67, 67, 67, 1, 56, 56, 67, 1, 56, 1, 1, 1, 56, 1, 67, 56, 56, 1, 67, 67, 56, 67, 56, 1, 1, 56, 67, 56, 67, 67, 1, 56, 56, 67, 1, 56, 1, 1, 1, 56, 1, 67, 56, 56, 1, 67, 67, 67, 67, 56, 1, 1, 56, 67, 56, 67, 67, 67, 56, 56, 67, 1, 56, 1, 1, 1, 1, 1, 67, 56, 56, 1, 67, 67, 67, 56, 67, 1, 1, 56, 67, 56, 67, 67, 67, 1, 56, 67, 1, 56, 1, 1, 1, 1, 56, 1, 56, 56, 1, 67, 67, 67, 56, 67, 56, 1, 56, 67, 56, 67, 67, 67, 1, 56, 56, 67, 56, 1, 1, 1, 1, 56, 1, 67, 56, 1, 67, 67, 67, 56, 67, 56, 1, 1, 56, 56, 67, 67, 67, 1, 56, 56, 67, 1, 56, 1, 1, 1, 56, 1, 67, 56, 56, 1, 67, 67, 56, 67, 56, 1, 1, 56, 67, 56, 67, 67, 1, 56, 56, 67, 1, 56, 1, 1]

Ex: 用 numpy.where() 的指令, 抓出 使得 $a^{341-1}\equiv 1(mod341)$ 成立的所有 a (小於 341).
Ans:

小心 , 由於 NumPy 類似 C 或是 Matlab 的科學計算, 整數太大與會溢位, 不像原生的 Python 可以任意長整數, 所以最好先用 原生的 Python 計算大的次方, 再轉成 numpy.array 的格式

>>> list341 = [ i**340 % 341 for i in range(2,341)]
>>> import numpy
>>> listNp341 = numpy.array(list341)
>>>> numpy.where(listNp341==1)
(array([  0,   2,   6,  13,  14,  21,  25,  27,  28,  30,  33,  37,  44,
        45,  52,  56,  58,  59,  61,  62,  68,  76,  83,  87,  89,  90,
        92,  93,  95,  99, 106, 107, 114, 118, 120, 121, 123, 124, 126,
       137, 138, 145, 149, 151, 154, 155, 157, 161, 168, 169, 176, 180,
       182, 183, 186, 188, 192, 199, 200, 211, 213, 214, 216, 217, 219,
       223, 230, 231, 238, 242, 244, 245, 247, 248, 250, 254, 261, 269,
       275, 276, 278, 279, 281, 285, 292, 293, 300, 304, 307, 309, 310,
       312, 316, 323, 324, 331, 335, 337, 338], dtype=int64),)
>>>index341 = numpy.where(listNp341==1)
>>>index341
(array([  0,   2,   6,  13,  14,  21,  25,  27,  28,  30,  33,  37,  44,
        45,  52,  56,  58,  59,  61,  62,  68,  76,  83,  87,  89,  90,
        92,  93,  95,  99, 106, 107, 114, 118, 120, 121, 123, 124, 126,
       137, 138, 145, 149, 151, 154, 155, 157, 161, 168, 169, 176, 180,
       182, 183, 186, 188, 192, 199, 200, 211, 213, 214, 216, 217, 219,
       223, 230, 231, 238, 242, 244, 245, 247, 248, 250, 254, 261, 269,
       275, 276, 278, 279, 281, 285, 292, 293, 300, 304, 307, 309, 310,
       312, 316, 323, 324, 331, 335, 337, 338], dtype=int64),)
>>>index341[0]
array([  0,   2,   6,  13,  14,  21,  25,  27,  28,  30,  33,  37,  44,
        45,  52,  56,  58,  59,  61,  62,  68,  76,  83,  87,  89,  90,
        92,  93,  95,  99, 106, 107, 114, 118, 120, 121, 123, 124, 126,
       137, 138, 145, 149, 151, 154, 155, 157, 161, 168, 169, 176, 180,
       182, 183, 186, 188, 192, 199, 200, 211, 213, 214, 216, 217, 219,
       223, 230, 231, 238, 242, 244, 245, 247, 248, 250, 254, 261, 269,
       275, 276, 278, 279, 281, 285, 292, 293, 300, 304, 307, 309, 310,
       312, 316, 323, 324, 331, 335, 337, 338], dtype=int64)
>>>len(index341[0])
99
>>> index341+ numpy.ones((1, len(index341)))*2
array([[  2.,   4.,   8.,  15.,  16.,  23.,  27.,  29.,  30.,  32.,  35.,
         39.,  46.,  47.,  54.,  58.,  60.,  61.,  63.,  64.,  70.,  78.,
         85.,  89.,  91.,  92.,  94.,  95.,  97., 101., 108., 109., 116.,
        120., 122., 123., 125., 126., 128., 139., 140., 147., 151., 153.,
        156., 157., 159., 163., 170., 171., 178., 182., 184., 185., 188.,
        190., 194., 201., 202., 213., 215., 216., 218., 219., 221., 225.,
        232., 233., 240., 244., 246., 247., 249., 250., 252., 256., 263.,
        271., 277., 278., 280., 281., 283., 287., 294., 295., 302., 306.,
        309., 311., 312., 314., 318., 325., 326., 333., 337., 339., 340.]])
# 以下寫法會出錯, 因為有廣播功能, 所以結果還是對,  但是容易混淆
index341+ numpy.ones((1, len(index341)))*2
array([[  2.,   4.,   8.,  15.,  16.,  23.,  27.,  29.,  30.,  32.,  35.,
         39.,  46.,  47.,  54.,  58.,  60.,  61.,  63.,  64.,  70.,  78.,
         85.,  89.,  91.,  92.,  94.,  95.,  97., 101., 108., 109., 116.,
        120., 122., 123., 125., 126., 128., 139., 140., 147., 151., 153.,
        156., 157., 159., 163., 170., 171., 178., 182., 184., 185., 188.,
        190., 194., 201., 202., 213., 215., 216., 218., 219., 221., 225.,
        232., 233., 240., 244., 246., 247., 249., 250., 252., 256., 263.,
        271., 277., 278., 280., 281., 283., 287., 294., 295., 302., 306.,
        309., 311., 312., 314., 318., 325., 326., 333., 337., 339., 340.]])
# 以下才是較穩健的寫法:
>>>index341[0]+ numpy.ones((1, len(index341[0])))*2
array([[  2.,   4.,   8.,  15.,  16.,  23.,  27.,  29.,  30.,  32.,  35.,
         39.,  46.,  47.,  54.,  58.,  60.,  61.,  63.,  64.,  70.,  78.,
         85.,  89.,  91.,  92.,  94.,  95.,  97., 101., 108., 109., 116.,
        120., 122., 123., 125., 126., 128., 139., 140., 147., 151., 153.,
        156., 157., 159., 163., 170., 171., 178., 182., 184., 185., 188.,
        190., 194., 201., 202., 213., 215., 216., 218., 219., 221., 225.,
        232., 233., 240., 244., 246., 247., 249., 250., 252., 256., 263.,
        271., 277., 278., 280., 281., 283., 287., 294., 295., 302., 306.,
        309., 311., 312., 314., 318., 325., 326., 333., 337., 339., 340.]])

>>>numpy.ones((1, len(index341)))*2
array([[2.]])
>>>numpy.ones((1, len(index341[0])))*2
array([[2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,
        2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,
        2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,
        2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,
        2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,
        2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2.,
        2., 2., 2.]])

以上列出的 2,4,8,15,16,就是使得 $a^{341-1}\equiv 1(mod341)$ 的 a.
以上是屬於第1種的情形, 就是合成數 341, 會同時有 費馬騙子數, 與費馬證人數.

第2種情形, 稱為 Carmichael number 卡米歇爾數, 則是所有與 n 互質的數 a, 都是費馬證人數!
此種 Carmichael number 卡米歇爾數, 儘管數量很稀少, 還是使得費馬質數檢測的使用不如米勒拉賓質數檢測等其他機率型的質數檢測演算法,
但是 OpenPFGW, PGP 是使用費馬質數檢測.

最小的 卡米歇爾數 為 $561=3\times 11\times 7$,

接下來的 6 個卡米歇爾數 (sequence A002997 in the OEIS):

1105 = 5×13×17 (4 | 1104, 12 | 1104, 16 | 1104)
1729 = 7×13×19 (6 | 1728, 12 | 1728, 18 | 1728)
2465 = 5×17×29 (4 | 2464, 16 | 2464, 28 | 2464)
2821 = 7×13×31 (6 | 2820, 12 | 2820, 30 | 2820)
6601 = 7×23×41 (6 | 6600, 22 | 6600, 40 | 6600)
8911 = 7×19×67 (6 | 8910, 18 | 8910, 66 | 8910)

>>> [a**560 % 561  for a in range(2,561)]
[1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 34, 
375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528,
 34, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 
 1, 408, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 
 528, 1, 34, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 34, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 408, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 34, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 34, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 408, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 34, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 187, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 408, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 34, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 34, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 408, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 34, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 34, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 408, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 34, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 34, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 408, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 187, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 34, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 408, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 34, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 34, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 408, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 34, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 34, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 408, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 34, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 34, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1]

Ex: 檢驗以上 Wiki 列出的數 561, 1105 等符合 卡米歇爾數的定義. (參考: 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 14.3.1, P.231, 指出, 比 $2.5\times 10^{10}$ 小的 卡米歇爾數只有 2163 個.)

以下檢驗 561 的計算狀況
列表解析式 list comprehension 的語法:
[a**560 % 561 for a in range(2,561)]
底下列出的結果, 由於只有部分是與 561 互質, 所以難以一目瞭然,

>>> [a**560 % 561  for a in range(2,561)]
[1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 34,
 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528,
  34, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1,
   1, 408, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1,
    528, 1, 34, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 34, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 408, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 34, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 34, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 408, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 34, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 187, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 408, 1, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 34, 375, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 34, 1, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 408, 1, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 34, 375, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 34, 1, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 408, 1, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 34, 375, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 34, 154, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 408, 1, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 187, 375, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 34, 1, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 408, 1, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 34, 375, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 34, 1, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 408, 1, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 34, 375, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 34, 1, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 408, 1, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 34, 528, 1, 1, 375, 1, 1, 375, 1, 1, 375, 1, 154, 375, 1, 1, 375, 34, 1, 375, 1, 1, 375, 154, 1, 375, 1, 1, 375, 1, 1, 375, 1, 1]

因為只需檢驗與 n=561 互質的 a, 我們使用 math.gcd(),
結合 列表解析式 list comprehension 的 if 語法:

[a**560 % 561 for a in range(2,561) if math.gcd(a,561)==1]

將不互質的 a 去掉, 則結果相當清楚, 一目了然!

>>> import math
>>> [a**560 % 561  for a in range(2,561) if math.gcd(a,561)==1]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

發現, 所有與 n=561 互質的數 a, 確實都會 $a^{n-1}\equiv 1(modn)$! 所以 n=561 確實是 米歇爾數.

Ex: 用 Python 找出 不在以上 Wiki 列出的 卡米歇爾數.

費馬質數檢測的 flowchart

利用 費馬小定理的反向敘述來檢測質數, 就是, 假設你想檢測 n 是否為質數,
則隨機挑幾個與 n 互質的數 a, 看 $a^{n-1}\equiv 1(\mathrm{mod}n)$ 是否成立,
如果有遇到一個 a 使得 $a^{n-1}\not\equiv 1(\mathrm{mod}n)$, 則可以立刻得到 $n$ 不是質數,
如果嘗試夠多個 a, 都得到 $a^{n-1}\equiv 1(\mathrm{mod}n)$, 此時只能說 n 是質數的可能性極大.

例如 341, 有 費馬騙子數 2,4,8,15,16, 使得
$a^{341-1}\equiv 1(\mathrm{mod}341)$ 成立.
另外也有費馬證人數 Fermat witness, 例如 a=3, 6, 9 等等, 使得 $a^{341-1}\not\equiv 1(\mathrm{mod}341)$ 成立,
這樣的數, 可以使我們立刻斷定 341 是合成數, 此種數, 數量也很多.

根據 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 14.3.1, P.231, 指出,
隨機選出 a, “n 不是以 a 為底的偽質數” 的機率至少是 $50\%$,
故 , 如果只做一次檢驗, 成功的機率至少是 $50\%$,
做兩次檢驗, 成功的機率至少是 $75\%$,
如果做10次檢驗, 成功的機率至少是 $99.9\%$.

# Fermat primality test
# By Prof. P-J Lai MATH NKNU 20220302
##>>> help('random.randint')
##Help on method randint in random:
##
##random.randint = randint(a, b) method of random.Random instance
##    Return random integer in range [a, b], including both end points.

# a**(n-1)%n 
# 等價於 pow(a, n-1, n)

##The algorithm can be written as follows:
##
##Inputs: n: a value to test for primality, n>3; k: a parameter that determines the number of times to test for primality
##Output: composite if n is composite, otherwise probably prime
##Repeat k times:
##Pick a randomly in the range [2, n − 2]
##If {\displaystyle a^{n-1}\not \equiv 1{\pmod {n}}}a^{{n-1}}\not \equiv 1{\pmod  n}, then return composite
##If composite is never returned: return probably prime
##The a values 1 and n-1 are not used as the equality holds for all n and all odd n respectively, hence testing them adds no value.


import random

def FermatPrimalityTest(n,k):

    times = k
    while times >= 1:
        a = random.randint(2,n-2)
        if a**(n-1)%n != 1:
        # 等價於 pow(a, n-1, n)
            return False

        times = times - 1

    return True

##輸出
##>>> 
##============= RESTART: C:/Users/user/Desktop/FermatPrimalityTest.py ============

##>>> FermatPrimalityTest(101,3)
##True
##>>> FermatPrimalityTest(195,3)
##False
##>>> FermatPrimalityTest(421,3)
##True
##>>> FermatPrimalityTest(91,3)
##False
##>>> FermatPrimalityTest(91,3)
##False
##>>> FermatPrimalityTest(91,3)
##False          

Ex: 寫一個 Python 程式搜尋例如 十萬100000或百萬1000000以內的所有卡米歇爾數.

2.7 Miller-Rabin test 米勒拉賓質數檢測

理論介紹可以參考:

• 李永樂老師, 如何快速篩選質數？費馬素性檢驗和米勒-拉賓測試, https://youtu.be/E5vtn_OIe-c link

• wiki 米拉賓檢測, link

以下主要參考 Wiki 的說明:

Lemma1: If p is prime, then the root of $X^2\equiv 1\mathrm{mod}n$ should be either $X\equiv 1\mathrm{mod}n$ or $X\equiv -1\mathrm{mod}n$.
Pf: 請參考 Wiki.

由 Lemma1 可以得到

Lemma2: If n is prime, then for any a ∈ { 2 , 3 , 4 , , , , n − 1 } a \in \{2,3,4,,,,n-1\} a∈{2,3,4,,,,n−1}, we have that either
$a^d\equiv 1\mathrm{mod}n$ for some odd number $d>0.-----(1)$
or
$a^{2^{r}d}\equiv -1\mathrm{mod}n$ for some $r\ge 0.---------(2)$
Pf:
假設 n 為質數
取 a ∈ { 2 , 3 , 4 , , , , n − 1 } a \in \{2,3,4,,,,n-1\} a∈{2,3,4,,,,n−1}
此時由 費馬小定理, 有
$a^{n-1}\equiv 1\mathrm{mod}n$
n 為質數,
如果 n=2, 則 (1) 一定成立,
如果 n>2, 則為奇數,
則有 n-1 為偶數, 可以表為 $n-1=2^{s}d$
則
$a^{n-1}\equiv 1\mathrm{mod}n$
就可以寫成
$a^{2^{s}d}\equiv 1\mathrm{mod}n$
用冪次律整理, 抽出 2 次方:
$(a^d{)}^{2^s}\equiv 1\mathrm{mod}n$,
$((a^d{)}^{2^{s-1}2}\equiv 1\mathrm{mod}n$,
$((a^d{)}^{2^{s-1}}{)}^2\equiv 1\mathrm{mod}n$.
根據前面的小引理Lemma1, 得到平方根只可能是 $\pm 1$,
也就是
$(a^d{)}^{2^{s-1}}\equiv 1\mathrm{mod}n$
or
$(a^d{)}^{2^{s-1}}\equiv -1\mathrm{mod}n$

或者 更清楚一點表達
Let $X:=(a^d{)}^{2^{s-1}}$
則有
$X^2\equiv 1\mathrm{mod}n$
根據前面的 Lemma1, 得到
$X\equiv 1\mathrm{mod}n$
or
$X\equiv -1\mathrm{mod}n$,
也就是 平方根 X 只可能是 $\pm 1$,

case1 : 假設我們得到的是 $X\equiv 1\mathrm{mod}n$
如果 $s>1$, 則 $X:=(a^d{)}^{2^{s-1}}$ 就可以再用冪次律整理, 抽出 2 次方:
$X=((a^d{)}^{2^{s-2}}{)}^2\equiv 1\mathrm{mod}n$

可以令 Let $X1=((a^d{)}^{2^{s-2}})$, 則有
$(X1{)}^2\equiv 1\mathrm{mod}n$
則可以再用一次 Lemma1,
得到平方根 X1 只可能是 $\pm 1$,

所以, 如果每次都 mod 等於 1, 則可以把 $2^s$ 抽光,
最後剩下
$a^d\equiv 1\mathrm{mod}n$
則得到 (1) 的狀況.

case2 : 假設我們得到的是 $X\equiv -1\mathrm{mod}n$
另一方面, 不管 $s=1$ 或 $s>1$,
如果是mod 等於 -1 的情況: $X\equiv -1\mathrm{mod}n$.
則得到 (2) 的狀況: $a^{2^{r}d}\equiv -1\mathrm{mod}n$ for some $r\ge 0$.

經過上面討論, 就知道一定會落入 (1) 或 (2) 的狀況.
得證.$■$

根據 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 14.3.2, P.231, 指出, 米勒拉賓質數檢測是在 費馬質數檢測的基礎上又添加了另外的判斷條件, 令 卡米歇爾數更難通過檢測!

程式碼可以看:

• 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 14.3.3, P.233.
• 雲帝, Python Miller Rabin 米勒-拉宾素性检验,
  link

以下參考 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 14.3.3, P.233 的程式碼:

 # 以下參考 王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 14.3.3, P.233 的程式碼
 # By Prof. P-J Lai MATH NKNU 20220403
 
from random import randrange

# 以 a 為底數, 檢查 n 是否為合數, n-1 = 2^k*d 

def checkComposite(a,k,d,n):
    x = pow(a, d, n)
    if x == 1 or x == n-1:
        return False

    for _ in range(k-1):
        x = pow(x, 2, n)
        if x == n-1:
            return False

    return True

# n 是被檢驗的數, m 是精準度

# 判斷 n 是否不是合數

def miller_rabin(n, m):
    if n == 2 or n == 3:
        return True
    if n % 2 == 0:
        return False
    k = 0
    d = n-1
    # 計算 k  與 d
    while d % 2 ==0:  
        k += 1
        d = d//2
        
    for _ in range(m):
        a = randrange(2, n-1)
        if checkComposite(a,k,d,n):
            return False
        return True
  

3. 用 sympy 求因數與質數

sympy 求因數的指令

sympy.factorint 列出因數與因數之次方(重數)
sympy.primefactors 列出質因數

>>> import sympy
>>> sympy.primefactors(120)
[2, 3, 5]
>>> sympy.primefactors(17894)
[2, 23, 389]
>>> sympy.factorint(120)
{2: 3, 3: 1, 5: 1}
>>> sympy.factorint(17894)
{2: 1, 23: 1, 389: 1}

sympy才有與質數與密碼有關的指令

prevprime, nextprime, isprime, randprime, prime(), sieve(), primerange(a, b)
primerange(a, b) 範圍包含 a, 不包含 b.
prime(k) 產生第k個質數

>>> import sympy
>>> print([i for i in sympy.primerange(1,30)])
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
>>> print([i for i in sympy.sieve.primerange(1,30)])
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

P.310 The Sieve method, primerange, is generally faster but it will
occupy more memory as the sieve stores values.

>>> sympy.isprime(100003)
True
>>> sympy.randprime(1,1000)
331

>>> sympy.factorint(100001)
{11: 1, 9091: 1}

4. primePy PyPI: Python 的其他質數模組

5. 測量以上各種不同求質數或質數檢測之寫法, 使用 time.time(), 比較它們的效能

Ex: 測量以上各種不同求質數或質數檢測之寫法, 使用 time.time(), 及 timeit.timeit 比較它們的效能

6. 展現 Python 呈現超長整數的能力: 孿生質數猜想: 用Python呈現 2009 為止最大的孿生質數 (Using Python to show the biggest twin prime pairs up to 2009)

Ex: 寫出 Python 程式, 列出給定範圍之所有孿生質數.

展現 Python 呈現超長整數的能力
Pyhton的整數位數在內存記憶允許之下, 真的整數幾乎無限制,

整數太長時, 在IDLE無法顯示,會卡住, 他會顯示建議你複製或用view, 用view還是當住,
我是按右鍵copy, 再貼入純文字檔, 就可以了

應數一整數論: Burton 課本 sec3.3 The Goldbach conjecture P. 50
20200515
用Python算2009發現最大之孿生質數,
孿生質數猜想
孿生質數就是相隔2的兩個質數,
2009年發現最大的是
65516468355*pow(2,333333)-1
65516468355*pow(2,333333)+1
是否存在無限多組, 目前還未證明出, 但是電腦搜尋似乎是正確的猜想
證明Pyhton的整數位數在內存記憶允許之下, 真的幾乎無限制,
在IDLE無法顯示,會卡住, 他會顯示建議你複製或用view, 用view還是當住,
我是按右鍵copy, 再貼入純文字檔, 就可以了
C++ 就要引入大(任意長)整數的class,或是單純C就自己寫一個, 河西朝雄那本有教, 基本上是用拼接的方式,
2009最大孿生質數 65516468355*pow(2,333333)-1
Python執行
下方數字結果為貼入純文字檔之後段數字( 太長了! 前段沒貼出)

>>> 65516468355*2**333333–1

,
0486204337045830221697947493623893208668843236296432330579337007454613523827007270017022742840808830512021099751772089661025726082430148634950529172853905398252455950441387165528658094024941734837889747048797176590531811034084024893902785968263923457852989284900876894451850607250159224932244297921003913913871058528427574279544850405389475352929766783627081676782738253911305468497669721299043953090587286157634369169664899480688857366800375578871448984709935499333259890644537584669287564834661609820559944905794803947018011664859306262616169855607254851038466500959334652515316978420118362309783716550790213164548090179110078222154290489633403899510017232142570096708316232747641485642658928268494424639302394298903475238312901914372524456364615166666035851403790970205218045734933905677045844011921952784183957697550021803201409719944091070402482721436819252332247832126191716837863031239249901676116836904951120854350189972990029539733826953721723922362125737807974748278081705412230891338699189206372317216443616890716159

>>> 65516468355*2**333333+1

,
0486204337045830221697947493623893208668843236296432330579337007454613523827007270017022742840808830512021099751772089661025726082430148634950529172853905398252455950441387165528658094024941734837889747048797176590531811034084024893902785968263923457852989284900876894451850607250159224932244297921003913913871058528427574279544850405389475352929766783627081676782738253911305468497669721299043953090587286157634369169664899480688857366800375578871448984709935499333259890644537584669287564834661609820559944905794803947018011664859306262616169855607254851038466500959334652515316978420118362309783716550790213164548090179110078222154290489633403899510017232142570096708316232747641485642658928268494424639302394298903475238312901914372524456364615166666035851403790970205218045734933905677045844011921952784183957697550021803201409719944091070402482721436819252332247832126191716837863031239249901676116836904951120854350189972990029539733826953721723922362125737807974748278081705412230891338699189206372317216443616890716161

Appendix: 比較程式執行時間之語法

0.1 練習用 time 模組量時間程式執行時間

例如:

import time
start = time.time()
Factor(10000000)
print(time.time() - start)


start = time.time()
Factor_half(10000000)
print(time.time() - start)

0.2 練習用 timeit 模組 timeit.timeit()量時間

另一種測量時間之語法:

針對時間很短暫的程式之測量用
預設是執行1000000，一百萬次

引數說明：
import timeit
timeit.timeit(stmt, setup,timer, number)
stmt : statement的縮寫，你要測試的程式碼或者語句，純文字，預設值是 “pass”
setup : 在執行 stmt 前的配置語句，純文字，預設值也是 “pass”
timer : 計時器，一般忽略這個引數
number : stmt 執行的次數，預設是1000000，一百萬

from timeit import timeit
print(timeit('Factor(100)',  setup='from __main__ import Factor'))
print(timeit('Factor_half(100)', setup='from __main__ import Factor_half'))

Reference

• 快速篩選質數
  理論介紹可以參考: 李永樂老師, 如何快速篩選質數？費馬素性檢驗和米勒-拉賓測試, https://youtu.be/E5vtn_OIe-c link

• Fermat primality test, Wiki,
  https://en.wikipedia.org/wiki/Fermat_primality_test link.

• Carmichael number, Wiki, https://en.wikipedia.org/wiki/Carmichael_number link

• wiki 米拉賓檢測, https://zh.wikipedia.org/wiki/%E7%B1%B3%E5%8B%92-%E6%8B%89%E5%AE%BE%E6%A3%80%E9%AA%8C link

• Burton, Elementary Number Theory, McGraw Hill 2011.

• AI Sweigart, Cracking Codes With Python

• 博主推薦: 以下王碩等作者之書, 重點講述, 讓讀者立刻抓住核心概念, 不多廢話, 是難得好書!
  王碩, 董文馨, 張舒行, 張ㄗㄜˊ, 李秉倫, Python 算法設計與分析, 人民郵電, 2020.

• Python turtle 海龜繪圖 初學者可以參考
  官網手冊海龜繪圖指令:
  https://docs.python.org/3/library/turtle.html#turtle.write link

• Python turtle 海龜繪圖 初學者可以參考本文作者的網誌:
  從turtle海龜動畫學習Python-高中彈性課程1 http://t.csdn.cn/O7NVW link