카테고리 없음
(정수론/파이썬) 수론 관련 함수 모음 (number_util.py)
미친토끼
2021. 9. 16. 15:07
정수론 함수 목록 2021-9-25 기준
isPrime(num): 소수인지 판정
isPrime2(num): 큰수에 대해 페르마 정리와 카마이클수를 이용해 소수 판정
isPrimeWilson(n): 차례곱을 이용해 소수 판정. 숫자가 조금만 커져도 시간이 많이 걸림.
isPrimeGuess(n): 큰수에 대해 페르마 정리와 카마이클수를 이용해 소수 판정
gcd(a, b): 유클리드 호제법으로 최대공약수 구하기
gcd2(a, b): 최대공약수 정의대로 삽질로 구하는 함수
lcm(a, b): 최소공배수 구하기. 각기 배수들의 교집합에서 제일 작은 수
lcd(a, b): 최소공약수. 교집합 방식으로
lcm2(a, b): 최소공배수. 중학교 1학년 때 배운 식대로
getSquarePrime(): N^2 + 1 형태 소수 구하기
findTwinPrimes(): 쌍둥이 소수 구하기
factorize(num): 소인수분해. 구버전
factorize2(value): 소인수분해. 개선 버전
factorize_list(num): 소인수분해 결과를 목록으로. 지수는 빼고 [2, 3, 7]
goldbach(num): 골드바흐의 추정 검증하기. 4 이상의 모든 짝수는 소수 2개의 합으로 표현된다.
linear_equation_compute(a, b, c): 선형방정식 해 구하기. 내 방식대로.
linear_equation(a, b): 선형방정식 해 구하기, 책에 나온 대로
congruence_equation(a, c, m): 합동방정식 해 구하기
congruence_equation_poly(poly, m): 다항식 합동식의 해를 구한다
dirichlet(a, m, max_num): 디리클레의 등차수열 정리를 확인
euler_phi(num): 서로소인 숫자를 찾아 그 갯수를 반환하는 함수
euler_phi_list(num): 서로소인 숫자를 찾아 그 목록을 반환하는 함수
euler_phi2(m): 공식을 이용해 오일러 파이값을 구하는 함수
euler_phi3(num): 인수분해하여 오일러 파이값을 구하는 함수
crt(b, m, c, n): 중국인의 나머지 정리
crt2(*args): 중국인의 나머지 정리. 가변적인 개수의 인자를 받는다.
getDivisor(num): 일일이 나눠보는 방식으로 약수 구하는 함수, 속도가 느림
getDivisor2(num): 소인수분해를 해서 약수 구하는 함수: 향상된 버전
divisorsCount(num): 약수의 개수만 구하는 함수. 소인수들의 지수 이용
sigma_fun(num): 시그마 함수: 자신을 포함한 모든 약수의 총합을 돌려주는 함수:
예) 12의 시그마 함수값: 1, 2, 3, 4, 6, 12 을 모두 더한 값
isPerfectNumber(num): 자신을 제외한 약수들을 모두 더해서(진약수 합) 돌려준다
get_mersennePrime(): 메르센 소수 구하기, 2^N - 1 = 소수..인 경우 구하기
euclidPerfectNumber(): 유클리드 완전수 공식을 이용해 완전수 찾기
isProductPerfect(num): 곱셈완전수인지 판정하는 함수
isAmicablePair(m, n): 두 수가 친화수인지 확인
isSociableNumbers(numbers, info=False): 숫자 리스트가 사교수인지 확인
getSociableNumbers(start, end, max_try=30): 사교수의 범위를 정해주면 사교수를 찾는다
successiveSquaring(a, k, m): 연속제곱법. 진행 과정을 그대로 구현. 메시지 다수 출력
successiveSquaring2(a, k, m): 연속제곱법. 향상된 버전. 책에서 나온 것.
my_successiveSquaring(a, k, m): 연속제곱법 : 쑤튜브에서 힌트 얻어서 작성
high_congruence_equation(k, b, m): 고차합동방정식 푼다. phi(m)구하기, 일차 합동식 풀기, 연속제곱법 사용해서 합동식 풀기 등 3가지 기법 사용
cipher_book: RSA 암호에 사용하는 암호 딕셔너리
RSA_incode(msg, m, k): RSA 암호화 알고리즘
RSA_decode(ciphered_list, m, k, p, q): RSA 복호화 알고리즘
isCarmichael(m): 카마이클 수 판정. 확인 위해 일일이 나눠보는 함수
isCarmichael2(n): 카마이클 수 판정. 확인 위해 소인수들만 나눠보는 함수
compositeWitness(m): 카마이클 수 관련, m이 합성수라는 증인 목록을 만들어준다
# '친절한 수론 길라잡이' 공부를 하면서 만든 함수들
# 파일명: number_util.py
# 이 파일을 같은 폴더에 넣고, 호출할 때 from number_util import * 로 불러옴
import math
import random
# 소수 판정 함수. 소수면 True, 아니면 False 반환
def isPrime(num):
if num <= 1: return False
if num == 2 or num == 3: return True
if num % 2 == 0: return False # 짝수라면 합성수
root_val = round(num ** 0.5) + 1 # 제곱근 값을 구한다
# 이미 2를 제외했으므로 3부터 2씩 증가시키며 홀수로만 나눈다.
for i in range(3, root_val, 2): # num이 5인 경우부터 시작한다
if num % i == 0: # 중간에 나눠지면 소수가 아니다
return False
return True
# 연습문제 16.4 : 큰수에 대해 페르마 정리와 카마이클수를 이용해 소수임을 판정
# True: 소수일 확률이 높음, False: 합성수
def isPrime2(n):
MAX_NUM = 10
#n = 243034209889348945908458903908324908449804908434389747894389438943
# 2<= a < n, 10개의 수를 고른다.
a_s = []
for i in range(MAX_NUM):
a = random.randrange(2, n)
#if gcd(a, n) != 1: continue # 큰수에 대해서는 gcd를 구하는 행위 자체가 엄청난 부담
a_s.append(a) # 2 ~ n - 1
#print(a_s)
maybe_prime = True
for i in range(MAX_NUM): # 숫자 10개를 채택해서 페르마 정리를 이용해 n을 검사함
if successiveSquaring2(a_s[i], n-1, n) != 1:
maybe_prime = False # composite number
break
#if maybe_prime:
# print(f"{n} 은 아마도 소수")
#else:
# print(f"{n} 은 합성수")
return maybe_prime
# Wilson's Theorem
# p is prime if (p-1)! = -1 (mod p)
def isPrimeWilson(n):
m = n
gop = 1
while m >= 3:
gop = (gop * (m - 1)) % n
m -= 1
return gop == n - 1
# 큰수가 소수인지 아닌지 판정한다. 페르마 소정리를 이용.
# '아마도 소수'에는 아주 가끔 합성수가 끼어 있는 경우가 있다.
def isPrimeGuess(n):
MAX_NUM = 1000 # 난수를 이 수치만큼 발생
# 2와 N - 1 사이에서 10개의 수를 고른다.
# n과 서로소인 수를 고르는 게 좋겠지만 그 자체에 자원이 꽤 소모되므로 단순한 방식으로.
a_s = []
for i in range(MAX_NUM):
a = random.randrange(2, n)
if gcd(a, n) != 1: continue
a_s.append(a) # 2 ~ n - 1
#print(a_s) # 생성된 10개의 난수
maybe_prime = True
for i in range(MAX_NUM):
if successiveSquaring2(a_s[i], n-1, n) != 1: # ai(n-1) mod n 계산
maybe_prime = False # composite number
break
if maybe_prime:
print(f"{n} 은 아마도 소수")
else:
print(f"{n} 은 합성수")
# 두 수 a와 b (단 a > b) 의 최대공약수 구하기: 유클리드 호제법
def gcd(a, b):
if a == b:
return a
if a < b: # a > b 조건을 맞추기 위해 두 값을 맞바꾸기
tmp = a
a = b
b = tmp
r0 = a
r1 = b
for i in range(b): # 두 수 증 작은 수만큼 최대 반복
r2 = r0 % r1 # 몫은 필요 없으므로 나머지만 구함
if r2 == 0: # 나머지가 0이라면
return r1 # 직전의 나머지가 최대공약수임
# 값들을 한칸씩 shift 시켜서 다시 반복
r0 = r1
r1 = r2
# 최대공약수를 삽질로 구하는 함수
def gcd2(a, b):
a_set = getDivisor(a) # a의 약수를 모두 구한다
b_set = getDivisor(b) # b의 약수를 모두 구한다
#print(a_set)
#print(b_set)
c_set = set(a_set) & set(b_set) # 두 약수들의 교집합을 구한다. 공약수들이 된다.
return max(list(c_set)) # 리스트로 바꾸어 최대값을 구한다.
# 최소공배수 구하기
# lcm(a, b)라고 할 경우, 최소공배수는 a*b와 같거나 작다. a, b 가 서로소이면 최소공배수는
# a*b 이고, 공약수가 있을 경우, lcm(a, b) < a*b 이다.
# a의 배수를 a*b 이하에서 찾아서 나열하고, b의 약수를 a*b 이하에서 찾아서 나열해서,
# 교집합을 구하면 그것이 a*b 이하 공배수들인데, 그중에서 가장 작은 수가 최소공배수이다.
def lcm(a, b):
ab = a * b
a_multi = []
for i in range(1, b+1):
a_multi.append(a * i)
b_multi = []
for i in range(1, a+1):
b_multi.append(b * i)
common = set(a_multi) & set(b_multi) # 배수들의 교집합 즉, 공배수 찾기
return min(list(common)) # 리스트형으로 바꿔 최소값을 리턴한다.
# 최소공약수를 구하는 함수 (1 제외)
def lcd(a, b):
a_ = getDivisor(a) # a의 약수를 모두 구한다
b_ = getDivisor(b) # b의 약수를 모두 구한다
#print(a_)
#print(b_)
c_ = set(a_) & set(b_) # 두 약수들의 교집합을 구한다. 공약수들이 된다.
if len(c_) == 1: # 공약수가 1뿐이라면, 1을 반환: 최소공약수
return min(c_)
elif len(c_) > 1: # 공약수가 1외에도 있다면 1을 제거하고 가장 작은 수를 반환한다: 최소공약수
return min(list(c_ - {1}))
# 최소공배수 구하기
# a, b 의 최소공배수로 나눠서 두 수가 서로소일 때까지 나누고,
# 그 제수들(나누는 수)와 최종 몫들을 모두 곱하면 최소공배수가 된다. 중학교 1학년때 배우는...
def lcm2(a, b):
if gcd(a, b) == 1: # a와 b가 서로소라면, a*b를 최소공배수로 반환
return a * b
# 이제 a, b의 공약수가 있다는 말이므로, 최소공약수를 구해서 a, b를 나눈다.
divisors = []
while True:
lcd_val = lcd(a, b) # 최소공약수를 구한다.
if lcd_val != 1: # 최소공약수가 2 이상이라면
a = a / lcd_val # 그 수로 나눠서 몫을 새로 a라고 놓는다
b = b / lcd_val # 그 수로 나눠서 몫을 새로 b라고 놓는다
divisors.append(lcd_val) # 나누는 수(최소공약수)를 배열에 저장한다.
else: # 최소공약수가 1이라면 더 이상 나눌 수 없으므로
divisors.append(a) # 최후의 몫들을 배열에 저장한다.
divisors.append(b)
total = 1
for n in divisors: # 배열의 모든 원소들을 곱하여 리턴한다
total *= n
return int(total)
#print(lcm(8, 12))
#print(lcm2(32383800, 8055))
# N^2 + 1 이 소수인 경우가 무한할까? (N = 1, 2, 3...)
# N^2 + 1 형태의 소수. N = 1, 2, 3...을 대입하여 얻은 N2 + 1 형태의 수를 나열해보면,
# 몇몇 소수가 발견된다. N이 홀수이면 N^2 + 1은 짝수이므로 N=1일 때를 제외하고는 N^2 + 1이
# 소수일 리 없다. N에 짝수를 대입하는 경우만 생각해보자.
def getSquarePrime():
i = 2 # 2부터 2씩 증가시켜 대입한다
count = 0 # 소수 등장 횟수
while True:
num = i ** 2 + 1
if isPrime(num):
count += 1
print(f"{i} ^2 + 1 = {num} ({count})")
i += 2 # 홀수에는 관심없고 짝수만 취급한다.
#getSquarePrime()
# 쌍둥이 소수: 인접한 두 수(차이가 2)가 소수인 경우를 찾는다
def findTwinPrimes():
i = 5 # 현재 수
prev_prime = 3 # 이전에 등장한 소수
count = 0 # 쌍둥이 소수 등장 횟수
while True:
if isPrime(i): # 현재 수가 소수라면
if i - prev_prime == 2: # 이전의 소수(홀수)가 인접한 소수(홀수)라면
count += 1
print(f"({prev_prime}, {i}) {count}")
prev_prime = i
i += 2 # 짝수에 대해서는 굳이 소수 여부 조사를 할 필요가 없음. 홀수에 대해서만 조사함
#findTwinPrimes()
# 소인수분해 해주는 함수. 구형 버전. 속도가 다소 느림.
# 가령 인수분해 결과가 2^3 * 3^2 * 7 이라면
# {2:3, 3:2, 7:1} 반환
def factorize(num):
divisor_primes = [] # 분해된 인수를 차곡차곡 저장하는 리스트
factorize_dict = {} # 리스트의 원소 갯수를 세어 넣는 용도
# num이 소수라면 (num, 1)을 딕셔너리에 넣어 반환한다.
# 소수 판정은 숫자의 제곱근까지만 조사하면 되므로 간단하다.
if isPrime(num):
factorize_dict[num] = 1
return factorize_dict
# 인수분해 루틴은 소수 판정 루틴과는 조금 다르게 접근할 필요가 있다.
i = 2
while True:
# i가 소수인지 검사해서 소수라면 num을 i로 나눈다
# i로 나눠서 나머지가 0이라면, 인수분해가 한번 된 것이다.
if isPrime(i) and num % i == 0:
divisor_primes.append(i) # 나누는 수를 저장
num = num // i # 몫, 가령 24//2 = 12
# 12는 소수가 아니므로 반복문 계속. 소수가 나오면 그 소수를 저장하고 반복문 종료.
if isPrime(num):
divisor_primes.append(num)
break
i = 2 # 인수분해를 한번 했으므로 숫자를 다시 2부터 나눠보기로 설정
else:
i += 1 # 인수분해가 되지 않으면 젯수(나누는 수)를 1증가
for n in divisor_primes:
if n in factorize_dict:
factorize_dict[n] += 1
else:
factorize_dict[n] = 1
return factorize_dict
# 소인수분해 해주는 함수, 개선 버전
def factorize2(value):
num = value # 원래값을 보존한다.
divisor_primes = [] # 분해된 인수를 차곡차곡 저장하는 리스트
factorize_dict = {} # 리스트의 원소 갯수를 세어 넣는 용도
root_val = int(num ** 0.5)+1
i = 2
found = False
while True:
# 소수 검증 때처럼 num까지 모두 나누어볼 필요는 없고 제곱근값 까지만 나누어보면 된다.
if i >= root_val: break
# i가 소수인지 검사해서 소수라면 num을 i로 나눈다
# i로 나눠서 나머지가 0이라면, 인수분해가 한번 된 것이다.
if isPrime(i) and num % i == 0:
divisor_primes.append(i) # 나누는 수를 저장
num = num // i # 몫, 가령 24//2 = 12
root_val = int(num ** 0.5)+1
if found == False:
found = True # 찾았다는 표시를 함
# 12는 소수가 아니므로 반복문 계속. 소수가 나오면 그 소수를 저장하고 반복문 종료.
if isPrime(num):
divisor_primes.append(num)
break
else:
i += 1 # 인수분해가 되지 않으면 젯수(나누는 수)를 1증가
# 소인수분해 결과를 딕셔너리 형태로 바꿈
for n in divisor_primes:
if n in factorize_dict:
factorize_dict[n] += 1
else:
factorize_dict[n] = 1
if found == False: # 소인수를 못 찾았다면 소수이므로
factorize_dict[value] = 1
return factorize_dict
# 소인수분해하여 그 목록을 돌려준다
# 가령 소인수분해 결과가 2^3*3^2*7 이라면,
# [2, 3, 7]을 돌려준다
def factorize_list(num):
divisor_primes = [] # 분해된 인수를 차곡차곡 저장하는 리스트
if isPrime(num): # 소수 판정
return divisor_primes
# 인수분해 루틴은 소수 판정 루틴과는 조금 다르게 접근할 필요가 있다.
i = 2
while True:
# i가 소수인지 검사해서 소수라면 num을 i로 나눈다
# i로 나눠서 나머지가 0이라면, 인수분해가 한번 된 것이다.
if isPrime(i) and num % i == 0:
divisor_primes.append(i) # 나누는 수를 저장 (가령, num이 12라면 i = 2)
num = num // i # 몫, 가령 24//2 = 12
#print(f"{i}, {num}")
# 12는 소수가 아니므로 반복문 계속. 소수가 나오면 그 소수를 저장하고 반복문 종료.
if isPrime(num):
divisor_primes.append(num)
break # return divisor_primes
i = 2 # 인수분해를 한번 했으므로 숫자를 2로 다시 나눠보기로 설정
else:
i += 1 # 인수분해가 되지 않으면 젯수(나누는 수)를 1증가
return sorted(set(divisor_primes)) # 중복 제거, 정렬하여 돌려주기
# 골드바흐의 추정 : Goldbach's Conjecture
# 4 이상의 모든 짝수는 소수 두 개의 합으로 표현된다.
# ' 친절한 수론 길라잡이' 95쪽
# 4 = 2 + 2, 6 = 3 + 3, 8 = 3 + 5, 10 = 3 + 7, 12 = 5 + 7,
# 14 = 3 + 11, 16 = 3 + 13, 18 = 5 + 13, 20 = 7 + 13 ...
# 숫자 num을 인자로 받는다.
# 2(=i)부터 시작해서 num/2 까지 소수에 해당되는 걸로만 해서 num - i를 시도한다.
# num - i가 소수라면 num이 i와 num-i 라는 두 개의 소수로 표현되었으므로 이 둘을 반환한다
def goldbach(num):
# 이 변수 계산을 for문의 range 안에서 여러 수천, 수만 번 하지 않고 여기서 한번만 하는 것이 훨씬 효율적이다.
until = int(num/2) + 1
for i in range(2, until):
# 2를 제외하고 짝수면 소수가 아니기 때문에 다음으로 넘어간다.
# 함수 호출해서 검사하면 오버헤드 발생하기 때문에 그냥 여기서...
if i != 2 and i % 2 == 0:
continue
if i != 3 and i % 3 == 0:
continue
if isPrime(i) == False: # 2의 배수도, 3의 배수도 아닌 합성수이면 다음으로 넘어간다
continue
if isPrime(num - i): # 이제 i 는 소수이므로 num - i 의 결과가 소수이면 두 값을 리턴한다
return i, num - i
#print(goldbach(20))
#print(goldbach(200))
# my 선형방정식 해 구하기
# 893 * x - 2432* y = 266
# a * x - b * y = c
def linear_equation_compute(a, b, c):
g = gcd(a, b)
lst = []
if c % g != 0:
#print("정수해가 없습니다.")
return lst
count = 0
for i in range(1, b):
r = (a * i - c) % b
if r == 0:
y = (a * i - c) // b
lst.append((i, y))
#print(i, y)
count += 1
if count >= g:
break
return lst
# 선형방정식 해 구하기, 책에서 시키는 대로
# 3*x + 5*y = 1 (gcd(3, 5) = 1)
# (a, b) = (3, 5)
def linear_equation(a, b): # (3, 5)
g, w = a, b # g = 3, w = 5
x, v = 1, 0
while True: # w == 5 (w != 0)
if w == 0:
while x < 0:
x += b
y = (g - a*x) / b
#print(g, x, y) # g = 1, x = 2, y = (1 - 3 * 2) / 5 = -1 (1, 2, -1)
return x, y
break
# 나머지 구하는 항목
# t = 3 % 5 = 3
t = g % w # t = 3 % 5 = 3; t = 5 % 3 = 2; t = 3 % 2 = 1; t = 2 % 1 = 0
# 몫을 구하는 부분
# q = 5 // 3 = 1
q = g // w # q = 3 // 5 = 0; q = 5 // 3 = 1; q = 3 // 2 = 1; q = 2 // 1 = 2
# ?? 1 -몫*0
s = x - q*v # s = 1 - 0*0 = 1; s = 0 - 1 * 1 = -1; s = 1 - 1*(-1) = 2; s = -1 - 2*2 = -5
# x = v, g = w; x = 0, g = 5
x, g = v, w # x = 0, g = 5; x = 1, g = 3; x = -1, g = 2; x = 2, g = 1
v, w = s, t # v = 1, w = 3; v = -1, w = 2; v = 2, w = 1; v = -5, w = 0
#linear_equation(42, 30) # 6 -2 3
#linear_equation(10, 35) # 5 -3 1
# 합동방정식 해 구하기
# '친절한 수론 길라잡이' 65쪽 연습문제 8.7
# ax = c (mod m)
# 합동방정식의 해가 있다면 gcd(a, m) | c 이므로,
# gcd로 세 값을 모두 나눈 것을 da, gm, dc라 하면,
# da*x = dc (mod dm) 은 단 하나의 해를 가진다.
# 이 해를 x1 이라 하면, 해의 갯수는 모두 g개가 된다.
# 모든 해 = {x | x1 + dm * k, k는 0, 1, 2, ..., g-1}
# linear_equation_compute 도 같은 역할을 하는데, 단순하게 대입해서 계산하므로, 수가 크지면
# congruence_equation 함수가 훨씬 빠르게 작동한다.
def congruence_equation(a, c, m):
g = gcd(a, m) # 7
lst = []
if c % g != 0:
#print(f"gcd({a}, {m}) = {g} : {c}을 나눌 수 없습니다. 정수해가 없습니다.")
return lst
da = a // g # 3
dc = c // g # 2
dm = m // g # 13
for i in range(0, dm):
r = (da * i - dc) % dm
if r ==0:
# 정답은 i
#y = (a * i - c) // m # 정답은 (i, y) => (x, y)
break
for k in range(0, g):
x = i + dm*k
lst.append(x)
return lst
#lst = congruence_equation(943, 381, 2576)
#print("{}개 : {}".format(len(lst), lst))
# 다항식 합동식의 해를 함수를 이용해 구한다. 연습문제 8.8
# f(x) 가 x^4 + 5*x^3 + 4*x^2 - 6*x^1 - 4*x^0 이라면 차수와 계수를 딕셔너리 형태로 넘겨준다.
# m도 같이 넘겨주면, 아래 함수에서는 f(x)가 mod m 으로 0일 때를 계산해서 그 x값 리스트를 돌려준다.
# {4:1, 3:5, 2:4, 1:-6, 0:-4}
# 1* (x ** 4) + 5*(x**3) + 4*(x**2) + -6*(x**1) + -4*(x**0)
def congruence_equation_poly(poly, m):
# poly = {차수:계수, 차수:계수, 차수:계수....}
lst = []
for i in range(0, m): # 0 <= i < m 에서 ...
fx = 0 #함수값은 0으로 초기화하고 높은 차수부터 차례로 누적한다.
for d, c in poly.items(): # 딕셔너리 크기만큼 반복 d: 차수, c: 계수
fx += c * (i ** d) # fx 함수값을 계속 누적해서 구한다.
if fx % m == 0: # 함수값이 구해졌으면 mod m으로 합동인지 검사한다.
lst.append(i) # 합동이면 정답이므로, 목록에 추가한다.
return lst
# X^4 + 5X^3 + 4X^2 -6X -4 = 0 (mod 11) 0 <= X <11
#poly = {4:1, 3:5, 2:4, 1:-6, 0:-4}
#m = 11
#res = congruence_equation_poly(poly, m)
#print("m = {}: {}".format(m, res))
# 디리클레의 등차수열 정리를 확인하기 위해,
# 해당 정리를 만족하는 무수히 많은 소수가 있는지 찾아보는 함수
def dirichlet(a, m, max_num):
# p ≡ a (mod m) 인 p(소수)를 찾기
# 일단 소수를 먼저 찾아서 a (mod m)인지 확인하고 리스트에 저장
# max_num은 한도를 정함
lst = []
for i in range(2, max_num): # 2부터 max_num까지에서 소수를 찾음
if isPrime(i):
if i % m == a: # 소수가 나오면 m으로 나눠서 나머지가 a라면
lst.append(i) # 리스트에 추가
return lst
#res = dirichlet(3, 4, 200)
#print(res)
# 서로소인 숫자를 찾아 그 갯수를 반환하는 함수
def euler_phi(num):
lst = [1]
# 서로소인 숫자를 담는 리스트, 1은 기본
for i in range(2, num): # 2 ~ num-1
if gcd(i, num) == 1:
lst.append(i)
return len(lst)
# num 이하의 서로소인 숫자를 찾아 그 목록을 반환하는 함수
def euler_phi_list(num):
lst = [1]
# 서로소인 숫자를 담는 리스트, 1은 기본
for i in range(2, num): # 2 ~ num-1
if gcd(i, num) == 1:
lst.append(i)
return lst
# phi(m) = m(1-1/p1)(1-1/p2)...(1-1/pr) 공식을 이용해 오일러 파이값을 구하는 함수
# p1, p2... pr 은 m의 소인수
# 24 = 2^3 * 3^1 이고 p1 = 2, p2 = 3
def euler_phi2(m):
factorized_dict = factorize2(m) # num = 24일 경우, {2: 3, 3: 1}
phi_m = m
for key, _ in factorized_dict.items(): # key만 사용함
phi_m *=1 - 1/key
return int(phi_m)
# 인수분해하여 오일러 파이값을 구하는 함수
def euler_phi3(num):
factorized_dict = factorize2(num)
total = 1 # 파이값 즉, 서로소인 숫자들 누적 초기화: 계속 곱해야 하기 때문에 초기값은 1
for key, value in factorized_dict.items():
total *= key ** value - key ** (value - 1)
return total
def crt(b, m, c, n): # Chinese Remainder Theorem
# 예시)
# x = 1 (mod 3)
# x = 2 (mod 5)
mn = m * n
b_lst = []
for i in range(b, mn, m):
b_lst.append(i)
c_lst = []
for i in range(c, mn, n):
c_lst.append(i)
return set(b_lst) & set(c_lst)
#print(crt(1, 3, 2, 5))
#print(crt(3, 7, 5, 9))
#print(crt(3, 37, 1, 87))
#=============================================================
# 정해지지 않은 개수의 인자들을 받아 중국인 나머지 정리의 해를 구한다.
# 예: 5, 7, 2, 12, 8, 13
def crt2(*args): # 가변인자를 받는다
length = len(args) # 인자의 개수를 구한다 # 6
if length % 2 != 0:
print("인자의 개수는 짝수여야 합니다.")
return
r = []
m = []
multi_mn = 1 # 각 m을 모두 곱한 값이다. 이 값 이내에 해가 하나 있다.
for i in range(0, length, 2):
r.append(args[i])
m.append(args[i+1])
multi_mn *= args[i+1]
# mul_mn (= m1*m2*m3...)을 먼저 구해야 각 식의 해를 구할 수 있음. 위에서 구한 이유
answers = [] # 각 방정식마다 해를 구해 모은 집합
for i in range(len(r)):
answer = [i for i in range(r[i], multi_mn, m[i])] # 첫번째 방정식의 해 집합
answers.append(answer)
solution = set(answers[0].copy())
# 교집합을 누적으로 구해야 해서, 0번째 교집합을 해집합1로 둔다.
for i in range(len(answers)-1): # 교집합을 구한다. 해 집합이 2개일 경우, 그 교집합은 1개
solution = solution & set(answers[i+1])
return sorted(solution)
# 연습문제 11.5
#print(crt2(5, 7, 2, 12, 8, 13))
# 연습문제 11.6 <손자산경> 풀이
# x = 2 (mod 3), x = 3 (mod 5), x = 2 (mod 7)
#print(crt2(2, 3, 3, 5, 2, 7))
# 연습문제 11.7 계란을 실은 트럭
# x = 1 (mod 2), x = 1 (mod 3), x = 1 (mod 4), x = 1 (mod 5), x = 1 (mod 6), x = 0 (mod 7)
#print(crt2(1, 2, 1, 3, 1, 4, 1, 5, 1, 6, 0, 7))
# 일일이 나눠보는 방식으로 약수 구하는 함수, 속도가 느림
def getDivisor(num):
divisors = []
root_num = round(num ** 0.5)
for i in range(1, root_num+1):
if num % i == 0:
divisors.append(i)
divisors.append(num // i)
return sorted(list(set(divisors)))
# 소인수분해를 해서 약수 구하는 함수: 향상된 버전
# 숫자가 클 경우, 속도 면에서 이전 버전1 보다 훨씬 빠름
# 50234902430546664의 약수를 구했을 때 개선 버전이 1초 걸리고, 구버전이 35초 걸림.
def getDivisor2(num):
divisors = []
res = factorize2(num) # 예) {2:3, 3:2, 7:1} 형태
for key, value in res.items():
divisor = []
for i in range(value+1):
divisor.append(key**i)
divisors.append(divisor)
# divisors : [[2**0, 2**1, 2**2, 2**3], [3**0, 3**1, 3**2], [7**0, 7**1]]
# divisors : [[1, 2, 4, 8], [1, 3, 9], [1, 7]]
while True:
tmp = [] # 리스트 2개 원소를 서로 곱해서 약수를 만들어 임시 저장
# 리스트를 하나씩 줄여나가는데 리스트가 하나만 남았으면 그것이 약수 집합이므로 종료한다
if len(divisors) < 2:
break
# 리스트 2개의 각 원소를 서로 곱해서 약수를 만들어 저장
for n in divisors[0]:
for m in divisors[1]:
tmp.append(n * m)
# divisors 안의 리스트는 예에서 3개이고, 리스트[0]과 [1]을 곱해서 약수를 만들고(tmp),
# 기존 리스트[0]과 리스트 [1]을 삭제하고 1차 약수 리스트인 tmp를 리스트[0] 자리에 끼워넣는다.
# 한번 실행을 거치면 이제 리스트 목록은 총 2개가 되었다.
del divisors[0] # 원래 divisors[0] 삭제
del divisors[0] # 원래 divisors[1] 삭제
divisors.insert(0, tmp) # 생성된 tmp를 divisors 제일 앞에 끼워넣기
return sorted(divisors[0]) # 리스트의 리스트이므로 바깥 대괄호를 제거하고 정렬해서 반환
# 시그마 함수: 자신을 포함한 모든 약수의 총합을 돌려주는 함수
# 예) 12의 시그마 함수값: 1, 2, 3, 4, 6, 12 을 모두 더한 값
# '친절한 수론 길라잡이' 105~106에서 소개
def sigma_fun(num):
divisors = getDivisor2(num) # # 소인수분해 방식으로 약수를 모두 구한다
total = 0 # 합 초기화
for d in divisors: # 약수들을 모두 합한다
total += d
return total
# 자신을 제외한 약수들을 모두 더해서(진약수 합) 돌려준다
# 예) isPerfectNumber(28) : (True, 28)
def isPerfectNumber(num):
divisors = getDivisor2(num) # 소인수분방해 식으로 약수를 모두 구한다
divisors.remove(num) # 전체 약수에서 자기 자신을 제외
total = 0 # 합 초기화
for d in divisors: # 약수들을 모두 더한다
total += d
if total != num: # 약수들의 합이 자신이 아니라면 완전수가 아니다
return False, total
else: # 약수들의 합이 자신이라면 완전수다
return True, total
# 메르센 소수 구하기
# 2 ^ N - 1 = 소수 (인 경우 구하기)
def get_mersennePrime():
primes = [] # 소수들
m_primes = [] # 메르센 소수들
MAX_NUM = 200 # 여기까지 소수 구하기
# 일단 소수들 구하기
for i in range(2, MAX_NUM):
if isPrime(i):
primes.append(i)
print(primes, "\n")
for p in primes:
num = 2**p - 1
if isPrime(num):
m_primes.append((p, num))
print(p, num)
return m_primes
#get_mersennePrime()
# 유클리드 완전수 공식을 이용해 완전수 찾기
# 2^p - 1 이 소수면, 2^(p-1)(2^p-1)은 완전수다
# 2^p - 1 은 메르센 소수
def euclidPerfectNumber():
mersenne_prime = [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941]
for p in mersenne_prime:
num = (2**(p-1)) * (2**p-1)
print(num)
# 숫자가 곱셈완전수인지 판정하는 함수
# 소인수분해 방식으로 약수를 모두 구해 자신을 제외하고 모두 곱한다.
# 그 값이 자기 자신과 동일하면 곱셈완전수로 판정한다.
def isProductPerfect(num):
# 약수를 모두 구하기
divisors = getDivisor2(num) # 약수를 모두 구하기 (소인수분해 방식)
divisors.remove(num) # 전체 약수에서 자기 자신을 제외
gop = 1
for d in divisors: # 약수를 모두 곱한다.
gop *= d
if gop != num: # 원 수와 다르면
return False, gop
else: # 원 수와 같으면
return True, gop
# 약수의 개수만 구하는 함수
# 소인수분해해서 지수를 가지고 갯수를 구한다.
def divisorsCount(num):
divisors = factorize2(num) #{2:3, 3:2, 7:1} 형태
# 인수에는 관심없고 각기 지수에 1을 더해 모두 곱한다.
gop = 1
for _, e in divisors.items():
gop *= e + 1
return gop
# 두 수가 친화수인지 확인하는 함수
# m의 진약수의 합이 n이고 m의 진약수의 합이 n인 두 수 m, n을 친화수라 부른다.
# 예) 220의 진약수의 합: 284
# 284의 진약수의 합: 220
def isAmicablePair(m, n):
hap = m + n
# m의 전체 약수 합이 m+n이고, n의 전체 약수 합이 n + m이면 m ,n은 친화수이다.
return sigma_fun(m) == hap and sigma_fun(n) == hap
# 숫자 리스트를 건네주면 사교수인지 판정하는 함수. info 값을 True로 주면, 진행과정을 출력한다
def isSociableNumbers(numbers, info=False):
length = len(numbers) # 순환하는 이 수들의 개수가 위수(order)다, 물론 사교수라면.
isSociable = True # 사교수 여부를 체크하는 flag
for i, n in enumerate(numbers):
if i == length - 1: # 마지막 항목에 이르렀다면
i = -1 # 대조 항목의 인덱스인 i+1을 0(numbers의 첫번째 항목)으로 만들기 위해.
jinhap = sigma_fun(n) - n # 진합은 전체 약수 합에서 자신을 빼면 된다
if jinhap == numbers[i+1]: # 진합과 다음 항목값을 비교해서 동일하다면
if info:
print(f"{n}의 진약수의 합은 {numbers[i+1]}. 조건 만족!")
else: # 진합과 다음 항목값이 다르다면, 더 검사해볼 필요도 없이 루프를 중단해도 됨
if info:
print(f"{n}의 진약수의 합 {jinhap}가 {numbers[i+1]}와 다릅니다. 중단합니다!")
isSociable = False
break
if info:
if isSociable:
print("상기 수들은 사교수입니다.", f"위수는 {length}")
else: # False면
print("상기 수들은 사교수가 아닙니다.")
return isSociable
# 사교수의 범위를 정해주면 사교수를 찾는다
# max_try: 하나의 시작 수에 대해 진약수 합을 연이어 찾는 횟수
def getSociableNumbers(start, end, max_try=30):
for num in range(start, end):
jinhap = num # num을 다치지 않게 복사해서 사용한다.
lst = [] # 연속적인 진약수의 합(진합)을 저장할 리스트
lst.append(num) # num을 첫번째 원소로 저장한다
while True:
prev = jinhap # 첫 실행이면 jinhap은 num 이므로 prev 는 num 이 되고,
# 두 번째부터는 jinhap은 앞서 계산된 진약수의 합이므로,
# prev에는 앞서의 진합이 저장된다.
jinhap = sigma_fun(jinhap) - jinhap # 이전 진합(jinhap)을 넣어 새 진합(jinhap)을 구한다.
# 진합이 1(막다른 골목)이거나 이전 진합과 같은 값(그러면 완전수)이면 사교수가 아니므로 다음 수 검사로 넘어감
if jinhap < 2 and jinhap == prev:
break
# 이제 목록에 같은 값이 있는지 확인한다. 있다면
# 1) 완전수인 경우: 직전 진합과 진합이 같다. 이 경우는 위의 조건문에서 걸러졌음
# 2) 친화수인 경우: [220, 284]의 예라면, 284의 진합은 220이므로 220이 이미 목록(리스트)에 있다.
# 이때 목록에서 220의 인덱스는 0이므로 목록[0:]으로 슬라이싱하면 [220, 284]를
# 최종 예비 사교수 목록으로 얻을 수 있고 이의 위수는 len(목록)으로 2가 된다.
# 위수가 2라는 소리는, (넓은 의미의) 사교수에서 위수가 2이면 친화수가 된다는 뜻이고,
# 본격적인 사교수는 위수가 3일 때부터 취급한다는 뜻이다.
# 3) 사교수(위수 3이상)인 경우: 가령 [A, B, C, D, E, F, G](A, B,...G는 모두 숫자~)라는
# 진약수 합 목록이 있고 G의 진약수 합이 D라면, D의 인덱스는 3이다.
# 목록[3:]으로 슬라이싱하면 [D, E, F, G]라는 새로운 진합 목록을 얻는데, 이는 순환하는
# 사교수 목록이다. D->E->F->G->D->E.... 가 되니까.. 위수는 목록의 원소 개수로 4가 된다.
if jinhap in lst: # 구한 진약수 합이 목록에 이미 있다면
idx = lst.index(jinhap) # 해당 숫자의 인덱스를 찾아서
lst = lst[idx:] # 순환하는 부분만 잘라낸다
if len(lst) >= 3: # 위수가 해당값 이상이면 목적 사교수로 보고 출력한다.
print(num, lst)
break
else:
lst.append(jinhap) # 목록에 새로운 진약수 합과 같은 값이 없으므로 목록에 추가한다.
# 목록에 진합을 추가하다가 한참 뒤에 사교수(순환)의 첫 원소가 생길 수 있다.
# 가령 진합 목록이 [A, B, C, D, ....., Z, a, b, c, d, e, f, g](모두 숫자~) 이렇게 되어가는데,
# g의 진합이 d가 되어 d의 인덱스를 찾아내어 목록[인덱스:]로 [d, e, f, g]를 슬라이싱할 수 있다.
# 사교수 순환의 첫 원소인 d가 이렇게 한참 뒤에 출현할 수 있기 때문에 목록의 갯수에 약간
# 여유를 두는 게 좋다(찾고자 하는 위수의 2~3배 정도?). 목록의 개수를 지나치게 크게 하면 엉뚱한 데서 한없이 찾을 수 있으므로 주의.
if len(lst) >= max_try: break # 충분히 찾았는데도 사교수를 찾지 못했다면 그 숫자는 건너뛴다
# 연속제곱법
# '친절한 수론 길라잡이' 16장 연습문제 16.2
def successiveSquaring(a, k, m):
expo = [] # 지수 보관 예: [1, 2, 4, 8, 16, 32, 64, 128, 256]
c = [] # 계수 보관 예: [1, 1, 1, 0, 0, 0, 1, 0, 1]
count = int(math.log2(k)) + 1 # 2의 거듭제곱값 중에서 k값을 넘지 않는 마지막 수까지의 갯수
for n in range(0, count + 1):
if 2**n > k:
break
expo.append(2**n)
print(expo)
# [1, 2, 4, 8, 16, 32, 64, 128, 256]
length = len(expo) # 9
# bin_str = format(k, 'b') # 한줄로 표현할 수 있음.
bin_str = bin(k) # 0b101000111
bin_str = bin_str[2:] # 얖의 이진수 표시문자 '0b'를 제거. '0b101000111' -> '101000111'
bin_len = len(bin_str) # 9
if length != bin_len:
print("2진수 개수와 지수 개수가 다릅니다.")
print(f"count = {length}, bin_len = {bin_len}")
c = [0] * bin_len # 계수를 저장할 리스트 초기화 # [0, 0, 0, 0, 0, 0, 0, 0, 0]
for i in range(bin_len): # 0 ~ 8
# bin_str의 값들을 뒤집어 정수로 변환해서 c에 저장하기
# 뒤집는 이유는 지수 순서와 일치시켜 주 기 위해
c[bin_len-1-i] = int(bin_str[i])
print(c) # [1, 1, 1, 0, 0, 0, 1, 0, 1]
# c와 expo의 원소들을 한 쌍씩 곱한 다음, 차례로 더하면 k의 값이 나와야 한다. 점검.
total = 0
for i in range(bin_len):
total += c[i] * expo[i]
print(total) # 327
if total != k:
print("계수와 2의 거듭제곱을 곱해서 모두 더했는데 k의 값과 다릅니다.")
mod_ks = []
# (mod k)를 했을 때 나온 값들(=나머지)을 저장.
# 예에서는 [7, 49, 695, 227, 349, 675, 123, 628, 298]
mod_k = a % m # 첫번째 원소로 '7 % m' 을 저장
mod_ks.append(mod_k)
for i in range(length - 1): # 총 8회만 mod_k 계산하면 됨
mod_k = mod_k ** 2 % m
mod_ks.append(mod_k)
print(mod_ks) # [7, 49, 695, 227, 349, 675, 123, 628, 298]
values = []
# 반환할 값들을 계산한다
for i in range(length): # 0~8
if c[i]: # 계수가 0이 아닌 (즉 1인) 지수의 mod_k만 뽑아낸다.
values.append(mod_ks[i])
print(values) # [7, 49, 695, 123, 298]
# 이 값들을 모두 곱해서 (mod m)을 해주면 원하는 값이 나온다.
val_len = len(values) # 5
# 총 val_len - 1 회만큼 계산한다. 한꺼번에 전부 곱하면 숫자가 너무 커질 수 있기 때문에
# 차례로 둘씩 곱하면서 (mod m)을 해준다.
res = values[0]
for i in range(1, val_len):
print(res, values[i])
res = res * values[i] % m
return res # 최종 값: a**k (mod m)
# 연속제곱법, 향상된 버전
# '친절한 수론 길라잡이 연습문제 16.2.b에 나온 코드
def successiveSquaring2(a, k, m):
b = 1
while k >= 1:
if k % 2 == 1:
b = a*b % m
a = a**2 % m
k = k // 2
return b # a^k (mod m)
# 연속제곱법, 쑤튜브 공식 따라서
def my_successiveSquaring(a, k, m):
# 일단 지수의 이진수 표현을 얻기 k => '100011011'
bin_str = format(k, 'b')
val = 1
for i in bin_str:
c = int(i) # 문자열 '1'이나 '0'을 정수로 변환
val = (val**2 * a**c) % m # a^0 혹은 a^1을 곱함
return val
# 고차 합동식을 푼다 (17장: phi_m 구하면서 3단계식 풀이)
# x^k = b (mod m) x의 값은?
def high_congruence_equation(k, b, m):
# 1. phi_m을 구한다. 소인수분해가 안 되면 이 값을 못 구한다.
phi_m = euler_phi2(m)
# 2. 일차 합동방정식을 푼다.
# k*u - phi_m*v = 1
u, _ = linear_equation(k, phi_m)
# 3. 연속제곱법으로 합동식을 푼다.
# b^u (mod m) = x
x = my_successiveSquaring(b, u, m)
return x
# k, b, m = 329, 452, 1147
# print(high_congruence_equation(k, b, m))
# k, b, m = 113, 347, 463
# print(high_congruence_equation(k, b, m))
# k, b, m = 275, 139, 588
# print(high_congruence_equation(k, b, m))
# 암호책: 문자를 숫자로 바꾼다. 원래의 것에 공백과 문장부호 추가
cipher_book = {'A':11, 'B':12, 'C':13, 'D':14, 'E': 15, 'F':16, 'G':17, 'H':18, 'I':19, 'J':20,
'K':21, 'L': 22, 'M':23, 'N':24, 'O':25, 'P':26, 'Q':27, 'R':28, 'S':29, 'T':30,
'U':31, 'V':32, 'W':33, 'X':34, 'Y':35, 'Z':36, ' ':37, '.':38, ',':39, "'":40, '!':41, '~':42}
# key:value 인 암호북으로 value: key 암호북을 새로 생성
cipher_book2 = {}
for key, value in cipher_book.items():
cipher_book2[value] = key
#print(cipher_book2)
def RSA_incode(msg, m, k):
msg = msg.upper() # 대문자로 바꾼다
len_m = len(str(m)) # 9
len_ = len_m - 1 # 8
new_msg = "" # 메시지를 숫자로 바꾼 문자열
# 문자를 숫자로 변환
for c in msg: #'TOBEORNOTTOBE'
new_msg += str(cipher_book[c])
#print(new_msg) # "30251215252824253030251215"
len_new_msg = len(new_msg) # 26
howmany = math.ceil(len_new_msg / (len_m - 1)) # 배열을 쪼갠 갯수 round(26/(9-1))+1 = round(3.25)+1 = 4
div_new_msg = []
for i in range(howmany): # 0~3
if i == howmany - 1: # 마지막이라면
div_new_msg.append(int(new_msg[i*len_:]))
else:
div_new_msg.append(int(new_msg[i*len_:i*len_+len_]))
#print(div_new_msg) #[30251215, 25282425, 30302512, 15]
ciphered_list = []
for i in range(howmany): # 0~3
#ciphered = div_new_msg[i] ** k % m # 단순한 계산 방식. 수가 클 경우 시간이 많이 걸림.
ciphered = my_successiveSquaring(div_new_msg[i], k, m)
ciphered_list.append(ciphered)
return ciphered_list # [149419241, 62721998, 118084566, 40481382]
# 암호 풀기
def RSA_decode(ciphered_list, m, k, p, q):
# 1. phi_m 을 구한다. 소수 p와 q값을 알고 있기에 쉽게 구함
phi_m = (p-1)*(q-1)
# 2. k*u - phi_m*v = 1 1차 방정식을 푼다
u, _ = linear_equation(k, phi_m)
# 3. b^u (mod m) = x : 연속제곱법으로 구한다
answers = ""
for ciphered in ciphered_list:
answer = my_successiveSquaring(ciphered, u, m) # 어떤 숫자
answers += str(answer)
#print("---------", answers) # '30251215252824253030251215'
count = len(answers) // 2
restored = "" # 복원된 암호문 저장
for i in range(count):
tmp = int(answers[i*2:i*2+2])
#print(i, tmp)
restored += cipher_book2[tmp]
return restored # 복원된 암호문
'''
msg = "I loved her. I didn't know that at that time. Now I regret so much."
ciphered_lst = RSA_incode(msg, m, k)
print(ciphered_lst)
word = RSA_decode(ciphered_lst, m, k, p, q)
print(word)
'''
# 카마이클 수인지 정의대로 판정하는 함수
# 일일이 검사하므로 시간이 많이 걸림
def isCarmichael(m):
if isPrime(m) == True: # 소수는 카마이클 수가 아님.
return False
for a in range(m):
if a ** m % m != a:
return False #print("a는 합성수이지만 카마이클 수는 아닙니다.")
return True
# 소인수분해하여 소인수들에 대해서만 mod n을 검사한다. 따라서 속도가 빠르다.
def isCarmichael2(n):
if isPrime(n) == True: # 소수는 카마이클 수가 아님.
return False
res = factorize2(n) # 소인수분해 결과를 딕셔너리로 받는다
lst = []
for key, value in res.items(): # {7: 1, 19: 1, 67: 1}
if value != 1:
return False # 소인수의 지수가 1이 아니라면 카마이클 수가 아님.
lst.append(key)
for a in lst: # [7, 19, 67]
if a ** n % n != a: # 7 ^ n % n = 7 임을 확인.
return False
return True # 조건을 만족시켰으므로 카마이클 수임
# 카마이클 수 관련해서 m이 합성수라는 증인 목록을 만든다.
# 증인 갯수, 증인 비율, (증인 목록)을 반환한다.
def compositeWitness(m):
witnesses = []
if isPrime(m) == True: # 소수는 카마이클 수가 아님.
return witnesses
for a in range(2, m):
if a ** m % m != a:
witnesses.append(a)
return len(witnesses), len(witnesses)/m, witnesses
'''
for i in range(3, 21):
print(i, compositeWitness(i))
numbers = [287, 190, 314, 586, 935, 808, 728, 291]
for n in numbers:
print(n, compositeWitness(n)[0:2]) # 증인 목록이 너무 길어서 목록 출력은 제외
num = [561, 1105, 1729, 2465, 2821, 6601, 8911, 62745, 423899]
for n in num:
print(n, isCarmichael2(n))
'''