사용자:토끼군/정수스크립트

수 목록의 정수 페이지들에서 "수학적 성질" 부분을 채워 넣을 때 쓸 수 있는 간단한 파이썬 스크립트입니다. 이 스크립트는 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