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
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.
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:
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.
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)
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))
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.
math.log(x)
to see how close to the actual value you got.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:
exp(x) = (exp(x/2n))2n
(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
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))