def my_sigma(n, k): prod = 1 for p, e in factor(n): psum = 0 for i in range(0, e + 1): psum += p^(k * i) prod *= psum return prod def my_xgcd(n, m): if n == 0: return (m, 0, 1) if m >= 0 else (-m, 0, -1) else: g, x, y = my_xgcd(m % n, n) return (g, y - (m // n) * x, x) def my_euclidean(n, m): if abs(m) < abs(n): m, n = n, m r = m % n while r != 0: print "%d = %d * %d + %d" % (m, n, m // n, r) m, n = n, r r = m % n print "%d = %d * %d" % (m, n, m // n) def my_euler_phi(n): res = n for p, e in factor(n): res *= (1 - 1/p) return res def my_kronecker(a, n): prod = 1 for p, e in factor(a): if e % 2 == 0: continue else: prod *= kronecker(p, n) return prod def solve_pell(N, numTry = 100): cf = continued_fraction(sqrt(N)) for i in range(numTry): denom = cf.denominator(i) numer = cf.numerator(i) if numer^2 - N * denom^2 == 1: return numer, denom return None, None if __name__ == "__main__": # Chapter1 print "my_sigma test" print "my_sigma(12, 3) == sigma(12, 3): ", my_sigma(12, 3) == sigma(12, 3) print # Chapter2 print "my_xgcd test" g, a, b = my_xgcd(123, 54) print "my_xgcd(123, 54) = %d, %d, %d" % (g, a, b) print "a * 123 + b * 54 == g: ", a * 123 + b * 54 == g print print "my_euclidean(741, 715)" my_euclidean(741, 715) print # Chapter3 print "my_euler_phi test" print "my_euler_phi(24) == euler_phi(24)", my_euler_phi(24) == euler_phi(24) print # Chapter5 print "kronecker_multiplicity test" print "my_kronecker(48, 31) == kronecker(48, 31)",\ my_kronecker(48, 31) == kronecker(48, 31) print # Chapter8 print "solve_pell test" a, b = solve_pell(21) print "solve_pell(21) = (%d, %d)" % (a, b) print "%d^2 - 21 * %d^2 == 1" % (a, b), a^2 - 21 * b^2 == 1 print