사용자:토끼군/정수스크립트
위키백과 ― 우리 모두의 백과사전.
< 사용자:토끼군
수 목록의 정수 페이지들에서 "수학적 성질" 부분을 채워 넣을 때 쓸 수 있는 간단한 파이썬 스크립트입니다. 이 스크립트는 GNU LGPL 라이센스로 공개됩니다. --토끼군
[편집] numtheory.py
"""numtheory.py -- python implementation of number theory functions.
Kang Seonghoon (Tokigun) @ TokigunStudio 2005.
Distributable under the GNU LGPL.
"""
def issprp(n, a):
"""Returns True if n is a-SPRP, False otherwise."""
if n < 2 or n & 1 == 0: return False
d = n - 1
while d & 1 == 0:
d >>= 1
if pow(a, d, n) + 1 == n: return True
return pow(a, d, n) == 1
def isprime(n):
"""Returns True if n < 341,550,071,728,321 is prime, False otherwise."""
if n < 2: return False
elif n < 4: return True
elif not issprp(n, 2): return False
elif n < 2047: return True
elif not issprp(n, 3): return False
elif n < 1373653: return True
elif not issprp(n, 5): return False
elif n < 25326001: return True
elif n == 3215031751 or not issprp(n, 7): return False
elif n < 118670087467: return True
elif not issprp(n, 11): return False
elif n < 2152302898747: return True
elif not issprp(n, 13): return False
elif n < 3474749660383: return True
elif not issprp(n, 17): return False
elif n < 341550071728321: return True
else: raise ValueError, "cannot test primality of n >= 341,550,071,728,321."
def factor_g():
"""Generator for factorization."""
yield 2; yield 3; k = 6
while True: yield k-1; yield k+1; k += 6
def factorize(n):
"""Factorize given number, returns {p1:k1, p2:k2, ...} for n=p1^k1*p2^k2*..."""
result = {}
gen = factor_g()
k = gen.next()
while n > 1:
while k < n and n % k != 0: k = gen.next()
try: result[k] += 1
except: result[k] = 1
n /= k
return result
def unfactorize(n):
"""Returns original value of factorized dictionary."""
result = 1
for p, k in n.iteritems():
result *= p ** k
return result
def phi(n):
"""Euler's Totient Function (Phi Function)."""
if not isinstance(n, dict): n = factorize(n)
result = 1
for p, k in n.iteritems():
result *= (p - 1) * p**(k-1)
return result
def divisor(n, a):
"""Divisor Function."""
if not isinstance(n, dict): n = factorize(n)
result = 1
if a == 0:
for p, k in n.iteritems():
result *= k + 1
else:
for p, k in n.iteritems():
result *= (p**(a*(k+1)) - 1) / (p**a - 1)
return result
def tau(n):
"""Tau Function, special case of divisor function for a=0."""
return divisor(n, 0)
def sigma(n):
"""Sigma Function, special case of divisor function for a=1."""
return divisor(n, 1)
def unitarydivisor(n, a=1):
"""Unitary Divisor Function."""
if not isinstance(n, dict): n = factorize(n)
result = 1
for p, k in n.iteritems():
result *= 1 + p**(a*k)
return result
def moebius(n):
"""Moebius' Mu Function."""
if not isinstance(n, dict): n = factorize(n)
if [x for x in n.values() if x>1] != []: return 0
return len(n) % 2 == 0 and 1 or -1
def mertens(n):
"""Mertens Function."""
if isinstance(n, dict): n = unfactorize(n)
try:
if n < len(mertens.cache):
return mertens.cache[n]
except: mertens.cache = [0]
ii = len(mertens.cache)
result = mertens.cache[ii-1]
for i in xrange(ii, n+1):
result += moebius(i)
mertens.cache.append(result)
return result
if __name__ == '__main__':
import sys
num = []
for x in sys.argv[1:]:
if '..' in x:
start, done = x.split('..')
num += range(int(start), int(done)+1)
else:
num.append(int(x))
for x in num:
if isprime(x):
y = {x: 1}
print x, "is prime."
else:
y = factorize(x)
print x, "is composite:", x, "=",
print " * ".join([k>1 and "%d^%d"%(p,k) or str(p) for p,k in y.iteritems()])
print "\tphi: %d, tau(sigma0): %d, sigma(sigma1): %d" % (phi(y), tau(y), sigma(y))
print "\tmoebius: %d, unitarydivisor1: %d, mertens: %d" % (moebius(y), unitarydivisor(y), mertens(x))
[편집] 사용법
다음과 같이 뒤에 계산하고자 하는 숫자를 넣으면 실행할 수 있습니다. ".."가 들어 가면 범위를 나타냅니다.
$ python numtheory.py 30..35 40..42 50
30 is composite: 30 = 2 * 3 * 5
phi: 8, tau(sigma0): 8, sigma(sigma1): 72
moebius: -1, unitarydivisor1: 72, mertens: -3
31 is prime.
phi: 30, tau(sigma0): 2, sigma(sigma1): 32
moebius: -1, unitarydivisor1: 32, mertens: -4
32 is composite: 32 = 2^5
phi: 16, tau(sigma0): 6, sigma(sigma1): 63
moebius: 0, unitarydivisor1: 33, mertens: -4
33 is composite: 33 = 11 * 3
phi: 20, tau(sigma0): 4, sigma(sigma1): 48
moebius: 1, unitarydivisor1: 48, mertens: -3
34 is composite: 34 = 17 * 2
phi: 16, tau(sigma0): 4, sigma(sigma1): 54
moebius: 1, unitarydivisor1: 54, mertens: -2
35 is composite: 35 = 5 * 7
phi: 24, tau(sigma0): 4, sigma(sigma1): 48
moebius: 1, unitarydivisor1: 48, mertens: -1
40 is composite: 40 = 2^3 * 5
phi: 16, tau(sigma0): 8, sigma(sigma1): 90
moebius: 0, unitarydivisor1: 54, mertens: 0
41 is prime.
phi: 40, tau(sigma0): 2, sigma(sigma1): 42
moebius: -1, unitarydivisor1: 42, mertens: -1
42 is composite: 42 = 2 * 3 * 7
phi: 12, tau(sigma0): 8, sigma(sigma1): 96
moebius: -1, unitarydivisor1: 96, mertens: -2
50 is composite: 50 = 2 * 5^2
phi: 20, tau(sigma0): 6, sigma(sigma1): 93
moebius: 0, unitarydivisor1: 78, mertens: -3

