Simple Recipes in Python

William Park

March 1999

This article describes few elementary numerical routines
that I wrote in Python. They are loosely modelled after *Numerical
Recipes in C* [1] because I needed, at the time,
actual source codes which I can examine instead of just wrappers
around Fortran libraries like **NumPy** [2] and
**Octave** [3]. As evident from the
documentations, the routines were written with emphasis on clarity
rather than on runtime efficiency.

I use TeX notation whenever HTML markup becomes too
complicated or impossible to do; so, this article should be readable
in both graphical and text browsers. Also, I use my own
**HTMLtag** package to generate LaTeX-like document structures,
such as title, section numbers, footnotes, table of contents, and
other HTML tags. To get an overview of this article, see *Table of
Contents* at the bottom.

Since all the source codes are included, you can save this
to a file, remove HTML, and use it as your own module --- it's called
**emath.py** on my machine. In the **Usage** section,
`int` denotes integer; `real` denotes integer or
float, especially when there is possibility of complex;
`number` denotes any integer, float, or complex.

**Summary**- Sinc function

**Usage**`real`= sinc(`real`)

The sinc() function usually comes up when doing Fourier
Transform, since sinc(*f*) is the fourier transform of
rect(*t*) = {1, if |*t*| < 1/2;
0, if |*t*| > 1/2}.

**Summary**- Cubic root

**Usage**`real`= cbrt(`real`)

GNU C **libm** library has cubic root, hyperbolic,
gamma, bessel, expontial, and other higher functions, but standard
Python **math** module does not include them. So, I implement
cbrt() using pow() which is defined as pow(*x*, n) = *x*^n
for only *x* >= 0.

**Summary**- Conversion between rectangular and polar coordinates

**Usage**-
`real`,`real`= rect(`real`,`real`[, deg=1])

`real`,`real`= polar(`real`,`real`[, deg=1])

Normally, rect() and polar() uses radian for angle; but,
if *deg*=1 specified (any non-zero will do), degree is used
instead.

**Summary**- Evaluate n'th degree polynomial

**Usage**`number`= polyeval(`list`,`number`)

Evaluating n'th degree polynomial is simple loop, starting with highest coefficient a[n]. I wrote this because the polynomial module supplied with previous Python releases was brain-dead.

**Summary**- Find the first derivative of polynomial

**Usage**`list`= polyderiv(`list`)

**Summary**- Factor out a root from n'th degree polynomial, and return the remaining (n-1)'th degree polynomial.

**Usage**`list`= polyreduce(`list`,`number`)

**Summary**- Solve quadratic equation with real coefficients

**Usage**`number`,`number`= quadratic(`real`,`real`[,`real`])

Normally, *x*^2 + a*x* + b = 0 is assumed with
the 2 coefficients as arguments; but, if 3 arguments are present, then
a*x*^2 + b*x* + c = 0 is assumed.

**Summary**- Solve cubic equation with real coefficients

**Usage**`real`,`number`,`number`= cubic(`real`,`real`,`real`[,`real`])

Normally, *x*^3 + a*x*^2 + b*x* + c = 0 is
assumed with the 3 coefficients as arguments; but, if 4 arguments are
present, then a*x*^3 + b*x*^2 + c*x* + d = 0 is
assumed.

Even though both quadratic() and cubic() functions take real arguments, they can be modified to accept any real or complex coefficients because the method of solution does not make any assumptions. [4]

**Summary**- Solve for a zero of function using Newton-Raphson method

**Usage**-
`real`= func(`real`)

`real`= funcd(`real`)

`real`= newton(func, funcd,`real`[, TOL=`real`])

Even though it converges quadratically once a root has
been "sighted", it does not guarantee global convergence. So, I use
**print** statement to see intermediate results.

**Summary**- Solve for a zero of function using Secant method

**Usage**-
`real`= func(`real`)

`real`= secant(func,`real`,`real`[, TOL=`real`])

This routine is used when f'() is difficult or impossible to evaluate. It requires 2 initial points to start the iteration, and the order of convergence is 1.6, which is better than linear method like bisection but not as fast as quadratic method like Newton.

**Summary**- Closed integration of function using extended Simpson's rule

**Usage**-
`real`= func(`real`)

`real`= closedpoints(func,`real`,`real`[, TOL=`real`])

**Summary**- Open integration of function using extended Simpson's rule

**Usage**-
`real`= func(`real`)

`real`= openpoints(func,`real`,`real`[, TOL=`real`])

Both closedpoints() and openpoints() routines use Trapezoid formula to calculate Simpson's formula, as opposed to more direct summation. And, because the discrete integration is over uniformly sampled points, it can take advantage of already sampled data taken at previous iterations.

**Summary**- Find 2^n that is equal to or greater than

**Usage**`int`= nextpow2(`int`)

This is internal function used by fft(), because the FFT routine requires that the data size be a power of 2.

**Summary**- Return bit-reversed array for FFT

**Usage**`list`= bitrev(`list`)

**Summary**- FFT and inverse FFT for 2^n data size [5]

**Usage**-
`list`= fft(`list`[, sign=1])

`list`= ifft(`list`)

I calculate N points along the unit circle at the
beginning instead of calculating as needed; so, this should be more
efficient than the FFT routine given in *Numerical Recipes in
C*, but not by much.

**Summary**- DFT and inverse DFT using direct summation

**Usage**-
`list`= dft(`list`[, sign=1])

`list`= idft(`list`)

I wrote these functions to test FFT routines. DFT should only be used for small data size (less than 64).

**Summary**- Discrete convolution and correlation [5]

**Usage**-
`list`= conv(`list`,`list`)

`list`= corr(`list`,`list`)

**Summary**- FFT convolution and correlation [5]

**Usage**-
`list`= fftconv(`list`,`list`)

`list`= fftcorr(`list`,`list`)

Due to the nature of FFT routines, fftconv() and fftcorr() assume that input data have been zero-padded to a length that is big enough to hold the output size P+Q-1. However, conv() and corr() can accept any length, because they use direction summation.

**Summary**- Vector dot

**Usage**`real`= vdot(`list`,`list`)

**Summary**- Vector norms

**Usage**-
`real`= vnorm(`list`, normtype=1)

`real`= vnorm(`list`)

`real`= vnorm(`list`, normtype=3)

**Summary**- Mean and standard deviation of data

**Usage**`real`,`real`= meanstdv(`list`)

**Summary**- Read vector from stdin or file

**Usage**-
`list`= getv([`string`])

`list`,`list`= getv2([`string`])

**Summary**- Print vector to stdout or file

**Usage**- printv([
`string`])

**Summary**- Linear regression of y = ax + b

**Usage**`real`,`real`= linreg(`list`,`list`)

Only the coefficients of regression line are returned, since they are usually what I want. Other informations are sent to stdout to be read later.

- [1]
- W H Press, B P Flannery, S A Teukolsky, W T
Vetterling,
*Numerical Recipes in C: The Art Of Scientific Computing*, 1988. http://www.nr.com/ - [2]
- NumPy is numerical array extension to Python. ftp://ftp-icf.llnl.gov/pub/python/
- [3]
- Octave is GNU clone of Matlab. http://www.che.wisc.edu/octave/
- [4]
- W Gellert, H Kustner, M Hellwich, H Kastner,
*The VNR Concise Encyclopedia of Mathematics*, 1975, sections 4.5-4.6 - [5]
- E O Brigham,
*The Fast Fourier Transform*, 1974

- 1 Introduction
- 2 Elementary Funtions
- 3 Polynomials
- 4 Zeros of Functions
- 5 Integration
- 6 FFT
- 7 Convolution and Correlation
- 8 Vector Operations
- 9 Footnotes
- 10 Table of Contents

This file was generated using

William Park opengeometry@yahoo.ca

March 1999