#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
additions to math standard library
"""
__author__ = "Philippe Guglielmetti"
__copyright__ = "Copyright 2012, Philippe Guglielmetti"
__credits__ = ["https://github.com/tokland/pyeuler/blob/master/pyeuler/toolset.py",]
__license__ = "LGPL"
import operator
from math import sqrt, log, log10, ceil
from itertools import count,izip_longest, ifilter
from itertools import combinations, permutations, product as cartesian_product
from itertools2 import drop, ireduce, groupby, ilen, compact, flatten
from decorators import memoize
import fractions
[docs]def lcm(a,b):
"""least common multiple"""
return abs(a * b) / fractions.gcd(a,b) if a and b else 0
#vector operations
[docs]def accsum(it):
"""Yield accumulated sums of iterable: accsum(count(1)) -> 1,3,6,10,..."""
return drop(1, ireduce(operator.add, it, 0))
cumsum=accsum #numpy alias
[docs]def dot(a,b):
"""dot product"""
try: #vector*vector
return sum(map( operator.mul, a, b))
except:
pass
try: #matrix*vector
return [dot(line,b) for line in a]
except:
pass
#matrix*matrix
res=[dot(a,col) for col in zip(*b)]
return map(list,zip(*res))
[docs]def product(nums):
"""Product of nums"""
return reduce(operator.mul, nums, 1)
def _longer(a,b,fillvalue=0):
'''makes a as long as b by appending fillvalues'''
n=len(b)-len(a)
if n>0:
a.extend([fillvalue]*n)
[docs]def vecadd(a,b,fillvalue=0):
"""addition of vectors of inequal lengths"""
args=izip_longest(a,b,fillvalue=fillvalue)
return [a+b for a,b in args]
[docs]def vecsub(a,b,fillvalue=0):
"""substraction of vectors of inequal lengths"""
return [ai-bi for ai,bi in izip_longest(a,b,fillvalue=fillvalue)]
[docs]def vecmul(a,b):
"""product of vectors of inequal lengths"""
return [a[i]*b[i] for i in range(min([len(a),len(b)]))]
[docs]def vecdiv(a,b):
"""quotient of vectors of inequal lengths"""
return [a[i]/b[i] for i in range(min([len(a),len(b)]))]
[docs]def veccompare(a,b):
"""compare values in 2 lists. returns triple number of paris where [a<b, a==b, a==c]"""
res=[0,0,0]
for i in range(min([len(a),len(b)])):
if a[i]<b[i]:res[0]+=1
elif a[i]==b[i]:res[1]+=1
else:res[2]+=1
return res
#stats
[docs]def mean(data):
"""mean of data"""
return sum(data)/len(data)
[docs]def variance(data,avg=None):
"""variance of data"""
if avg==None:
avg=mean(data)
s = sum(((value - avg)**2) for value in data)
var = s/(len(data) - 1)
return var
[docs]def stats(l):
"""returns min,max,sum,sum2,avg,var of a list"""
lo=float("inf")
hi=float("-inf")
n=0
sum=0
sum2=0
for i in l:
if i is not None:
n+=1
sum+=i
sum2+=i*i
if i<lo:lo=i
if i>hi:hi=i
if n>0:
avg=float(sum)/n
var=float(sum2)/n-avg*avg #mean of square minus square of mean
else:
avg=None
var=None
return lo,hi,sum,sum2,avg,var
# numbers functions
# mostly from https://github.com/tokland/pyeuler/blob/master/pyeuler/toolset.py
[docs]def fibonacci():
"""Generate fibonnacci serie"""
get_next = lambda (a, b), _: (b, a+b)
return (b for (a, b) in ireduce(get_next, count(), (0, 1)))
[docs]def factorial(num):
"""Return factorial value of num (num!)"""
return product(xrange(2, num+1))
[docs]def is_integer(x, epsilon=1e-6):
"""Return True if the float x "seems" an integer"""
return (abs(round(x) - x) < epsilon)
[docs]def int_or_float(x, epsilon=1e-6):
return int(x) if is_integer(x, epsilon) else x
[docs]def divisors(n):
"""Return all divisors of n: divisors(12) -> 1,2,3,6,12"""
all_factors = [[f**p for p in range(fp+1)] for (f, fp) in factorize(n)]
return (product(ns) for ns in cartesian_product(*all_factors))
[docs]def proper_divisors(n):
"""Return all divisors of n except n itself."""
return (divisor for divisor in divisors(n) if divisor != n)
[docs]def is_prime(n):
"""Return True if n is a prime number (1 is not considered prime)."""
if n < 3:
return (n == 2)
elif n % 2 == 0:
return False
elif any(((n % x) == 0) for x in xrange(3, int(sqrt(n))+1, 2)):
return False
return True
[docs]def get_primes(start=2, memoized=False):
"""Yield prime numbers from 'start'"""
is_prime_fun = (memoize(is_prime) if memoized else is_prime)
return ifilter(is_prime_fun, count(start))
[docs]def digits_from_num_fast(num):
"""Get digits from num in base 10 (fast implementation)"""
return map(int, str(num))
[docs]def digits_from_num(num, base=10):
"""Get digits from num in base 'base'"""
def recursive(num, base, current):
if num < base:
return current+[num]
return recursive(num/base, base, current + [num%base])
return list(reversed(recursive(num, base, [])))
[docs]def str_base(num, base, numerals = '0123456789abcdefghijklmnopqrstuvwxyz'):
"""string representation of an ordinal in given base"""
if base < 2 or base > len(numerals):
raise ValueError("str_base: base must be between 2 and %i" % len(numerals))
if num == 0:
return '0'
if num < 0:
sign = '-'
num = -num
else:
sign = ''
result = ''
while num:
result = numerals[num % (base)] + result
num //= base
return sign + result
[docs]def num_from_digits(digits, base=10):
"""Get digits from num in base 'base'"""
return sum(x*(base**n) for (n, x) in enumerate(reversed(list(digits))) if x)
[docs]def is_palindromic(num, base=10):
"""Check if 'num' in base 'base' is a palindrome, that's it, if it can be
read equally from left to right and right to left."""
digitslst = digits_from_num(num, base)
return digitslst == list(reversed(digitslst))
[docs]def prime_factors(num, start=2):
"""Return all prime factors (ordered) of num in a list"""
candidates = xrange(start, int(sqrt(num)) + 1)
factor = next((x for x in candidates if (num % x == 0)), None)
return ([factor] + prime_factors(num / factor, factor) if factor else [num])
[docs]def factorize(num):
"""Factorize a number returning occurrences of its prime factors"""
return ((factor, ilen(fs)) for (factor, fs) in groupby(prime_factors(num)))
[docs]def greatest_common_divisor(a, b):
"""Return greatest common divisor of a and b"""
return (greatest_common_divisor(b, a % b) if b else a)
[docs]def least_common_multiple(a, b):
"""Return least common multiples of a and b"""
return (a * b) / greatest_common_divisor(a, b)
[docs]def triangle(x):
"""The nth triangle number is defined as the sum of [1,n] values. http://en.wikipedia.org/wiki/Triangular_number"""
return (x*(x+1))/2
[docs]def is_triangle(x):
return is_integer((-1 + sqrt(1 + 8*x)) / 2)
[docs]def pentagonal(n):
return n*(3*n - 1)/2
[docs]def is_pentagonal(n):
return (n >= 1) and is_integer((1+sqrt(1+24*n))/6.0)
[docs]def hexagonal(n):
return n*(2*n - 1)
[docs]def get_cardinal_name(num):
"""Get cardinal name for number (0 to 1 million)"""
numbers = {
0: "zero", 1: "one", 2: "two", 3: "three", 4: "four", 5: "five",
6: "six", 7: "seven", 8: "eight", 9: "nine", 10: "ten",
11: "eleven", 12: "twelve", 13: "thirteen", 14: "fourteen",
15: "fifteen", 16: "sixteen", 17: "seventeen", 18: "eighteen",
19: "nineteen", 20: "twenty", 30: "thirty", 40: "forty",
50: "fifty", 60: "sixty", 70: "seventy", 80: "eighty", 90: "ninety",
}
def _get_tens(n):
a, b = divmod(n, 10)
return (numbers[n] if (n in numbers) else "%s-%s" % (numbers[10*a], numbers[b]))
def _get_hundreds(n):
tens = n % 100
hundreds = (n / 100) % 10
return list(compact([
hundreds > 0 and numbers[hundreds],
hundreds > 0 and "hundred",
hundreds > 0 and tens and "and",
(not hundreds or tens > 0) and _get_tens(tens),
]))
# This needs some refactoring
if not (0 <= num < 1e6):
raise ValueError, "value not supported: %s" % num
thousands = (num / 1000) % 1000
strings = compact([
thousands and (_get_hundreds(thousands) + ["thousand"]),
(num % 1000 or not thousands) and _get_hundreds(num % 1000),
])
return " ".join(flatten(strings))
[docs]def is_perfect(num):
"""Return -1 if num is deficient, 0 if perfect, 1 if abundant"""
return cmp(sum(proper_divisors(num)), num)
[docs]def number_of_digits(num, base=10):
"""Return number of digits of num (expressed in base 'base')"""
return int(log(num)/log(base)) + 1
[docs]def is_pandigital(digits, through=range(1, 10)):
"""Return True if digits form a pandigital number"""
return (sorted(digits) == through)
#norms and distances
[docs]def sets_dist(a,b):
"""http://stackoverflow.com/questions/11316539/calculating-the-distance-between-two-unordered-sets"""
c = a.intersection(b)
return sqrt(len(a-c)*2 + len(b-c)*2)
[docs]def sets_levenshtein(a,b):
"""levenshtein distance on sets
@see: http://en.wikipedia.org/wiki/Levenshtein_distance
"""
c = a.intersection(b)
return len(a-c)+len(b-c)
[docs]def levenshtein(seq1, seq2):
"""return http://en.wikipedia.org/wiki/Levenshtein_distance distance between 2 iterables
http://en.wikibooks.org/wiki/Algorithm_Implementation/Strings/Levenshtein_distance#Python """
oneago = None
thisrow = range(1, len(seq2) + 1) + [0]
for x in xrange(len(seq1)):
twoago, oneago, thisrow = oneago, thisrow, [0] * len(seq2) + [x + 1]
for y in xrange(len(seq2)):
delcost = oneago[y] + 1
addcost = thisrow[y - 1] + 1
subcost = oneago[y - 1] + (seq1[x] != seq2[y])
thisrow[y] = min(delcost, addcost, subcost)
return thisrow[len(seq2) - 1]
#combinatorics
[docs]def ncombinations(n, k):
"""Combinations of k elements from a group of n"""
return cartesian_product(xrange(n-k+1, n+1)) / factorial(k)
[docs]def combinations_with_replacement(iterable, r):
"""combinations_with_replacement('ABC', 2) --> AA AB AC BB BC CC"""
pool = tuple(iterable)
n = len(pool)
for indices in cartesian_product(range(n), repeat=r):
if sorted(indices) == list(indices):
yield tuple(pool[i] for i in indices)