You should be aware that Python is slow compared to native code like Java (which is compiled at runtime). The factor is always more than 10 and can be more than 100. So be warned!

Often, it is possible to improve the runtime by switching the algorithm. Let us study the Collatz problem an example. There is no need to explain this problem here. You can easily find a description on Wikipedia. For a start, we want to find the longest sequence for any starting point less than a million, and later for hundred millions.

The following is a recursive function to get the length of the Collatz sequence starting at n.

def get_collatz_length (n):
    if n==1:
        return 0
    elif n%2:
        return get_collatz_length(3*n+1)+1
    return get_collatz_length(n//2)+1

To compute the lengths up to N=1’000’000 takes around 20 seconds. We added a code to compute the 10 worst starting points.

import time
start = time.process_time()
N=1_000_000
lengths = [(n,get_collatz_length(n)) for n in range(1,N)]
print(f"{time.process_time()-start:5.4} seconds")

records = sorted(lengths,key=lambda x:x[1],reverse=True)
for k in range(0,11):
    print(f"{records[k][0]}: {records[k][1]}")

18.72 seconds
837799: 524
626331: 508
939497: 506
704623: 503
910107: 475
927003: 475
511935: 469
767903: 467
796095: 467
970599: 457
546681: 451

A comparable Java code is the following.

public class Main 
{
	
	public static int get_collatz_length (long n)
	{
		int res = 0;
		if (n>1)
		{
			if (n%2 == 1)
				res = get_collatz_length(3*n+1)+1;
			else
				res = get_collatz_length(n/2)+1;
		}
		return res;
	}

	public static void main(String[] args) 
	{
		long time = System.currentTimeMillis();
		int N=1000000;
		int L[] = new int[N];
		for (int k=0; k<N; k++)
			L[k] = get_collatz_length(k);
		int max=0,imax=0;
		for (int k=1; k<N; k++)
			if (L[k]>max)
			{
				max=L[k]; imax=k;
			}
		System.out.println("Maximal value " + max + " at "+imax);
		time = System.currentTimeMillis()-time;
		System.out.println((time/1000.0) + " seconds.");
		
	}

}

Maximal value 524 at 837799
0.498 seconds.

This takes about half a second. The factor is 40.

What we also see is that the Python code looks more elegant and mighty. E.g., sorting the array of lengths and printing the 10 highest is done with three lines of code. In Java, the same can only be achieved by writing a Comparator class.

Consequently, I see Python as an interactive script language to do math. If speed is needed, native code must be used, as in Numpy. Unfortunately, the user cannot easily write a function in C or Java.

Can we speed up the process? We could try to avoid the doubled computations by remembering and using the already known lengths in a dictionary. The following code does this.

collatz_lengths = {}

def get_collatz_length (n, nmax):
    res = collatz_lengths.get(n)
    if res:
        return res 
    elif n==1:
        res = 0
    elif n%2:
        res = get_collatz_length(3*n+1,nmax)+1
    else: 
        res = get_collatz_length(n//2,nmax)+1
    if n<nmax:
        collatz_lengths[n] = res
    return res

The storage of the lengths is only done if n is not too high. For the call, we have set this maximal value equal to N=1’000’000. This code reduces the time to about 1 second. We gained a factor of 20 just by this simple trick. Obviously, an idea might be worth more than a compiled code.

Maybe, we can do oven better, if we do not use a dictionary. Instead, we can use a list and store the known lengths in it. Let us try.

def get_collatz_length (n, L, N):
    if n<N and L[n]>=0:
        return L[n]
    if n==1:
        res = 0
    elif n%2:
        res = get_collatz_length(3*n+1,L,N)+1
    else: 
        res = get_collatz_length(n//2,L,N)+1
    if n<N:
        L[n]=res
    return res

import time
start = time.process_time()
N=1_000_000
L = list(range(N))
for k in range(N):
    L[k] = -1
lengths = [(n,get_collatz_length(n,L,N)) for n in range(1,N)]
records = sorted(lengths,key=lambda x:x[1],reverse=True)
for k in range(11):
    print(f"{records[k][0]}: {records[k][1]}")
print(f"{time.process_time()-start:5.4} seconds")

This reduces the time only slightly. It is still one second.

However, implementing the same trick in Java, we get the lengths up to one million in 0.033 seconds, i.e., 30 times faster. The reason for this is partly the infinite integer type of Python, and also the use of a dictionary instead of an array. Here is the Java code.

public class Main 
{
	
	public static int get_collatz_length (long n, int L[], int N)
	{
		if (n<N && L[(int)n]>=0)
			return L[(int)n];
		int res = 0;
		if (n>1)
		{
			if (n%2 == 1)
				res = get_collatz_length(3*n+1,L,N)+1;
			else
				res = get_collatz_length(n/2,L,N)+1;
		}
		if (n<N) L[(int)n] = res;
		return res;
	}

	public static void main(String[] args) 
	{
		long time = System.currentTimeMillis();
		int N=100000000;
		int L[] = new int[N];
		for (int k=0; k<N; k++)
			L[k] = -1;
		for (int k=0; k<N; k++)
			L[k] = get_collatz_length(k,L,N);
		int max=0,imax=0;
		for (int k=1; k<N; k++)
			if (L[k]>max)
			{
				max=L[k]; imax=k;
			}
		System.out.println("Maximal value " + max + " at "+imax);
		time = System.currentTimeMillis()-time;
		System.out.println((time/1000.0) + " seconds.");
		
	}

}

Increasing N to 10’000’000, the Java computation time rises to 0.15 seconds. The more clever Python code was faster now and took 9 seconds, i.e., it was 60 times slower.

Going further to N=100’000’000 made Java run for 1.44 seconds. The more clever code in Python, however, broke down and took more than 2 minutes, again a factor of over 50.

Finally, I wanted a plot of the frequencies of the lengths. I had to write a procedure for the count, not trusting the performance of the Counter class.

lengths1 = [x[1] for x in lengths]
M = max(lengths1)
def get_count (L,M):
    count = [0 for k in range(M)]
    for x in L:
        if x<M:
            count[x] += 1
    return count
count = get_count(lengths1,M+1)

This allows to plot the frequencies. The vertical axis is the number of times a length occurs among the first 100 million numbers. The same plot can be found on Wikipedia.

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

fig,ax = plt.subplots(figsize=(14,6))
ax.grid(True)
ax.plot(count);

In summary, Python did well. But it seems to be so flexible that it is not predictable in terms of performance and memory consumption. That is so with many modern software systems who develop a mind of their own from the perspective of a user.

Continuing my efforts to learn Python for Math, I wrote a function for adaptive integration. It uses 10-point Gauss quadrature. To check the accuracy, the algorithm uses a bipartition of the interval. So each step takes 30 evaluations of the function, one for the complete step and two for the check. If the accuracy is not high enough, the step size will decreased. Otherwise, the function will try to double the step size in the next step.

Here is the complete code.

def legendre (n):
    """ Get the Legendre Polynomials
    
    n : maximal degree
    
    Return: list of Lenegdre Polynomals
    """
    p0 = np.poly1d([1])
    if n<=0:
        return [p0]
    p1 = np.poly1d([1,0])
    px = p1
    if n==1:
        return [p0,p1]
    v=[p0,p1]
    for i in range(1,n):
        p = ((2*i+1)*px*p1 - i*p0) / (i+1)
        p0 = p1
        p1 = p
        v.append(p1)
    return v

def make_gauss(n):
    """ Compute the Coefficiens for Gauss Integration
    
    n : Number of points
    
    return: x,a
    x : Gauss popints
    a : Gauss coefficients
    """
    v = legendre(n)
    x = np.roots(v[n])
    A = np.array([v[i](x) for i in range(n)])
    w = np.zeros(n)
    for i in range(n):
        p = np.polyint(v[i])
        w[i] = p(1)-p(-1)
    return x,np.linalg.solve(A,w)

gauss_x5,gauss_a5 = make_gauss(5)

def gauss5 (f,a,b):
    """ Gauss Integrate with 5 points
    
    f : function of one variable
    a,b : Inteterval bounds
    
    returns: Integral
    """
    return np.sum(f((a+b)/2+(b-a)/2*gauss_x5)*gauss_a5)*(b-a)/2

gauss_x10,gauss_a10 = make_gauss(10)

def gauss10 (f,a,b):
    """ Gauss Integrate with 10 points
    
    f : function of one variable
    a,b : Inteterval bounds
    
    returns: Integral
    """
    return np.sum(f((a+b)/2+(b-a)/2*gauss_x10)*gauss_a10)*(b-a)/2

def integrate (f,a,b,n=1,eps=1e-14):
    """ Adaptive Integral using 10 point Gauss
    
    f : function of one variable
    a,b : Interval bounds
    n=1 : Initial step size is (b-a)/n
    eps : Desired absolute accuracy
    
    returns: Integral
    """
    h = (b-a)/n
    x = a
    I = gauss10(f,x,x+h)
    I1 = gauss10(f,x,x+h/2)
    I2 = gauss10(f,x+h/2,x+h)
    res = 0.
    while x<b-eps:
        ## print(h)
        if np.abs(I-(I1+I2))>h/(b-a)*eps:
            h = h/2
            I = I1;
            I1 = gauss10(f,x,x+h/2)
            I2 = gauss10(f,x+h/2,x+h)
        else:
            x = x+h
            h=2*h
            if x+h>=b-eps:
                h=b-x
            res = res+I1+I2
            I = gauss10(f,x,x+h)
            I1 = gauss10(f,x,x+h/2)
            I2 = gauss10(f,x+h/2,x+h)
    return res

This algorithm is usually fast enough and very accurate. The implementation is by no means optimal. E.g., if a factor of 2 is used in the step size, a lot of computations are redone. Moreover, it might not be optimal to try to double the step size in each step. However, it works for interactive use, and that is what it is intended for.

Let me try to explain the code a bit. The function legendre() computes the Legendre polynomials using the recursion formula. These polynomials are orthogonal with respect to the weight w=1. The function returns a list of the polynomials. It makes use of the polynomial handling in Python (in old style).

The function make_gauss() computes the zeros of the Legendre polynomials using a function from Python, and coefficients a[k], such that the following is exact for all polynomials p up to degree n-1:

\(\sum_{k=1}^n a_k p(x_k) = \int\limits_{-1}^1 p(t) \,dt\)

It sets up n linear equations from this requirement using the Legendre polynomials we just computed as test cases and solves the equations. To compute the correct integral for polynomials, there is a function in Python.

Now, we can define the functions for Gauss integration with 5 and 10 points, and use it for our adaptive integration.

For the demonstration, I created a very easy to use function fplot() which plots a function adaptively using the code of the last posting.

def fplot (f,a,b,n=100,nmax=100000):
    """ Simple plot of a function
    
    a,b : interval bounds
    n=100 : initital step size (a-b)/n
    nmax=100000 : minimal step size (a-b)/nmax
    
    """
    x,y = adaptive_ev(f,a,b,n,nmax)
    fig,ax = plt.subplots(figsize=(8,8))
    ax.grid(True)
    ax.plot(x,y)

I now have put all this and the code of the last posting into a Python module. When using the Python notebook interface in a Browser, the file should be placed in the same directory where the „.ipynb“ file. The Python file can even be opened and edited in the Python notebook IDE.

Then it can be used as follows.

import renestools as rt
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

def f(x):
    return np.abs(x)

rt.fplot(f,-np.sqrt(2),1,n=1000)

print(rt.integrate(f,-np.sqrt(2),1))
1.5

This is not easy to integrate numerically, since the function is not smooth. The result is exact to the last digit. Even the plot is not trivial. With the arbitrary step size, the point x=0 might not be met exactly, and the corner may look round.

Here is another example using the Gaussian distribution. All digits agree. Note, that the Gauss function is almost zero for large x.

def f(x):
    return np.exp(-x**2)

print(rt.integrate(f,-20,20))
print(np.sqrt(np.pi))

1.7724538509055159
1.7724538509055159

You can even make the integral to a function and plot it adaptively.

def F(x):
    return rt.integrate(f,0,x)/np.sqrt(np.pi)

rt.fplot(F,0,5)

I could not find a solution for adaptive plots in Matplotlib. So, I wrote a small script that can evaluate a function in adaptive step sizes. „Adaptive“ means that the step size adapts to the output of the plot. This is to prevent coarse line segments with corners in the plot.

Here is the simple code for a function of one variable.

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

def adaptive_ev (f,a,b,n=100):
    xc = a
    yc = f(xc)
    x = [xc]
    y = [yc]
    h = (b-a)/n
    uplim = (2*h)**2
    lowlim = uplim/6
    eps = (b-a)/100000.
    while xc<b:
        xcnew=xc+h
        if xcnew>=b:
            xcnew=b
        ycnew=f(xcnew)
        d = h**2+(ycnew-yc)**2 
        if d>uplim and h>eps:
            h = h/2
        else:
            if d<lowlim/3:
                h = 2*h
            xc = xcnew
            yc = ycnew
            x.append(xc)
            y.append(yc)
    return np.array(x),np.array(y)

def func(x):
    return np.sin(1/x)

x,y=adaptive_ev(func,0.001,1)
print(x.size,' points generated')
fig,ax = plt.subplots(figsize=(14,10))
ax.grid(True)
ax.plot(x,y);

To apply the function adaptive_ev(), you need a function of one variable, an interval and the target number of plot points. The function will double or half the step size as needed to keep the output under control.

This is a screenshot of the graph generated by the code above. The function is difficult to plot near 0. Thus, adaptive evaluate is the only chance to do that. The routine generated 4270 data points to create this plot.

Euler Math Toolbox has the adaptive evaluation built into the plot2d() routine. So, a similar code would look like this.

>aspect(1.2); plot2d("sin(1/x)",0.001,1,thickness=1.5,color=blue):

Admittedly, the anti-aliasing of Matplotlib is outstanding. I like this library a lot and recommend using it.

I also did a routine for curves in the plane, which are parameterized by two functions for the x- and y-coordinate.

def adaptive_ev2 (fx,fy,a,b,size=1,n=1000):
    tc = a
    xc = fx(tc)
    yc = fy(tc)
    x = [xc]
    y = [yc]
    h = size/n
    uplim = (2*h)**2
    lowlim = uplim/6
    eps = size/100000.
    while tc<b:
        tcnew=tc+h
        if tcnew>=b:
            tcnew=b
        xcnew=fx(tcnew)
        ycnew=fy(tcnew)
        d = (xcnew-xc)**2+(ycnew-yc)**2 
        if d>uplim and h>eps:
            h = h/2
        else:
            if d<lowlim/3:
                h = 2*h
            xc = xcnew
            yc = ycnew
            tc = tcnew
            x.append(xc)
            y.append(yc)
        ## print((tc,xc,yc,h,d))
    return np.array(x),np.array(y)

def adaptive_ev2 (fx,fy,a,b,size=1,n=1000):
    tc = a
    xc = fx(tc)
    yc = fy(tc)
    x = [xc]
    y = [yc]
    h = size/n
    uplim = (2*h)**2
    lowlim = uplim/6
    eps = size/100000.
    while tc<b:
        tcnew=tc+h
        if tcnew>=b:
            tcnew=b
        xcnew=fx(tcnew)
        ycnew=fy(tcnew)
        d = (xcnew-xc)**2+(ycnew-yc)**2 
        if d>uplim and h>eps:
            h = h/2
        else:
            if d<lowlim/3:
                h = 2*h
            xc = xcnew
            yc = ycnew
            tc = tcnew
            x.append(xc)
            y.append(yc)
        ## print((tc,xc,yc,h,d))
    return np.array(x),np.array(y)

def funcx(t):
    return np.exp(-t/10)*np.cos(t)
def funcy(t):
    return np.exp(-t/10)*np.sin(t)

x,y=adaptive_ev2(funcx,funcy,0,40*np.pi)
print(x.size,' points generated')
fig,ax = plt.subplots(figsize=(10,10))
ax.grid(True)
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
ax.plot(x,y);

The code follows the same ideas. However, it uses a plot size to set the dimensions of the curve in the plane. You only need a guideline. The step size will be approximately 1/1000 of this size.

Close to 0, the plot should avoid generating too many points. But it should look smooth. This plot generates 14509 points. That is 14 times the given size of 1. But one should be aware that it winds 20 times around the origin!

Euler Math Toolbox does the same, if two functions are fed into plot2d(). It uses a function adaptive() to compute the data points, which works for vector valued functions. This function could be ported to Python easily.

function adaptive (f$: call, a: number, b: number, eps,  ..
dmin, dmax, dinit, factor)

## Default for eps : 0.01
## Default for dmin : 1e-05
## Default for dmax : 100000
## Default for dinit : 0.1
## Default for factor : 2

    sn=sign(b-a);
    t=a;
    sy=f$(t;args());
    if sn==0 then return {a,sy}; endif;
    if cols(sy)>1 then 
        error(f$|" does not produce a real or column vector"); 
    endif
    x=zeros(1,1000); y=zeros(rows(sy),1000); 
    n=1; h=dinit;
    x[:,1]=t;
    y[:,1]=sy;
    repeat
        repeat
            tn=t+sn*h;
            if sn*(tn-b)>=0 then tn=b; endif;
            syn=f$(tn;args());
            d=sqrt(colsum(abs(syn-sy)^2));
            if isNAN(d) then break; endif;
            if d>eps and h>dmin then h=max(h/factor,dmin);
            elseif d<eps/factor then h=min(h*factor,dmax); break;
            else break;
            endif;
        end;
        n=n+1;
        if n>cols(x) then
            x=x|x;
            y=y|y;
        endif;
        x[:,n]=tn; y[:,n]=syn;
        sy=syn;
        t=tn;
        if sn*(t-b)>=0 then break; endif;
    end;
    return {x[:,1:n],y[:,1:n]};
endfunction

I find Python very aspiring for mathematics. The next step will be to write an adaptive integration using Gauss formulas. Stand by.