Programmer-defined functions

We have already seen that an important part of Python programming is calling functions to perform calculations or carry out other useful tasks for us. Examples of functions we have used include print() and input(). Functions are so useful for programming that it is only natural that Python provide a way for us to create our own functions.

The basic form for a function definition in Python is

def <name>(<parameter list>):
  <optional doc string>
  <statements>
  <return statement>

Here is a simple example that illustrates all of the key features of a user-defined function.

def square(x):
  """Compute the square of the argument"""
  result = x*x
  return result

This example function computes and returns the square of an integer. Let's take a closer look at the components of this function.

The first line of the function definition sets up the interface to the function, including its name and the parameters it takes.

def square(x):

The interface declaration tells us that the name of this function is 'square' and that it takes a single parameter named 'x'. Functions can take any number of parameters. To set up multiple parameters we simply make a comma-separated list of parameters.

The body of the function is the indented portion. The body contains the statements that carry out the computation and ultimately return a result.

The first line of the body can be an optional doc string, which is some text that documents what the function does. Many development environments, like PyCharm, look for these doc strings and can show you that documentation when you hover over the name of a function. We always use the special Python triple quoted string format for doc strings, since this format allows us to construct strings that span more than one line.

def square(x):
  """square(x) computes the square of its argument, x.
     There is really nothing more to it than that.
     Please don't try to pass something other than a number
     to this function when you call it. The results will
     be unpredictable."""
  result = x*x
  return result

Next we will have one or more statements that carry out some computation. The inputs that come in to the function via the parameters play a role in these computations.

result = x*x

The final statement is the return statement, which is responsible for returning the result of the computation.

return result

Additional rules for functions

Besides the structure rules I described above for defining functions there are some additional rules you have to follow to use functions in a Python program.

  1. Functions have to be defined before they are used. If code in one function wants to use another function that other function definition has to appear before the definition of the using function.
  2. You can define functions inside the body of other functions, but those nested functions can only be used inside the other function they were defined in.
  3. When you call a function you have to make sure that you are passing the number of parameters required by that function.
  4. A function can return only one value, although you can have multiple return statements in a function. As soon you hit one of the return statements you will exit the body of the function.
  5. Function names have to be unique. You can not define two functions with the same name.

Newton's Method

Newton's method is a method for finding roots of functions. Given a differentiable function f(x) we would like find a root of f, that is, a value of x such that f(x) = 0.

The method consists of the following steps:

  1. Pick a point x0 close to a root. Find the corresponding point (x0 , f(x0)) on the curve.
  2. Draw the tangent line to the curve at that point, and see where it crosses the x-axis.
  3. The crossing point, x1, is your next guess. Repeat the process starting from that point until you are satisfied that you are close enough to the root.

Mathematical details

We pick a point

(x0, f(x0))

on the curve and find the equation of the tangent line at that point. The slope of the tangent line at that point is

The tangent line has general form

This line intersects the x-axis when y = 0.

To repeat the process we just have to feed the point x1 back into the last formula to generate the next approximation for the root:

Each time we generate another approximation we can easily test to see how close that approximation is to the root by simply plugging that back into the function f. Once |f(xn)| drops below a preset tolerance we can stop the iteration.

Program

Here is the code for a Python program that uses Newton's method to find a root for the function

f(x) = x4 + 2 x3 - 14 x2 + 2 x + 1

This function has a positive real root near about 0.4. To find this root we will write a function that uses the Newton iteration with the function f(x) and its derivative

To make computing the polynomials f(x) and f(x) more efficient we will use a special trick called Horner's rule to evaluate the polynomials. Can you see how Horner's rule works?

# Program to find the roots of f(x) = x^4+2 x^3-14 x^2+2 x+1
# by Newton's method

# Define the function and its derivative
def f(x):
    return (((x+2)*x-14)*x+2)*x + 1

def fp(x):
    return ((4*x+6)*x-28)*x+2

x = input("Enter a starting guess for the root: ")
x = float(x)

y = f(x)
tolerance = 10e-6

while y < -tolerance or y > tolerance:
    # Compute f'(x)
    yp = fp(x)
    # Do one round of Newton's method
    x = x - y/yp
    y = f(x)

print('The root estimate is ',x)

Importing functions from modules

Most programming systems include a library of useful functions, and Python is no exception. For example, to make use of the standard mathematical functions such as exp(x) or sin(x), we place an import statement at the top of our program to import the math module.

import math

We can use any of the functions in that module, such as exp(), sin(), or sqrt(), by prefacing their names with the name of the module we imported them from.

print('exp(3.3) = ',math.exp(3.3))

Alternatively, you can use a variant of the import statement that imports only specific functions from a module. You can use those function names without prefacing them with the module name if you use this approach.

from math import sin, cos

print('sin(1.43) = ',sin(1.43))

Programming Exercise

Newton's method has a wide range of applications. In this exercise we will explore one of those applications.

Many mathematical functions can be computed based on other functions. An example of this is the function f(x) = ln x. One way to define the natural log function is to express the natural log as the solution of an equation. Consider the equation

x = et

Mathematically, the obvious solution to this equation is t = ln x.

Construct a program that uses this equation combined with Newton's method to compute ln x. Here is what you program should do.

  1. Prompt the user to enter a value for x.
  2. Use the Newton's method code above to construct some code that computes the root of the equation x = et.
  3. Print the value of t that Newton's method converged to along with the value of math.log(x) to see how close to the actual value you got.

Constructing a useful function: computing exp(x)

As an exercise in function writing, let us see if we can construct our own function to compute exp(x). The key to doing this is the power series expression for exp(x):

Here is a first draft of a function that can compute this power series.

def exp1(x):
    sum = 1
    for n in range(1,16):
        sum += x**n/math.factorial(n)
    return sum

Since the sum in the power series is an infinite sum, we will have to settle for a finite approximation. I determined the cutoff of n = 15 by trial and error, since this produced the best possible results with the fewest terms in the sum. Here is some output from a simple test program that confirms that this works well enough.

math.exp(0.3)= 1.3498588075760032
    exp1(0.3)= 1.349858807576003
math.exp(0.5)= 2.0137527074704766
    exp1(0.5)= 2.0137527074704766

Next, we introduce a simple optimization. Since x**n and math.factorial(n) are both expensive to compute, we take advantage of a simple relationship between successive terms in the sum.

This means that we can easily compute each new term in the sum from the term that came before it.

Here is a version of the exp function that uses this optimization.

def exp2(x):
    """A more refined version of the exponential function."""
    sum = 1 + x
    term = x
    for n in range(2,16):
        term = term*x/n
        sum += term
    return sum

Here is a quick test run to confirm that this of course produces the same numerical results as the earlier, less efficient version.

math.exp(0.3)= 1.3498588075760032
    exp1(0.3)= 1.349858807576003
    exp2(0.3)= 1.349858807576003
math.exp(0.5)= 2.0137527074704766
    exp1(0.5)= 2.0137527074704766
    exp2(0.5)= 2.0137527074704766

Our exponential computing function seems to be working OK, but it still has one last, serious weakness. The problem with power series approximations is that they don't converge equally well for all values of x. Here is another test run that shows that when we plug in a larger x the approximation gets much worse.

math.exp(3.3)= 27.112638920657883
    exp2(3.3)= 27.11262722548732

We could work around this problem by increasing the number of terms in the sum as the input x gets larger. Instead, I am going to pursue an alternative strategy. Since exp2 seems to do a pretty good job with inputs x in the range from 0 to 1, the trick is to map any x we are given into that range. We can do that by the following strategy:

  1. Compute the smallest power n such that x/2n < 1.
  2. Compute exp(x/2n).
  3. Note that
  4. exp(x) = (exp(x/2n))2n

  5. We can compute the latter expression by repeated squaring:
  6. (exp(x/2n))2n = ((exp(x/2n))2)2n-1

Here now is the final version of the exponential computing function that makes use of this strategy.

def exp3(x):
    """Final version of exp - works for a broader range of x."""
    # Initially we work with the absolute value of x
    z = math.fabs(x)
    # Compute the smallest n such that z/2**n <= 1
    n = 1
    while z > 1:
        n += 1
        z /= 2.0
    # Since z is now in the range 0 <= z <= 1, exp2 will
    # do a good job with it
    result = exp2(z)
    # Since now result = exp(z/2**n), to compute exp(z)
    # we now need to compute result**(2**n). We do this
    # by repeated squaring of result, since
    # (result**2)**(2**(n-1)) = result**(2**n)
    while n > 1:
        result = result*result
        n -= 1
    # Finally, we handle the case of a negative x.
    if(x < 0):
        return 1/result
    else:
        return result

and a test run to confirm that this works much better than exp2:

math.exp(3.3)= 27.112638920657883
    exp2(3.3)= 27.11262722548732
    exp3(3.3)= 27.112638920657805

Creating and using a module

As programs get larger and more complex we will want to start organizing useful functions into modules that we can import into our programs. In the next example I am going to write a set of useful functions for working with prime numbers and then import those functions into a simple test program.

First, here is the module with the useful function definitions:

def isPrime(n,known_primes):
    """
    Test n to determine whether or not it is a prime number.
    The second parameter must be a list of prime numbers that
    includes all of the prime numbers smaller than the square
    root of n.
    """
    for p in known_primes:
        if p*p > n:
            return True
        if n % p == 0:
            return False
    return True

def primeFactors(n,primes):
    """Compute and return a list of the prime factors of n."""
    factors = []
    x = n
    for p in primes:
        while x % p == 0:
            factors.append(p)
            x = x//p

    if x > 1:
        factors.append(x)

    return factors

def findPrimes(m,n,primes):
    """
    Construct and return a list of all the prime numbers less
    than n. primes must be a seed list of prime numbers less
    than or equal to m.
    """
    myPrimes = primes.copy()
    for num in range(m,n):
        if isPrime(num,myPrimes) and not num in myPrimes:
            myPrimes.append(num)
    return myPrimes

Here is the test program that uses the functions in this module.

from prime_stuff import *

p = [2,3,5,7]
allPrimes = findPrimes(10,100,p)

x = int(input('Enter an integer between 2 and 10000: '))
if isPrime(x,allPrimes):
    print(str(x) + ' is a prime number.')
else:
    factors = primeFactors(x,allPrimes)
    print(str(x) + ' is not a prime number.')
    print("Its factors are " + str(factors))