"""Simple tools for manipulating types linear over the positive integers.
Example linear spaces over the positive integers:
the natural numbers
the integers, whether modulo some value or not
the polynomials with integer coefficients
the rational, real or complex numbers
any vector or linear space over any of the above
"""
# Modular division (where possible, e.g. prime base).
def dividemod(num, den, base):
"""Division modulo a base.
Returns a q for which num - q*den is a multiple of base, if possible: raises
ValueError on failure."""
__algorithm = """
First reduce num and den mod base.
If num is 0, q is 0; if den is 0, we have ValueError.
Thereafter, to solve for n-q*d = s*b:
wlog, q = q%b and s = s%d
rearrange to n%d - s * (b%d) = q*d
i.e. s = (n%d) / (b%d) mod d
use this in n - s*b = q*d to discover q = n / d mod b
This may be implemented recursively as
q, r = divmod(n - dividemod(n%d, b%d, d) * b, d)
assert r is 0
return q
However, we can unroll the recursion to a pair of while loops.
Compare and contrast with Euclid's algorithm (below). """
stack, q = [], 0
n, d, b = num % base, den % base, base
try:
while n:
q, r = divmod(n, d)
assert 0 < n < b > d > 0
if r: stack.append((n, d, b))
assert r <= n and b % d < d and d < b
n, d, b = r, b % d, d
except ZeroDivisionError:
raise ValueError, '%d / %d mod %d' % (num, den, base)
while stack:
n, d, b = stack.pop()
q, r = divmod(n - (q % d) * b, d)
assert r is 0
return q % base
def gcd(a, b):
"""Pair-wise highest common factor.
The value returned is, strictly, that value whose set of factors is the
intersection of the sets of factors of the two arguments, ignoring all
universal factors (e.g. 1, -1: values which are factors of everything).\n"""
# Any negative factor's matching positive is also a factor, and is
# greater than any negative.
if b < 0: b = -b
if a < 0: a = -a
elif a == 0: return b # gcd(0,b) = abs(b)
# gcd(a,0) falls out naturally in the following
# Euclid's algorithm:
while b > 0: a, b = b, a % b
return a
def hcf(*others):
"""The highest natural common factor of its arguments.
All arguments should be members of a linear space over the natural numbers,
eg (optionally long) integers or polynomials with such coefficients; there
may be arbitrarily many arguments.
Formally, the % operator, as defined for the given arguments, must be
guaranteed to yield a positive or zero value whenever its (two) arguments
are positive. This is true for integers.
Formally, at least one argument (to hcf) must be non-zero: if (there are no
arguments, or) the arguments are all zero, of which every value is a factor,
there is no highest common factor - all integers are factors of 0. In this
case, the value zero is returned.
To justify not raising an error when (there are no inputs, or) all inputs
are zero, we can re-cast the definition as: `that non-negative value which
has, as its factors, exactly those naturals which are factors of all the
arguments' (reading 'is n a factor of all arguments' as the negation of 'no
argument does not have n as a factor' for the case of no arguments). So the
result must be a multiple of every common factor of the arguments, but of
nothing else. This coincides with `highest common factor' when at least one
argument is non-zero: and gives zero when all arguments are zero. Note that
-1 is a factor of every value (just as is 1), hence the need to specify
`non-negative'.
With this (undeniably less catchy) re-definition, we also get: a
concatenation of lists of values has, as its hcf, the hcf of the values
obtained by taking the hcfs of the lists seperately; i.e. hcf is a
transitive binary operator.\n"""
this = 0
for other in others: this = gcd(this, other)
return this
def lcm(*others):
"""The smallest common multiple of its arguments.
All arguments should be members of a linear space over the natural numbers,
e.g. (optionally long) integers or polynomials with such coefficients.
There may be arbitrarily many arguments.
If any entry is zero, so is the result: zero is a multiple of everything,
and is smaller than any other value; furthermore, nothing else is a multiple
of zero, so it is the only candidate. The return when no arguments are
supplied is 1, to ensure that lcm is transitive.\n"""
this = 1
for other in others:
if not other: return other # in case it's an object whose class views it as zero
c = gcd(this, other) # > 0, as other != 0.
this = other * this / c
# this's sign is currently the product of the signs of the arguments.
if this < 0: return -this
return this
theorem = """Any rational whose square is an integer is, itself, an integer.
As a special case, this tells us that the square root of 2 is irrational.
Proof:
Suppose, for positive integers p, q, that the square of p/q is an integer, n.
Thus p.p = q.q.n and n is positive. There are positive integers m, r for
which m.r.r = n and m has no perfect square as a factor, yielding p.p =
(q.r).(q.r).m.
Expressing p, q.r and m in terms of their prime factors we now find that every
prime factor of m has odd multiplicity as a factor of p.p, all of whose
factors have even multiplicity; thus m cannot have any prime factors, so m is
1 and p = q.r has q as a factor so p/q = r is a positive integer. """
_early_primes = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
def factorsum(N):
"""Returns the sum of the proper factors of N.
Thus N is perfect precisely if N == factorsum(N).
"""
# Find (upper bound on) highest proper factor of N:
for p in _early_primes:
i, r = divmod(N, p)
if r == 0: break
S = 0 # sum proper factors:
while i > 0:
if N % i == 0:
S += i
i -= 1
return S
def perfect():
"""Returns an iterator over the perfect numbers.
All known perfect numbers are 2**n * (2**(1+n) -1) for some positive natural
n (for which the second factor, 2**(1+n) -1, is prime), so iterating over
these (which is *much* quicker, even pausing to check primality of the
second factor) shall (probably) yield the same result. The iterator yielded
by this implementation checks for other perfect numbers and prints a message
if it ever finds one.\n"""
i = 1
while True:
# functionally, call factorsum(i): but abort if > i.
for p in _early_primes:
j, r = divmod(i, p)
if r == 0: break
S = 0
while j > 0 and S <= i:
if i % j == 0: S += j
j -= 1
if i == S: # iff i == factorsum(i)
yield i
# Check to see if i fits the usual pattern:
n = 0 # max power of two that is a factor of i.
while not (i & (1L << n)): n += 1
if (i >> n) != ((1L << (1+n)) -1):
print 'Unusually perfect', hex(i), n, hex(i >> (n-1))
i += 1
def Collatz(n):
"""Iterator for the Collatz conjecture's sequence for n.
It is conjectured that, whatever positive integer n you give to this
function, the resulting sequence shall ultimately terminate (by yielding 1).
It is known (by experiment) that the conjecture is good up to n = 10 * 2**58
(and, hence, also good for any n which is just a power of 2 times some
positive integer less than this limit).
The function iterated is, with Z+ = {positive integers}, the union of (: n
← 2.n :Z+) and (: 6.j+4 ←2.j+1 :Z+). The conjecture effectively
says that its transitive closure subsumes ({1}:|Z+).\n"""
yield n
while n != 1:
if n % 2: n = 3 * n + 1
else: n = n / 2
yield n
raise StopIteration
def sqrt(val):
"""Returns the highest natural whose square does not exceed val"""
if val < 0: # Every natural's square exceeds val.
raise ValueError('Negative value has no square root', val)
v, bit = val, 0
while v:
v >>= 2
bit += 1
bb = 1 << (2 * bit)
while bit:
bb >>= 2
up = (v << bit) + bb
bit -= 1
if up <= val:
v |= 1 << bit
val -= up
return v
# and now for something python-2.2-specific:
class Naturals (list):
class Suc (dict):
def __init__(self, bok=None):
self.clear()
if bok is not None:
self.update(bok)
self[bok] = self
def suc(self): return self.__class__(self)
Suc.__hash__ = Suc.__len__
def __init__(self, zero=Suc()): self[:] = [zero]
del Suc
__upget = list.__getitem__
def __getitem__(self, ind):
while ind >= len(self): self.append(self[-1].suc())
return self.__upget(ind)
naturals = Naturals()
del Naturals
# NB: len(str(naturals[1+n])) = 13 * 3**n -2