#!/usr/bin/python -Wall

# ================================================================
# These are some handy utility routines which I always want to
# be available in my software environment.
#
# John Kerl
# kerl.john.r@gmail.com
# 2007-03-16
# ================================================================
# This software is released under the terms of the GNU GPL.
# Please see LICENSE.txt in the same directory as this file.
# ================================================================

from __future__ import division # 1/2 = 0.5, not 0.
import sys
import math
import cmath

# ----------------------------------------------------------------
# Factorial function.
def fact(n):
	"""fact(n): Computes n factorial."""
	if (n == 0):
		return 1
	if (n <= 1):
		return n
	else:
		#return n * fact(n-1)
		prod = 1
		while n > 1:
			prod *= n
			n -= 1
		return prod

# ----------------------------------------------------------------
# Double factorial function.
def dfact(n):
	"""dfact(n): Computes n factorial."""
	if (n == 0):
		return 1
	if (n == -1):
		return 1
	if (n == 1):
		return 1
	return n * dfact(n-2)

# ----------------------------------------------------------------
# Binomial coefficients.
def binc(n, k):
	"""binc(n, k): Computes n choose k."""
	if (k > n): return 0
	if (k < 0): return 0
	if (k > int(n/2)):
		k = n - k

	rv = 1
	for j in range(0, k):
		rv *= n - j
		rv /= j + 1
	return int(rv)

# ----------------------------------------------------------------
def printf_row(row, fmt="%11.7f"):
	"""printf_row(list, optional format)"""
	fmt_with_space = " " + fmt
	for e in row:
		print fmt_with_space % (e),
	print

# ----------------------------------------------------------------
# Similar to Matlab's [start:step:end].
def mfrange(start, step, end):
	"""mfrange(start, step, end):  Similar to Matlab's [start:step:end]."""

	# If the arguments are already float or complex, this won't hurt.
	# If they're int, they will become float.
	start *= 1.0
	step  *= 1.0
	end   *= 1.0

	# Correctly handle the roundoff cases.  Do this instead of looping
	# while (current < end) since, depending on step value, that may
	# slightly undershoot or overshoot the end.
	n = int(round((end - start) / step)) + 1

	list = [None] * n
	current = start
	for i in xrange(0, n):
		list[i] = current
		current += step

	return list

# ----------------------------------------------------------------
# For reals or complexes (hurrah for polymorphism!):  linear from a to b.
def frange(start, end, n):
	"""frange(start, end, n):  Python's range(), for floats and complexes."""

	# If the arguments are already float or complex, this won't hurt.
	# If they're int, they will become float.
	start *= 1.0
	end   *= 1.0

	step = (end - start) / n
	list = []
	current = start
	for k in range(0, n):
		list.append(current)
		current += step
	return list

# ----------------------------------------------------------------
# Multiplicative from a to b.
# Example: mulfrange(1, 64, 6) gives [1., 2., 4., 8., 16., 32.]
def mulfrange(start, end, n):
	"""mulfrange(start, end, n):  multiplicative range() for floats."""

	# If the arguments are already float or complex, this won't hurt.
	# If they're int, they will become float.
	start *= 1.0
	end   *= 1.0

	step = (end / start) ** (1./n)
	list = []
	current = start
	for k in range(0, n):
		list.append(current)
		current *= step
	return list

# ----------------------------------------------------------------
def cphase(c):
	"""cphase(c):  returns complex phase: (-pi, +pi]."""
	return math.atan2(c.imag, c.real)

# ----------------------------------------------------------------
# frange-ish for complexes:  Re^i theta for fixed R and varying theta.
def ztrange(R, start_theta, end_theta, n, center=0+0j):
	"""ztrange(R, start_theta, end_theta, n):  R exp(i theta) for fixed R
	and varying theta."""

	# If the arguments are already float or complex, this won't hurt.
	# If they're int, they will become float.
	start_theta *= 1.0
	end_theta   *= 1.0

	list = []
	thetas = frange(start_theta, end_theta, n)
	for theta in thetas:
		list.append(center + R * cmath.exp(1j * theta))
	return list

# ----------------------------------------------------------------
# frange-ish for complexes:  Re^i theta for varying R and fixed theta.
def zrrange(start_R, end_R, theta, n):
	"""ztrange(start_R, end_R, theta, n):  R exp(i theta) for varying R
	and fixed theta."""

	# If the arguments are already float or complex, this won't hurt.
	# If they're int, they will become float.
	start_R *= 1.0
	end_R   *= 1.0

	list = []
	Rs = frange(start_R, end_R, n)
	for R in Rs:
		list.append(R * cmath.exp(1j * theta))
	return list

# ----------------------------------------------------------------
# Example:  rnumint(sin, frange(0, 2*pi, 1000)).
# Again, this works fine for complexes as well.

# No adaptive quadrature -- just sums along the passed-in mesh.

def numint(f, xs):
	"""numint(f, xarray):  Naive numerical integration for float and complex."""
	sum = 0
	n = len(xs)
	ys = map(f, xs)
	for i in range(0, n-1):
		sum += 0.5 * (ys[i] + ys[i+1]) * (xs[i+1] - xs[i])
	sum += ys[n-1] * (xs[n-1]-xs[n-2])
	return sum

# ----------------------------------------------------------------
# Use the left-hand rule at the left endpoint, the right-hand rule at the
# right endpoint, and the midpoint rule at the interior points.
def numderiv(f, xs):
	"""numderiv(f, xarray):  Numerical differentiation."""
	n = len(xs)
	list = []

	# f'(x) ~= [ f(x+h) - f(x) ] / h
	x0 = xs[0]; x1 = xs[1]
	y0 = f(x0); y1 = f(x1)
	list.append((y1 - y0) / (x1 - x0))

	for i in range(1, n-1):
		# f'(x) ~= [ f(x+h) - f(x-h) ] / 2h
		x0 = xs[i-1]; x1 = xs[i+1]
		y0 = f(x0); y1 = f(x1)
		list.append((y1 - y0) / (x1 - x0))

	# f'(x) ~= [ f(x) - f(x-h) ] / h
	x0 = xs[n-2]; x1 = xs[n-1]
	y0 = f(x0); y1 = f(x1)
	list.append((y1 - y0) / (x1 - x0))
	return list

# ----------------------------------------------------------------
def clist_to_rlist(cs):
	list = []
	for c in cs:
		list.append(c.real)
		list.append(c.imag)
	return list