Teaching Computer Programming to High School students: An introductory course using Python as the high level language

The number system we're used to is called decimal. The base of the decimal system is 10. Number systems used in computing are

bit, wikipedia (http://en.wikipedia.org/wiki/Bit).

Computers crunch numbers. We (people) use the decimal number system. Computers only have bits and use these bits to make binary numbers. There are no decimal numbers inside a computer. Binary numbers are transformed by software to a decimal representation on the screen for us.

Before launching into binary numbers, lets refresh our memories on the positional nature of the representation of decimal numbers.

102=1*100  + 0*10   + 2*1
   =1*10^2 + 0*10^1 + 2*10^0

The value represented by each digit in the number 102 depends on it's position. The "1" represents "100". If the "1" was in the rightmost position it would represent "1". Each time a digit moves 1 place to the left, it increases the number it represents by a factor of the base of the number system, here 10. The left most digits represent the biggest part of the number.

binary numeral system, wikipedia (http://en.wikipedia.org/wiki/Binary_numeral_system).

In a binary number, the base is two, so the number prepresented by each position increases by a factor of 2 as you move left in the number. In binary, you need two numbers to represent all the available values. Here are the two numbers and their decimal equivalents.

binary decimal
   0  =  0
   1  =  1

Here's a binary number. What is the decimal representation of this number?

1101 = 1*2^3 + 1*2^2 + 0*2^1 + 1*2^1
     = 1*8   + 1*4   + 0*2   + 1*1
     = 13 decimal

As with decimal, the left most digit carries the most value (in 1101₂, the left most digit represents 8₁₀).

leading zeroes behave the same as for decimal numbers

00=0
01=1

Here's some more binary numbers with their decimal equivalent

  binary decimal
      10=2
      11=3  (2+1)
     100=4
     101=5  (4+1)
    0101=5
00000101=5
    1001=9  (8+1)
    1011=11 (9+2,8+3))
   11011=27 (16+11, 24+3)
   11111=31

what is 1011₂ in decimal ^[2] ?

Let's go the other way, decimal to binary

what is 7₁₀ in binary? Do it by stepwise division.

Start with the power of 2 just below your number.
Take a guess. 2^3=8. This is bigger than 7.
Try next one down. 2^2=4. This is less than 7. Start there

7/4=1, remainder=3

do it again, the power of 2 just below 3 is 2

3/2=1, remainder=1

with the remainder being 0 or 1, we're finished

7=1*2^2 + 1*2^1 + 1*2^0

7 decimal = 111 binary

In the above exercise, would anything have gone wrong if you'd started dividing by 8 rather than dividing by 4? No, you would have got the answer 7₁₀=0111₂. The leading zero has no effect on the number, so you would still have got a right answer, you just would have done an extra step.

what is 15₁₀ in binary? ^[3]

Here's python code to convert decimal to binary (http://www.daniweb.com/code/snippet285.html).

Here's bash code to convert binary to decimal. If you aren't already in a bash shell, start one up in your terminal with the command /bin/bash. On windows, click on the Cygwin icon to get a bash prompt. (comments in bash start with #). The code here is from Bash scripting Tutorial (http://www.linuxconfig.org/Bash_scripting_Tutorial) in the section on arithmetic operations. (Knowing that bash can convert binary to decimal, I found this code with google using the search terms "convert binary decimal bash").

declare -i result	#declare result to be an integer
result=2#1111		#give result the decimal value of 1111 to base 2
echo $result		#output the decimal value to the screen
15
# same code in one line
declare -i result;result=2#1111; echo $result
15

In the previous section, I showed bash code I found on the internet. When writing a new program, you'll often need some functionality that's already been coded up (after 50yrs of computing, there's a lot of code available) that's available in books or on the internet. Books on any computer language will have lots of worked examples. It's often faster to find debugged and working code, than it is to write it yourself from scratch. You're expected to do this and everyone does it. For a class exercise, you may have to nut it out yourself, but when it comes to writing working code, you borrow when you can. When you use someone else's code, you should document where you got it.

byte, wikipedia (http://en.wikipedia.org/wiki/Byte).

It turns out to be convenient (hardware wise) to manipulate bits 8 at a time. 8 bits = a byte(B). Most PCs (in 2008) are 32 bit machines, meaning that the CPU manipulates 4 bytes (32 bits) at a time. Since these machines are running at somewhere between 100MHz and 2GHz, then they are doing between 400 and 8000 million byte operations/sec.

some 1Byte numbers expressed in decimal

byte     decimal
00001000=8
00001111=15
00010000=16
00100000=32
01000000=64
10000001=129
11111111=255

1 Byte can represent integers from 0-255

Bytes are always 8 bits. However data is shifted around according to the bus width of a computer. A 32 bit computer has 32 bit registers and 32 lines for addressing and fetching data. It can transfer data and instructions 4 bytes at a time. A term from the HPC (high performance computing or supercomputers) world, where 64 bit computing has been the standard for 30 yrs, the term word describes the width/size of a piece of data/instruction. In the HPC world, there are words of all sorts of lengths, including 128-bit. A 32-bit computer has a 32-bit (4 byte) word size.

Let's do some binary addition: The rules are similar to decimal. Here's the addition rules.

	0+0=0
	1+0=1
	1+1=0 with carry, 1+1=10

here's the addition rules in table form.

addition                carry
+ | 0 1              + | 0 1
-------              -------
0 | 0 1              0 | 0 0
1 | 1 0              1 | 0 1

	Note
	1+1=0 with a carry. 1+1!=10 (!= not equal). The extra leftmost digit, the "1" (as in "1+1=10") becomes the carry digit. It's handled separately by the computer. If you want it, you have to go find it. The equivalent table in decimal for a one digit computer would show that 8+7=5 (and not 15)

worked example, working from right to left, one bit at a time

 111		(what demical numbers are being added here?)
 010 +
 ---
  01
 1   carry

 001
1    carry

1001

What is 1001+0011₂ ^[4] ? What decimal numbers are represented?

	Note
	End Lesson 1

Algorithm, wikipedia (http://en.wikipedia.org/wiki/Algorithm).

Algorithm: A predetermined sequence of instructions (events) that will lead to a result.

e.g.the algorithm of right-to-left addition will lead to the sum of two numbers.

e.g.putting your books into your school bag in the morning will lead to them all being at school with you the next morning.

What country do we get the word "algorithm"? What does "al" mean? What other words start with "al" used in the same way? ^[5]

Some algorithms are better (faster, more efficient, scale better) than others. While computers are fast enough that even a bad algorithm can process a small number of calculations, it takes a good algorithm to process large numbers of calculations. Computers are expensive and it's only worth using a computer if you're doing large numbers of calculations. We're always interested in how an algorithm scales with problem size. So we're always interested in the goodness of the algorithm being used.

The measure of goodness of an algorithm is the rate at which the time for it's execution increases as the problem size increases (i.e. how the algorithm scales). If you double the size of the problem, does the execution time not change (independant of problem size), double (linear in problem size) or quadruple (scales as the square of the problem size)?

This change in execution time with problem size is called the Order of the algorithm.

Speed (Order) of the addition operation:

addition one bit at a time takes 8 steps for addition of 2 bytes. The Order of addition for addition, one step at a time, is O(n) i.e.time taken is proportional to n=number of bits.

There's a faster way: add the bits in parallel and handle the carries in a 2nd step

 111
 010
 ---
 101 add
 10  carry
 ---
1001

Addition done in parallel takes 2 steps, no matter how long the numbers being added (you need to increase the amount of hardware, but you only have to pay for the hardware once and the extra cost can be amortised over many addition operations).

The Order of parallel addition is O(1) i.e.is proportional to 1 i.e.addition takes the same time independant of the number of bits being added. The parallel addition algorithm is much faster than the stepwise mechanism.

Note

amortise: If one process is faster/better than another, but the faster process is more expensive but only has to be paid for once, and you can use the process as many times as you want, then you are are said to be amortising the extra cost over a the life of the machine. e.g. cost of stepwize adder=$10. cost of parallel adder=$20. If you're going to be doing 1015 additions before you retire the machine, the amortised extra cost of the parallel adder is 10-14$/addition. Most people will accept this extra cost, because of the increased speed will save them money for each run of the program. The most common place that amortisation is used in computing is in the cost of writing a program. Writing a program is expensive; you have to pay the programmers salary, benefits, heating/AC, electricity and buy computers and lease a premises to do this. A program may cost $1k-$100M to write. However if the program is run millions of times, the cost/run for the end user may be insignificant compared to the cost of paying their staff to run the program. In this case the costs of the programmer's time is said to be amortised over the number of times the program will be run. Because writing programs is so expensive, you only write programs that are going to be run many times.

Here's the time/steps for the two types of addition for different problem sizes

	Stepwise	Parallel
1	2		2
2	4		2
3	6		2
4	8		2
.
.
16	32		2

We don't really care what the constant of proportion is, i.e.we don't care if each step takes 1usec, 1msec or 1sec, only how the time to completion scales with problem size. We know that if we scale the problem by some large number (e.g.10⁶), that the constant of proportionality will be swamped by the problem size. Let's say that parallel addition took 8 steps instead of 2. Here's the new times for addition.

bits	Stepwise	Parallel
1	2		8
2	4		8
3	6		8
4	8		8
.
.
16	32		8

We only need to get to numbers of length 8 bits to be ahead with parallel addition.

say we have a 4bit computer, what is

 1010
 0110+
 ----

^[6] ?

This addition results in a digit rolling off the left hand end. Loosing a digit off the left end of the register is called overflow (underflow is loosing a digit off the right end of the register). If you do this addition on a 4bit computer, you'll get an erroneous answer. In some circumstances you'll get an error message on overflow and in other situations you won't. Since overflow is part of the design of a computer, it is expected and is not neccessarily regarded as an error. (A drinking glass will only hold so much fluid. If you put in more fluid, it will overflow. People accept overflow as part of the design of a drinking glass.)

With addition, you have to anticipate getting a number that is 1 bit longer than the two numbers being added, i.e. you'll get a 1 bit overflow.

the rules are similar to decimal - here's the multiplication table (no carry is needed for binary multiplication)

* | 0 1
-------
0 | 0 0
1 | 0 1

5.9.1. left shifting (and overflow), multiplying by 10

Multiplication by 1,10,100 left shifts the digits by 0,1,2...n places (i.e. by adding a '0' on the righthand end). This is the same whether you're working in decimal or binary. Left shifting is a fast operation in a computer. The computer uses left shifting to multiply.

	Note
	both in binary and decimal 1*10= 10 11*10= 110 110*10=1100 what numbers are represented in each case?

Left shifting produces overflow. Assume a 4bit computer

1100* 10=1000 (not 11000) 1100*100=0000 (not 110000)

If a 0 overflows, you get the correct result.

 1010
   11x
 ----
 1010
1010 +
 ----
11110 addition
 0000 carry
-----
11110

    1010
    1101x
    ----
    1010
   0000
  1010
 1010   +
 -------
 1110010
   1     carry (a couple of rounds)
 -------
100000010

the carries are done in parallel
multiplication is fast and is O(1).

It is not neccessary to have a separate set of hardware for subtraction. You could in principle add the negative of the number, but we don't have negative numbers. Hardware only knows about 0-255 and has no concept of positive or negative numbers.

	Note
	Software can interpret certain values of a byte as being negative, but we haven't got that far yet.

However we don't need negative numbers for substraction. Instead we add the complement (for binary, it's the two's complement: the two's complement instruction is fast - one clock cycle). Let's find out about the complement.

5.10.1. Ten's complement

Let's look at subtraction by adding the complement using a decimal example.

Let's say we want to do

9-3=6

and we have a 1 decimal digit (dit) computer. Here's the addition table for a 1 dit computer.

+ | 0 1 2 3 4 5 6 7 8 9 ----------------------- 0 | 0 1 2 3 4 5 6 7 8 9 1 | 1 2 3 4 5 6 7 8 9 0 2 | 2 3 4 5 6 7 8 9 0 1 3 | 3 4 5 6 7 8 9 0 1 2 4 | 4 5 6 7 8 9 0 1 2 3 5 | 5 6 7 8 9 0 1 2 3 4 6 | 6 7 8 9 0 1 2 3 4 5 7 | 7 8 9 0 1 2 3 4 5 6 8 | 8 9 0 1 2 3 4 5 6 7 9 | 9 0 1 2 3 4 5 6 7 8

	Note
	remember, because we have a one digit machine, 9+9=8. (9+9!=18, the 1 is lost by overflow.)

If we're going to do the substraction 9-3=6, by adding some number to 9, what number do we add to get 6? Looking at the addition table, we find we have to add 7.

9-3=6 what we want 9+?=6 what we're looking for 9+7=6 the answer from looking up the addition table above

The ten's complement of 3 then is 7.

What if we want to subtract 3 from any other number, say 8? If we want to do 8-3=5, by adding a number to 8, on looking at the addition table, we have to add 7. So whenever we want to subtract 3, instead we add 7.

8-3=5 what we want 8+?=5 what we're looking for 8+7=5 the answer from lookup up the addition table above.

We find the ten's complement of 3 is 7 no matter what number we subtract 3 from. Making the complement of a number only depends on the number, not what we subtract it from.

What is the ten's complement of 8?

9-8=1 9+?=1 9+2=1

the ten's complement of 8 is 2.

What's the ten's complement of 9?

9-9=0 9+?=0 9+1=0

the ten's complement of 9 is 1.

What is the ten's complement of 0, of 4? ^[9]

Overflow isn't an advantage or a disadvantage; it's just part of the design of a computer. Since we have overflow, we can use it to do subtraction by addition of the complement, rather than having to build a subtractor into the hardware.

Here's the decimal (10's or ten's) complement table

number complement 0 0 1 9 2 8 3 7 4 6 5 5 6 4 7 3 8 2 9 1

The complement is the number you have to add to a number to get a sum of 0.

	Note
	the sum is not really 0; the sum is 10, but the left digit is lost through overflow.

The complement of a single digit number then is

complement=(base of the system, here 10)-(the number).

But in a 1 dit computer you don't have two digits to make a 10. Instead if you want the ten's complement of 7, you ask the computer to come in 3 places from the biggest number (here 9), i.e.8,7,6 giving 6 and then you add 1 giving 7.

	Note
	Subtracting from 10 or marching in from 9 are the same to you, but if you're wiring up a computer, you can't subtract from 10, but you can count in from 9.

Summary: if we're going to do 9-3, we add 10 to the -3, giving +7. The answer we'll get by doing 9+7 will be 10 more than what we want. However the extra 10 will be lost through overflow (subtracting 10 from the result), giving us the correct answer.

what we wanted 9-3=6 what the computer calculated (it added 10 to both sides) 9+(-3+10)=10+6 the answer the computer gave us 6

5.10.2. Two's complement

What is the two's complement of the 2 bit number 01₂?

find some number that when added to 01₂ gives 00₂ (just start poking in numbers till you get the required answer).

01 #decimal 1 11 #decimal 3 ---- 00 Note: overflow is needed to get 00

See any connection between 1 and its complement 3, in a 2 bit system ^[10] ?

A 4 bit number system has base 16. See any connection between the values of the number-complement pairs and a 4 bit number system in the following examples? Using brute force, what's the two's complement of the 4 bit numbers

• 0110 ^[11] ?
• 1011 ^[12] ?
• 1001 ^[13] ?
• 1100 ^[14] ?

let's try an example in a regular 8 bit byte. Using brute force, what's the two's complement of 01000110? (for labelling, let's use the words minuend, subtrahend and difference.)

01000101 subtrahend 10111011 complement -------- 00000000

How do you make the complement in binary?

Following the decimal examples (above), to get the complement, you count in from the end number (1 or 0) by 1 number (i.e. you flip the bits), shown here

01000101 original subtrahend 10111010 bit flipped subtrahend 10111011 known two's complement

	Note
	By looking at the bit flipped number and the two's complement, you can see that you have to add 1 (as is done for the decimal example).

10111011 bit flipped subtrahend + 1

binary complement=(bit flipped subtrahend + 1)

we've found the complement (the -ve + the base number of the system)

01000101 subtrahend (decimal 69) 10111011 complement, (decimal 187)

What's the sum of the 8 bit subtrahend and its complement ^[15]

with the complement, we can do the subtraction.

10000000 minuend (decimal 128) 01000101- subtrahend (decimal 69) -------- 10111010 bit flipped subtrahend 10111011 bit flipped subtrahend +1 = complement of 69

Do the subtraction by adding the two's complement

10000000 minuend (decimal 128) 10111011 complement of decimal 69 --------- 00111011 difference (left bit overflows on an 8 bit computer) result: binary decimal 10000000 128 01000101- 69- ------- --- 00111011 59

using the two's complement to do subtraction, what is

• 01101101 01100011- --------

  ^[16] ?

• 10110110 10001111- --------

  ^[17] ?

	Note
	End Lesson 2. Some kids didn't get the material on the complement and didn't complete the excercises. I added more exercises and started Lesson 3 at the beginning of binary subtraction.

Most people can only do 4bits of binary in their head. You either go to hexadecimal (below) or use a binary calculator. Fire up a terminal and try a few examples (you can recall the last line with the uparrow key and edit it without having to type in the whole line again).

bc (basic calculator?) is a general purpose calculator. Using a terminal, try some standard arithmetic e.g.(+-/*).

echo "3*4"  | bc
12

bc does all input and output in decimal, until you tell it otherwise.

Here's a few binary examples.

#input will be in binary, output is decimal since you haven't changed output
echo "ibase=2; 1100-101" | bc
7

#with obase explicitly set (not needed if obase is 10)
echo "obase=10;ibase=2; 1100-101" | bc
7

#same problem, output in binary
echo "obase=2;ibase=2; 1100-101" | bc
111

#convert decimal to binary
echo "obase=2; 17" | bc
10001

#other examples:
echo "obase=10;ibase=2; 1100+101" | bc
17

echo "obase=2;ibase=2; 1100+101" | bc
10001

Exercises: Hint - the number(s) you're processing are in the last instruction on the line. Before you run the instruction, figure out the base for the input and for the output and then decide whether you need to set obase and/or ibase.

The normal order is obase/ibase. What happens if you reverse the order of obase and ibase without changing their values?

normal order
echo "obase=10;ibase=2;01110001" |bc
113

inverted order
echo "ibase=2;obase=10;01110001" |bc
1110001

bc defaults to decimal input. In the normal order, bc interprets obase as 10₁₀ and ibase as 2₁₀ (i.e. binary). The input will be intepreted as binary and output will be in decimal. In the inverted order, obase says that all further input will be interpreted as base 2₁₀ (i.e. binary). Thus the obase value is 10₂ (2₁₀), i.e. the answer will be in binary.

As long as you know what you're doing, you can use obase,ibase in any order. To minimise suprises, use obase first, leaving the input decimal, then input the value for ibase in decimal.

For the length of numbers used in a computer, binary is cumbersome. Unless you really want to know the state of a particular bit, you use hexadecimal (a number system with base 16), which uses 1 symbol for 4 bits, and runs from 0..f (or 0..F)

binary  hex  decimal
0000 	0        0
0001 	1        1
0010 	2        2
.
.
1000 	8        8
1001 	9        9
1010 	a or A  10
1011 	b or B  11
1100 	c or C  12
1101 	d or D  13
1110 	e or E  14
1111 	f or F  15

When input to a computer is ambiguous as to its value, hexadecimal is represented as "OxF" or "Fh" (preceded by "Ox" or postceded by "h").

Here's conversion of hexadecimal to decimal using bash

declare -i result	#declare result to be an integer
result=16#ffff		#give result the value in decimal of hex ffff
echo $result		#echo $result to the screen
65535

#or all in one line
declare -i result;result=16#ffff; echo $result
65535

using bc

echo "obase=16;ibase=16; F+F" | bc
1E

	Note
	End Lesson 3. Spent some time in the first half of the class going through the two's complement exercises which I added after lesson 2. I asked the kids to try the following exercises for homework. They didn't do them so I started with the decimal/binary/hex table above and then worked them through the exercises below, at the start of the class.

Using any method

Programs don't usually do much arithmetic with integers. Integers are used as counters in loops and to keep track of the position in an executing program. Integers do come from digital sensors: e.g. images from digital cameras, digital audio, digital sensors. However most data, by the time it arrives at the computer, is reals.

In a 32 bit computer, an integer has a range of 0-4294967295 (2³², this number is referred to, somewhat inaccurately as 4G, but we've all accepted what it means - it's the 32 bit barrier).

#in bash
#binary
declare -i result;result=2#11111111111111111111111111111111; echo $result
4294967295

#hexadecimal
declare -i result;result=16#ffffffff; echo $result
4294967295

Numbers needing more bits than the machine's register size are called Long (or long), e.g. a 64 bit number on a 32 bit machine. Arithmetic on long numbers needs at least two steps, each of 32-bit numbers, and requires an "add with carry" (ADC) instruction (found on all general purpose computers designed to be able to load any code the user wants to install, but not found on special purpose computers exclusively running preinstalled code such as cell phones, routers, automotive computer chips). Here's how addition of long numbers works. Let's assume a 2bit computer and we want to add a 4bit number.

0010
1011+
----
????

First split the problem into pieces managable by the hardware (here 2 bits) giving us the right hand half (the least significant bits) and the left hand half (the most significant bits).

LH         RH       
00         10         
10+        11+
--         --
??         ??  

Next a word about addition and carry: When doing addition by hand, there is never a carry for the rightmost digit, but a computer has a carry bit for the rightmost bit which is set to 0 at the start of addition.

RH
10
11+
----
?? sum
?0 carry

step 1: right column, add two digits + carry digit. The carry to the 2nd column is 0.
RH
10
11+
--
?1 sum
00 carry

step 2: left column, add two digits + carry digit. There is overflow
RH
10
01+
--
01 sum
00 carry (with overflow)

The computer has a FLAGS register (32-bits in a 32 bit computer), which holds, in each bit, status information about the executing program, including whether the previous instruction overflowed, underflowed or set a carry.

The addition above overflowed, but the computer doesn't know if the bit is required for Long addition, in which case the overflow is really a carry. The computer stores the overflow bit in the carry bit in the flags register just in case. If the computer is doing a Long addition, the next step will ask for the carry bit. If the computer isn't doing a Long addition, then then the carry bit will be ignored (and will be lost).

Here's what the calculation looks like now (only the state of the carry bit is shown in the FLAGS register). The computer will first add the right most digits in its 2bit registers, using the regular add (ADD) instruction, which only adds the two numbers and the information setup in the carry input to the adder.

before 1st addition

LH         RH        FLAGS
00         10         ?
10+        11+
--         --
??         ?? sum 
?0         ?0 carry

after 1st addition

LH         RH        FLAGS
00         10         1
10+        11+
--         --
??         01 sum 
?0         00 carry

Because of the overflow, the FLAGS register is now 1. The computer has been told that it's doing the 2nd step in a Long addition. It uses the "add with carry" (ADC) instruction, which transfers the carry bit in the FLAGS register to the adder, and then does a normal addition.

2nd addition. first step, copy carry bit from FLAGS to carry input for LH

LH         RH        FLAGS
00         10         1
10+        11+
--         --
??         01 sum 
?1         00 carry

2nd step, add digits and carry digits for LH numbers

LH         RH        FLAGS
00         10         1
10+        11+
--         --
11         01 sum 
01         00 carry

we now read out the sum digits

11         01

giving the required answer of 1101

You can chain addition to any precision (on a 32-bit computer, to 64, 96, 128-bits...) Standard calculations rarely need more than 64 bits, but some people want to calculate π to billions of places and this is how they do it.

Long arithmetic is slower than regular arithmetic. You don't ask for Long operations unless you know you need them.

	Note
	End Lesson 4

If we wanted negative integers, how would we do it? Pretend you're a 1 byte computer and you need to represent -1. You can do this by finding out the number which added to 1 gives 0.

00000001
????????+
--------
00000000

The answer is 11111111₂, 255₁₀or FFh (note: computers need overflow to work).

00000001
11111111+
--------
00000000

You've seen the computer version of -ve numbers before. They're called what ^[30] ? They are the (-ve of the number + the base) (in a 1 byte computer the base is 256).

What is -2₁₀ in binary, hexadecimal? ^[31]

How do we know whether 255 should be interpreted as -1 or 255?

The level of primitive data types, is one level above a byte. Your program keeps a record of the primitive data type that each particular byte represents. When you write your program, your code will have descriptors stating whether this integer will have +ve only values (called an unsigned int) or both +ve and -ve values (called a signed int). Some programming languages will have already decided that you'll be using a signed int and you won't have any choice.

If you have a signed int, then integers with high order bit=1 are -ve while those with high order bit=0 are +ve.

binary     hexadecimal decimal
00000000       00         0
00000001       01         1
.
.
01111111       7F       127
10000000       80      -127
10000001       81      -126
.
.
11111100       FC        -4
11111101       FD        -3
11111110       FE        -2
11111111       FF        -1

Linux runs on 32 bit computers. What values is represented by 32x1's (or FFFFFFFF) in Linux bash. On a 32 bit machine we might expect this to be -1.

declare -i result;result=16#ffffffff; echo $result
4294967295

This is not a negative number. What's going on? Let's try a 64-bit number (just for reference, the biggest number that can be represented by 64 bits is 2⁶⁴=18,446,744,073,709,551,616=18.45*10¹⁸).

declare -i result;result=16#ffffffffffffffff; echo $result
-1
declare -i result;result=16#7fffffffffffffff; echo $result
9223372036854775807
declare -i result;result=16#8000000000000000; echo $result
-9223372036854775808

bash in Linux rolls over to -ve numbers half way through the 64 numbers: the integers in bash on Linux on 32 bit machines are 64-bit signed integers.

	Note
	see long_numbers for 64-bit math on 32 bit machines.

Here is the range of values that various sized integers (int) can represent. To give an idea of the relative sizes, the table shows the time (stored as seconds) represented by that integer.

			8bit		16bit		32-bit		64-bit	
unsigned		0-255		0-65535		0-4294967296	0-18446744073709551616
signed			-/+127		-/+32767	-/+2147483647	-/+9223372036854775808

unsigned time		4mins		18hrs		136yrs		584,942,417,355yrs
signed time		2mins		9hrs		68yrs		292,471,208,677yrs

You'll see these numbers often enough that you'll start to remember them. For the moment be prepared to see numbers pop up that you've not previously seen much of.

The Y2.038K problem. In 1999 the computer world was in a flurry: programmers who'd prepresented years in their dates using 2 digits, realised that the year 00 would represent the year 1900. Paychecks would not be processed, elevators would stop working at midnight, trapping thousands (if not millions) of innocent people, and planes would fall out of the sky (except Japanese planes). A bigger calamity could not be imagined. The world's bureaucrats heroically spent millions$ of taxpayer's and consumer's money to prevent certain disaster. On 01 Jan 2000, none of the predicted misfortunes occured, for which we must thank the selfless and unacknowledged taxpayers and consumers of the world.

Unix represents time as a signed 32-bit integer of seconds, starting 1 Jan 1970, this date itself a major blunder ^[32] . If in Jan 2038, you're still using a 32 bit OS (not likely for desktop machines, but quite possible for embedded devices sitting in computer controlled machinery, which rarely need 32 bits, much less 64 bits), Unix time will overflow in Jan 2038. If in Jan 2038, your computer controlled refrigerator stops ordering food, it will be because the refrigerator is asking for food to be delivered in 1970. Jan 2008 was a good time to take out a 30yr loan for 250,000$; your monthly payments would be -1,600$ (a comment from slashdot in Jan 2008. http://it.slashdot.org/article.pl?sid=08/01/15/1928213).

A 32 bit computer can generate how many different different integers ^[33] ? This computer then is capable of generating that many different integers. If you use the integers to be labels or addresses, you can label that many different items.

A 32 bit computer addresses its memory in bytes using address lines in the hardware. The computer has to read/write a byte as one unit; it can't address individual bits in a byte in memory - it has to read or write the whole byte. Once the computer has read the byte into a register, then each bit within the byte is separately addressable. What is the maximum number of bytes a 32 bit computer can address, and how much memory is this ^[34] ?

Not too long ago, microprocessors had 4kbytes of memory. Now for many applications, computers need more than 4Gbytes of memory. These applications all have to be run on 64-bit computers. But if you wanted to increase the amount of memory available to 32-bit computers, how could you do it ^[35] ?

In fact most data and instructions in a 32-bit computer are 32 bits and are fetched 32 bits at a time, so changing the addressing to 32 bits would not be a big change. A char would now have to be represented by the byte of interest, followed (or preceded) by 3 empty bytes. The instructions that work on chars would have to mask off the 3 empty bytes. Since chars (in most programs) are a small fraction of the data held in memory, these unused bytes would not cause too much wasted space in memory.

I don't know why 32 bit computer manufacturer's didn't go to 32 bit addressing. Possible reasons might be

Harddisks read/write data to/from independantly addressable blocks (i.e. the computer cannot address individual bits or bytes within a block; it has to address and then read or write the block as one indivisible unit). Let's assume a blocksize of 512bytes. If the computer wants 1 bit off the harddisk, the computer has to know the address of the block, read the whole 512bytes into registers, manipulate the 1 bit and the write the 512byte block back to disk. What's the maximum number of blocks a 32 bit computer can address on a harddisk and how much can be stored on a disk with blocksize=512bytes ^[36] ?

If you wanted to put a bigger disk on a 32 bit computer, how would a harddisk manufacturer do it ^[37] ?

Harddisk manufacturers originally started with block size of 512 bytes and have incremented the size of the blocks continuously over the years. I don't know why hard disk manufacturers can change the blocksize whenever they want and not have programs fail, while at the same time microprocessor manufacturers have not been able change from 8bit addressing to 32 bit addressing. Possible reasons might be

My guess as why memory addressing stayed constant at the 1 byte granularity over decades, while disk block size (granularity) increased from 512 bytes to 8192 bytes, is that going to larger block size forgoes the chance of addressing a 512 byte block, which doesn't cause any problems (you don't ship disk blocks off the local machine), but choosing 32-bit addressing forgoes the chance to address a byte, which you do want to be able to move around (including sending to other machines or peripherals).

IPv4 internet addressing uses a 32 bit unsigned integer (IP) to uniquely identify each network device (e.g. no two devices/computers can have the came IP). How many devices can be simultaneously on the internet using an IPv4 IP ^[38] ?

The internet game "World of Warcraft" (reported in slashdot, http://games.slashdot.org/games/08/01/19/1321241.shtml) uses a 32 bit integer to store the amount of gold that players earn. Some players have reached that limit and can no longer receive gold from other players. If these players have 1c of gold for each bit, what is their wealth in $ ^[39] ?

The world's telephone system carries voice data in 8bit form i.e. it converts your voice into bytes, each byte representing the amplitude of your voice at any instant. How many different amplitude levels can be expressed using 1byte ^[40] ? Since the phone system uses 8bits, it's simple to send bytes from computer data across phone lines. Hifi audio is usually 12bits (how many levels is this? ^[41] ) which has less noise than 8bit audio.

Digital cameras, and computer monitors use 8 bits to represent the intensity of each of the 3 primary lights of a picture; red, green and blue. This turns out to be more levels than the eye can differentiate, but not much more (the eye doesn't see the edge between one intensity and the next and only sees a continuous change in color).

A color picture from an 8 bit digital camera can have how many different colors ^[42] ?

Matte print (e.g. books, newspapers) can only reproduce about 100 levels of light, but glossy print (e.g. in fine books and magazines) can reproduce about 180 levels, which is why expensive advertisements are run in glossy magazines.

The human eye can accomodate a range of light from nightime to midday on a cloudless day, a range of 10⁷ in intensity (I think). The eye can see features in a face in shadow and in the face of a person standing next to them in full sun, but an 8 bit digital camera will, according to the exposure, only see the face of one, while the other will be washed out (either dark, or light). To help in situations of high contrast, expensive digital cameras record 12bits, allowing range compression and expansion. These photos are post-processed to reduce the range to 8bits for display (or printing) but keeping constrast in both the light and dark areas (i.e. you can see the features of a face in shadow and a face in bright light in the same photo).

What other things could we use 32-bits for: How many -ve numbers could we have ^[43] ? How many prime numbers could we address ^[44] ?

Note

End Lesson 5. At the start of the next class, I went back over the number of items that a 32 bit computer could represent. The students had forgotten the number of integers that could be represented by 32 bits (not only the 4G value, but the concept of there being a limit associated with 32 bits). I went through the number of integers, computers on the internet, blocks on a hard disk etc again. The seemed to remember the concept after a few examples. My partner reminded me that you have to tell a student 3 times before you can expect them to start to remember a fact.

Now that we're at the level of primitive data types, we can use a language like python.

fire up python; you will be running python in immediate (interactive) mode, where it interprets each instruction one at a time, and returns to the prompt. (In normal operation a program keeps executing till it has to wait, say for a keystroke from you.)

You will get the python ">>>" prompt

Python 2.4.3 (#1, Apr 22 2006, 01:50:16) 
[GCC 2.95.3 20010315 (release)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> 

	Note
	The following examples are based on Chapter 1 of the LiveWires tutorial "Introducing Python".

try a few subtractions, multiplications and divisions on your own.

>>> 12 + 13
25

>>> 123456 * 3
370368

How about this one?

>>> 7/3
2	#7,3 are integers. You're asking to do integer arithmetic. You get an integer answer.
>>> 

>>> 7%3
1	#'%' is modulo (remainder)
>>> 

>>> 7.0/3
2.333333333333333 #as long as one of the numbers is real, the answer will be promoted to real

you'll learn about real numbers soon.

	Note
	This went over the kid's heads, so I skipped to the next section (they don't need to know this right now)

Most programs (including python) use the machine's native libraries (e.g. math, string) (which are usually written in C). (No-one writes a library when a well tested one is already available.) The size (number of bits) for various primitive types in python will then depend on the native libraries. The documentation for python says there are two types of integers (see Numeric Types http://docs.python.org/lib/typesnumeric.html).

The Python documentation doesn't tell you the sizes for these two types of integers for any particular platform: you're supposed to be able to work it out yourself. What's the largest plain integer that python can represent? (for likely numbers, look at the table in integer_range) i.e. is it 32 or 64 bit? You won't have to remember the range of integers in python, but you'll need to understand enough about a computer to figure it out. You also should not be surprised if numbers become Long when they become big enough. (In the following, remember 65536 fills 2 bytes. For compactness, I'll use hexademical to illustrate what's happening.)

Python has no trouble representing any size integers. Here are some integers from 16-256bits (Long integers end with "L").

>>> 65536-1			#16bit FFFF
65535
>>> 65536*65536-1		#32-bit FFFFFFFF
4294967295L
>>> 65536*65536*65536*65536-1	#64-bit FFFFFFFFFFFFFFFF
18446744073709551615L
>>> 65536*65536*65536*65536*65536*65536*65536*65536-1	#128-bit  (32 Fs)
340282366920938463463374607431768211455L
>>> 65536*65536*65536*65536*65536*65536*65536*65536*65536*65536*65536*65536*65536*65536*65536*65536-1 #256bit
115792089237316195423570985008687907853269984665640564039457584007913129639935L

From the above output, the 16bit number 65536 is a plain integer, but the 32 bit number is Long. To calculate 65536*65536-1, we would first have had to calculate the intermediate result 65536*65536 which would be a Long (needs 33bits). If you subtract 1 from a Long, you still have a Long (even if it could be represented as a plain), so FFFFFFFF could be a plain integer, but we wouldn't have found out doing it this way. Let's look around for a 32 bit limit. Remember that only half of the integer range is used in signed integers, so let's look at half of a 32-bit number.

#here we would need 33 bits to handle the multiplication overflow 
#so we already know that the answer will be a L number.
>>>65536*65536-1	# FFFFFFFF
4294967295L
#what's half of 65536?
>>> 65536/2 		# 8000       
32768
#The result will be 80000000H, which if a signed integer will be a -ve number
#It looks like python promotes the integer to a L
>>> 32768*65536		# 80000000
2147483648L
#Let's get below 80000000H to 7FFF0000H. 
Yes it's a plain integer
>>> 32767*65536
2147418112
#Let's try 7FFFFFFFH. Yes it's a plain integer.
>>> 32767*65536+65535
2147483647
#just checking 80000000H again, this time approaching from below.
#It's L
>>> 32767*65536+65536
2147483648L

The largest plain integer represented by python is 7FFFFFFFh or 2147483647.

This process above, of poking numbers (or different piece of code) into a program to see what it does is called noodling. It's a good way to learn.

What's the -ve most plain integer in python ^[45]

Python uses a signed 32-bit integer to represent plain integers.

We need to represent the characters in the alphabet, so the computer can type them on the screen and receive them from the keyboard. We need upper and lower case (52) + numbers (10), plus some punctuation and printer/screen control characters (move up/down/left/right, carriage return, line feed, end of file).

abcdefghijklmnopqrstuvwxyz
ABCDEFGHIJKLMNOPQRSTUVWXYZ
0123456789
,. !@#$%^&*()-_[]{};:'"<>;/?`~

This is more than 64, but less than 127. This number of characters requires 7 bits. A regular 8 bit byte is used with the top 127 generally unused. The mapping between the 256 possibilities in a byte and the symbols displayed above, as mandated by the USGovt, is called ASCII.

A table of ascii characters and their binary/decimal/hexadecimal equivalents is at wiki, ASCII (http://en.wikipedia.org/wiki/ASCII). The table of printable characters (http://en.wikipedia.org/wiki/ASCII#ASCII_printable_characters). shows that in ASCII, the characters are in alphabetical order.

	Note
	unlike some other character sets e.g.EBCDIC http://en.wikipedia.org/wiki/EBCDIC originally devised as an extension of Binary Coded Decimal (BCD) http://en.wikipedia.org/wiki/Binary-coded_decimal needed to handle money.

A table which better illustrates the hexadecimal organisation of ASCII is ASCII Chart and Other Resources (http://www.jimprice.com/jim-asc.shtml#table). (A slightly fancier table ASCII Table and Unicode Characters http://ascii-table.com/).

The numbers are 3hex+number. This allows easy conversion of a character representing a number into a number (you mask off the left 4 bits and add the right 4 bits into the output number).

bash converts the hexadecimal representation of a character to its ascii symbol using this command

echo $'\x41'
A

	Note
	The "41" is hex and the 'A' output is a char (not a hex number). The rest of the command is obscure bash magic.

although you're never likely to have the octal representation of a char, bash converts the octal representation of a character to its ascii symbol using this command

echo $'\41'
!

	Note
	418=84+1=3310=21h echo $'\x21' !

How many letters down the alphabet is the character 'B' (try it at the prompt) ^[46] ? How many letters down the alphabet is the character represented by hex '51' ^[47] ? Knowing that the hex for 'A' is 41h, figure out the hex for 'Z' and then try it ^[48]

To change between upper and lower case, the 6th bit in the byte is flipped. What change in value (as an integer) does flipping the 6th bit represent ^[49] ?

echo $'\x5A'
Z

echo $'\x7A'
z

In a program. to differentiate a character from a number or variable do this

'c'	char
 c      the variable named c #better to use a longer descriptive name, eg computer_name
 7      the number 7
'7'	the character 7

Computers can scan text to test which characters are letters (A-Z,a-z), which are numbers (0-9) and which are punctuation. The computer can match characters (e.g. is the character an 'A' or a '9'?).

Note

Every keystroke on your keyboard is a character. If you type "3.14159" on the keyboard, the computer accepts it as a series of characters. If you want this to be a real, then you have to explicitely tell the computer to convert the string of characters into a number. If the computer asks you to input a number at the keyboard, your keystrokes will first be put into a string buffer and later your program will have to convert the string to a number. If you have "3.14159" displayed on the screen and swipe it with your mouse and put it into a buffer, it will be in your buffer as a string of characters. All normal input and output on a computer is characters and strings e.g. keyboard, screen, printer. (Some programs exchange data as binary, but you have to set that up.)

	Note
	The representation of real numbers take a bit of explaining. You don't need to understand how the are represented to use them (we'll do that later - look in Real Numbers)

"real" numbers (also called a "floating point" numbers) in computing are numbers like 3.14159 - anything with a decimal point. You do arithmetic on them.

>>> 3.0*6.0
18.0

You can mix integers and reals - the computer handles it for you, promoting the integer to a real.

>>> 3.0*6
18.0

Be careful how you mix integers and reals. The computer first evaluates (7/2), not knowing that it will next do an evaluation of a real.

>>> 7/2*5.0
15.0

A minor rearrangement of this code gives

>>> 5.0*7/2
17.5

Minor editing of this code makes a big difference in the output (one is correct and one is not). Code where a minor edit (like rearranging the order of a multiplication, which should not change the result) gives a different answer. Such code is called fragile code. Someone (maybe you), years later, could be working on your code and not see the code bomb and will rearrange the line to trigger the bomb.

You should practice safe programming. When mixing integers and reals, explicitely promote the integers to reals, and don't expect the computer to do it for you. Don't rely on rules of precedence too much. Use brackets to make code clear for the reader. This is how the code should be written.

>>> (5.0*7.0)/2.0
17.5
>>> (7.0/2.0)*5.0
17.5

a string is a series of characters delineated by a pair of double or single quote characters; e.g. "happy birthday!", "1600 Pennsylvania Ave", "temperature=32F, pressure=29.90in", 'I am "late"'

In principle it's possible to operate on strings as arrays of characters, but strings are the dominant form of input and output on a computer (all computers and people can read them), so all languages have instructions to search, match, read and write strings.

In situations where enormous amounts of data are involved, which will never be read by a human and only ever read by another computer (mapping, photos, MRI data), then data is written in some compact binary (and likely platform specific) form. You'll still need a team of programmers to write the code to allow each new generation of computer to read and write that format.

	Note
	End Lesson 6

Until you get used to the rules, and familiar with the error messages, you will have to check each time what you have. Some interconversions are done without asking and others have to be invoked explicitely.

Here's some operations on numbers

>>> variable=3.14159
>>> 3*variable		#since you can multiply, it must be a number, the integer is automatically promoted to real
9.4247699999999988

			#how does python know 3.14159 is a number and not a string?
			#see below for variable=3.14159q	

>>> print variable	#you can print numbers
3.14159

			#the + operator joins strings
	                #but you can't print a string and number at the same time using a +
			#the error message is helpful
>>> print "the number is " + variable
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
TypeError: cannot concatenate 'str' and 'float' objects

			#you have to turn the number into a string
>>> print "the number is " + repr(variable)
the number is 3.1415899999999999

			#you print numbers using a ','
>>> print "the number is", variable
the number is 3.14159

Here's some operations on strings

>>> variable="3.14159"
>>> 3*variable 
'3.141593.141593.14159'	#same string 3 times. 
	                #not real useful, and probably not what you want.
	                #no other language has this. 
	                #if you really need to do this, 
	                #use a construction common to all languages
	                #or no-one will be able to maintain your code.
			#
			#what would have happened if you'd done the following?
			#variable="3.14159 " 
			#or
			#variable="3.14159q"
			#or
			#variable=3.14159q	

>>> 3.0*variable	#the error message is not helpful
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
TypeError: can't multiply sequence by non-int

>>> 3.0*float(variable) #float() converts a string to a real 
9.4247699999999988

>>> print variable      #you can print a string.
3.14159

Note the unhelpful error message resulting from the command 3.0*string. You should be prepared for error messages to be wrong and have to go to google to find out the real problem (don't expect the right answer to be available there either, it usually will be, but you'll still have to figure some of them out yourself). The interpreter knows that it has a string and not a number and could have told you. Unfortunately this is part of computing - it's always been this way. I wonder if the messages are designed to raise the cost of entry to being a programmer. The error messages from gcc, the GNU C compiler, have improved dramatically in the last 10yrs.

	Note
	One of the students commented that the error messages were like having life boats that don't work.

The data types described so far are found in most languages. Others are common, but on knowing these four, the new ones will be easy to use when we need them.

What primitive data types do we know now ^[53] ?

	Note
	If you've been using an editor to do your binary math examples above (rather than pencil and paper) and you're happy with it, then you can skip to saving_files.

From the NCSA discussion mailing list (and other places) comes these suggestions for windows/Mac editors other than vi and emacs

When I had to code on windows (worst 3 days of my life), I installed cygwin and could pretend I was working on a unix machine. This is a reasonable approach if you already know the unix tools.

For this class pick any editor you like. Be prepared to try a few editors before finding one you like. Everyone is different and everyone likes a different editor. I like mine really simple. It seems that others like theirs complicated.

If you find yourself programming on a long term basis in a unix environment, and you expect you'll be sitting in front of a machine that someone else has setup (e.g. at work), that you should learn one of vi or emacs. It's quite disappointing to realise that the program you spend most of your time with, the editor, has such a bad interface and that no-one has designed an editor for unix like the much nicer ones available on the much dispised Microsoft platforms.

pick a name for your python files directory e.g. class_files.

Your files need to be in a place where no-one else is likely to run them accidentally. Until your files are well tested and of general use, they should be kept in your directories.

Make the file executable.

chmod 755 hello_world.py

At the command prompt, run it

# python hello_world.py
hello world!

Congratulations, you are now officially a computer programmer.

The above command works because python will have been installed in your $PATH.

You can invoke the interpreter from within your program using the shebang convention. Here's the new code (the "#!" is called a shebang)

#! /usr/bin/python 
print "hello world!"

You run it from the command line like this

# ./hello_world.py
hello world!

In windows you'll be in a command box at \Documents and Settings\UserName\class_files (or \My Documents\class_files which is the same thing). The python install sets up the registry so that you can directly execute the python program. You only need to do

hello_world.py

Congratulations, you are now officially a computer programmer.

You can execute the program by clicking on the filename using windows explorer, but the output window will open and close too fast for you to see what happened.

The unix execution options are still available, but you don't need them in windows

All computer languages allow comparison of the contents of variables with the contents of other variables (or with numbers, strings...). Tests are

Following the test, code execution can branch in various ways e.g. (here shown in pseudo code)

Note: python doesn't have this problem:

Here is a common coding bug. The comparison operators are unique to comparisons, but the "==" and "=" operators are similar, at least visually.

x=0 means "assign x the value 0"

x==0 means "test if x has the value 0"

In a conditional you want

if (x==0) do something else do something else

If instead you do

if (x=0) do something else do something else

then x will be assigned the value 0. The language must know whether a command succeeded (allowing execution to continue, or to bail out with an error) and since assigning x=0 will always succeed, then execution will branch to the "do something" branch, independant of the original value of x.

Since you're unlikely to ever want to execute the instruction

if (x=0)

the language should trap this construct as a syntax error. None of them do, at least till python came along. This is why rockets continue to blow up, cars have different bumper heights and we had lead in paint for so long.

Here's the official Python Tutorial, section 4.1 If Statements (http://docs.python.org/tut/node6.html#SECTION006100000000000000000)

Note

I live in one of the few (the only?) country that uses Fahrenheit, miles, lbs... (and the country's currency is loosing value, while the economy is only surviving by exporting jobs to our competitors.) If you live in the rest of the world, please substitute other values in these examples. I also live in a part of the world where you need an airconditioner in the summer and heater in winter (I lived quite happily in another part of the world for the first 30yrs of my life without these things and had some trouble adjusting to living indoors). I still wonder why people settled in such places when they didn't have to.

Here's the python syntax for a single branching conditional statement.

if x < 0:
	print 'Negative number found'

fire up your editor to code up the file temperature_controller.py

Here's my code ^[59] .

Here's python code for a two branch conditional.

if x < 0:
       print "Negative number found"
    else:
       print "non negative number found"

starting with your current temperature_controller.py, add code to output the string "no action taken" if the temperature is 80 or below. Check your code with values of temp above and below 80°.

Here's my code ^[60] .

Here's another python conditional construct

if x < 0:
      print 'Negative number found'
elif x == 0:
      print 'Zero'
else:
      print 'Positive number found'

Use this construct to add a branch to temperature_controller.py which turns on the heater at 60° or below and outputs a message describing its action.

Here's my code ^[61] .

Change the preceding code to open the windows if the temperature is 61-80°. Test your code to confirm that all branches execute at the expected temperatures. Here's my code. ^[62] .

In our current version of temperature_controller.py, the string "the temperature is " is in all branches, so will be output no matter what. In this case the string could be output once, below the declaration of temp. Here's a print statement using a trailing ',' that doesn't put a carriage return at the end of the output line.

name = "Homer"
print "my name is", name,

Use this code fragment to only have one copy of the string "the temperature is" in the code. Here's my version ^[63] .

Note that you can't get the period correctly spaced after the temperature. Python thinks it knows what you want better than you do and outputs a gratuituous blank that you didn't ask it to output.

Formatted output allows you more control over your print statements. Use this code fragment to produce a better output for temperature_controller.py (the %d says to put the decimal value of the argument in that place in the string).

>>> temp = 65
>>> print "the temperature is %d." % temp
the temperature is 65.

Here's the improved code ^[64] . Note that python still outputs a gratuituous blank when it executes the comma.

	Note
	End Lesson 8

Computers are ideally suited to run the same piece of code over thousands, if not millions of pieces of data. Doing the same thing over and over again is called iteration and all computer languages have instructions for iteration. The part of the code that is iterated is called a loop. The philosophy behind python's iterative commands is different to all other languages.

In python if you create the list with range(), you have to create the whole list before you start the loop. In the 2nd method, you create one item at a time. In python, if your list is larger than can be held by the memory in your computer, then your program won't run. We'll find an example of this later when calculating π. There we'll use xrange() which generates one item at a time, rather than a list.

Here's example python code which iterates over a list of numbers (the '[' and ']' delineate a list).

#!/usr/bin/python

for x in [0,1,2,3,4,5]:
	print x,

Code this up as interation_1.py and check that you get reasonable output. What happens if you leave off the ',' at the end of the print statement ^[65] ?

The body of this loop has one instruction (the print() statement). The length of the loop isn't so obvious: no-one is pedantic about this - including the required blank line at the end, you could say that the loop is 3 lines long, or maybe 2 (the for statement and the print() statement) or maybe 1 (the print() statement).

Note

If you run this code at the python prompt (>>>), you will need a blank line after the print() statement, to let python know that the indentation has changed back to the previous level of nesting (here the left most column). Compare this piece of code #!/usr/bin/python for x in [0,1,2,3,4,5]: print x, print "hello" with this piece of code (the only change is the indentation for the print "hello" statement). Predict the output of the two pieces of code and then run the code. #!/usr/bin/python for x in [0,1,2,3,4,5]: print x, print "hello"

Change the code in iteration_1.py to

Here's my code ^[66] .

Warning

x (the variable whose value is being set by the for statement) is called the control or loop variable. It's only purpose is as input to the loop statements. The loop variable should be treated as a readonly variable. (i.e. you shouldn't try to change it) i.e. do not put the control variable on the left hand side of any instructions in the loop. DO NOT DO ANY OF THESE x = x + 2 x = 2 x = "hello" Changing the loop variable is unsafe programming. If you're debugging a long loop, you don't want code in the middle messing with the control variable. In some languages (fortran) you aren't allowed to change the loop variable, or it will change in unpredictable ways. (In python, in this loop above, changing x wouldn't cause any great problems, other than being unsafe programming.) If you want to do something with the control variable, make a copy of it (with say y = x) and do whatever you need on the copy y.

Using iteration_1.py as a template, write iteration_2.py to output and sum the squares of the numbers 0..5. Here's my first version ^[67] . Notice that the calculation of x*x is done twice. You shouldn't ever do a calculation twice - it takes too much time. If you need the result in two places, you do the calculation once and save the result in a variable. Here's a better version ^[68] .

	Note
	some new syntax: these are identical y = y + 3 y += 3 other operators which use the same syntax are -=, *=, /=, %=

Modify iteration_2.py to use this new syntax. Here's my new version ^[69] .

Most often, you aren't iterating over an irregular or random list of numbers. You usually iterate over a well behaved and ordered list. Rather than enter a list by hand, risking a mistake, python (and a few other languages) have functions to generate lists. Python's command is range() and true to python form, it behaves differently to the equivalent command in other languages.

Here's the way everyone else does it

range(start, end, [step]) #[] indicates an optional parameter, default value here is 1

range (0,10)
 0,1,2,3,4,5,6,7,8,9,10

Here's python's version

>>> range(10)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

The argument (parameter) 10, says to generate the default list (i.e. starts at 0, stepping/spacing of 1) and stop at the element before 10 (simple huh?). In this case range() produces 10 numbers, as you might expect, but see below for further developments.

For fun, try a few other sets of parameters for range().

Copy iteration_2.py to iteration_3.py and use range() to generate the list of numbers from 0..5. Here's my code ^[70] .

Well that was confusing, if you want to generate a list from 0..5, you need an argument of 6 for range().

Here's some more examples using range().

#as above, starting at 0
>>> range(0,10)		
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
#starting at 1, ending element before 10
>>> range(1,10)		
[1, 2, 3, 4, 5, 6, 7, 8, 9]
#starting at 1, in steps of 3, ending at element before 10.
>>> range(1,10,3)	
[1, 4, 7]

Using the previous code as a template, sum the squares of all the numbers from 0..100. Then sum the squares of all the numbers from 0..100 that are divisible by 7. Here's my code ^[71] .

	Note
	End Lesson 9

fencepost error (http://www.eps.mcgill.ca/jargon/jargon.html#fencepost%20error)

We've run into the fencepost error before. when I asked you to march 10h places down the alphabet from 'A'. The correct answer was 'Q', but someone got 'P'.

You want to build a fence 100ft long, in sections of 10'. How many sections of fence will you be building ^[72] and how many fenceposts do you need ^[73] ?

Same question again: if I want to sum the squares of the n consecutive integers, how many iterations (of the loop) do I need ^[74] ? If I want to sum the squares of integers 0..n, and how many iterations do I need ^[75] ?

This will fool you over and over and over, from the start of your programming career to the end. Anytime you iterate (go) from a start boundary condition, over a range, to an end boundary condition, you'll have to do a test count, to determine whether you're counting fenceposts or fences. (When you're programming, it's not always obvious whether a variable is a fence or a fencepost.) If you've got it wrong and you only have one iteration to do, the computer will think you only have to do 0 iterations. Your loop won't be executed and you'll think that there was something wrong inside the loop, that it didn't add (or whatever) your number.

First Steps to Programming (http://docs.python.org/tut/node5.html#SECTION005200000000000000000) shows a while loop.

A while loop has this structure

determine condition
while condition is true
	do a
	do b
	.
	.
	determine condition again

#code executed when condition becomes false

The while and for loops can be logically equivalent. These two produce the same output:

The while loop requires a little more setting up, but unlike the for loop which is good for regular input, the while loop can be setup to handle just about anything (here pseudo code accepts a name)

print "enter your name followed by a carriage return"
name = ""
char = getkeystroke()	#I made this function up
# != is not equal - the loop is entered for all chars except a carriage return
# if the char is a carriage return, which instruction is executed next?
while (char != carriage return):
	char = getkeystroke()
	name += char

print name

The essential features of the while loop are

Here's a program that has no functions. It gives a greeting. Call it greeting.py

#! /usr/bin/python
""" greeting.py
    greeting without any functions
"""

#---------
#main()

print "hello user"

# greeting.py ---------------------------------

The code includes documentation and white space to separate the logical sections (in this case, the documentation and the code). The dividing line and "#main" are to let people know where the code starts executing. This is optional, but some delineator is needed to help you (and others) read your code later. Remember if you're being paid to write code, or work in a group, others will not be impressed if they have to struggle to read it.

Modify this code to say "hello, this code was written by (your name)" ^[76]

If we want to call the greeting from several places in the code, we turn the printing of the greeting into a function. Call this file greeting_2.py

#! /usr/bin/python
""" greeting_2.py 
    greeting as a function
"""

#------
#functions

def print_greeting():
	"""prints greeting
	"""

	print "hello user"

# print_greeting---------

#main()

#call the function
print_greeting()

# greeting_2.py ---------------------------------

What's going on:

Indenting is a common coding practice to show the reader logical separation of code from the surrounding code. In all languages, except for python, the indenting is only for the reader and is ignored by the compiler/interpreter. In other languages, pairs of "{}" (braces, squigglies) are used to delineate blocks of code; pieces of code that have their own scope (variables created in a block don't exist after the block exits).

Python uses the indenting to show both the start/end of functions (and blocks), and for the reader to show logical sections of code. Proper indenting is required for the code to run, and is not just for clarity to the reader. The Style Guide for Python Code (http://www.python.org/dev/peps/pep-0008/) requires 4 spaces for python indenting. For compatibility with old code, python will accept a tab or 8 spaces (but not mixtures of a tab in one line and 4 spaces in the next).

The problem with python indenting is that python has syntax errors which are silent to inspection by eye (you can have mixtures of tabs and spaces which are the same on consecutive lines, but all you'll get from python is "syntax error" on the first non-white space character). If code has an error in it, you should be able to see it. (This sort of feature is called "broken by design".) This is a rocket waiting to blow up.

There are many style of indenting for readability and code nesting with curly braces Indent style (http://en.wikipedia.org/wiki/Indent_style). Each style has its adherants, and without any tests as to the best style, or any interest in running such tests, claims of superiority of any method are specious and cast doubts on the credibility of the claimant. There is no doubt that people can best read the style they're used to (just like they can best use the keyboard layout they're used to). People working in a team will usually be required to use the style that the leader is most familiar with.

It turns out that (without knowing it) I adopted the Allman style for curly braces, using a tab (1 key stroke) for indenting. I will be using a tab for indenting in the code here. You can use whatever your python interpreter will accept. Be aware that a tab when mouse swiped in X, will appear as 8 spaces in the target code (after a mouse swipe, I have my editor replace all occurrences of 8 spaces with a tab).

To help you read the code, you adopt an indenting style so that the indenting level indicates the block level. The compiler doesn't care if the indents and blocks don't match, the compiler looks at the block level. If your indenting and blocks don't match, you'll may have logical control problems (the program will do something, but it may not do what you think it should do). Some editors match pairs of "{}" so that you can check your blocks. In python, indenting is blocking. Editors don't show matching indents and until you're used to the python indenting, you'll run into apparently intractable problems. Here's an example:

#!/usr/bin/python

output_string = ""
sum = 0
for i in range(0,10):
	output_string += str(i)
	for j in range(0,i):
	        sum += j*j

	#lots of lines of code

	        #including multiple levels of endenting

	                #like these lines





	 print output_string

Here's the output

# ./indent.py
  File "./indent.py", line 18
    print output_string
	              ^
IndentationError: unindent does not match any outer indentation level

The problem is that the line "print output_string" is a tab followed by a space (2 indents) instead of a single tab.

Despite the advice of the python people, I think it's safer to use tabs for code indenting

If your code is nested to great depth (not a great idea), deeply nested code will be diplaced to the right across the page.

In python, function declaration (the tuple of def name (parameters):) and definition (the lines of code that do the work) occur together. Functions often call other functions. If this second function has not yet been declared, the language/interpreter/compiler doesn't know what to do with the name.

Compare these two programs.

#functions declared in order
def fn1:

def fn2:
	call fn1 #fn1 already declared, no problems

#main()
	fn2
#---------

#functions declared out of order
def fn1:	
	call fn2 #don't know fn2, will have to look for it

def fn2:

#main()
	fn1
#---------

All languages have mechanisms for handling calls to functions that aren't known. You can list all your functions in any order you like (alphabetical by name, the order you wrote them, grouped by functionality), at the top of the file (the usual place) or at the bottom of the file (some people do this).

In many languages, particularly those designed to write code with thousands of files in a large application, the code writer is required to declare functions in a separate header file, which is read before the code file is read, so that all functions are declared before any code is looked at. Separating declaration from definition allows the compiler/interpreter to process files singly. This helps when you're working on one file in a large project.

Interpreted languages (like python, perl) keep track of the function calls they don't know about and exhaustively search all files in their pervue - it will search the Module Search Path. (http://docs.python.org/tut/node8.html#SECTION008110000000000000000) - in the hopes of finding the declaration. This can take a while, but what the heck, it makes life easier for you (let the computer work hard - that's it's job - it's not your job to make life easier for the interpreter, it's the interpreter's job to make life easy for you).

However for small projects, worrying about this sort of thing is a hill of beans. You're probably looking at the screen for 90% of the time trying to figure something out and you don't care if the computer is scanning files for an extra few seconds. Some people have all their functions at the bottom of the file, so the interpreter collects many calls to unknown functions before it gets to the function declarations below. These people are just as productive as the people who have their function declarations in order at the top of their files.

Usually a function needs information from the calling code. The information is passed as parameter(s).

Here we're also going to get input from the user. (Getting data directly from the user is not easy. It's simpler to get it from a file. Here we're going to get data from the user.)

Note

Let's say we want to get the user's name and print it on the screen. Getting input from users is fraught with problems and makes life difficult for the programmer. See Conversing with the user (http://www.freenetpages.co.uk/hp/alan.gauld/tutinput.htm) and Saving Users From Themselves (http://linuxgazette.net/issue83/evans.html). they'll enter "y" or "n" or their phone number. they can enter a string which when printed to your screen will erase your hard disk, so you have to check that the the string they enter is safe (this is a lot of work). It's a waste of a computer to have it talk to one person. You should enter the person's name into a file with thousands of other names and let the computer process the data for everyone. Since you're doing this for the class and it's your own computer, let's assume you're a benign user who can reply with their name when asked for it.

Try this as greeting_3.py

#! /usr/bin/python 
""" greeting_3 
    greeting as a function
"""

#------
#functions

def print_greeting(user_name):
	"""prints greeting
	"""
	print "hello %s, nice to meet you" % user_name

# print_greeting---------

#main()

#get input
resp = raw_input ("What is your name? ")

#call the function
print_greeting(resp)

# greeting_3.py ---------------------------------

In main() the parameter passed is resp (a string). In print_greeting(), the parameter arrives as user_name. What's going on? In main() the instruction print_greeting(resp) pushes the value of resp onto the parameter stack; the instruction doesn't push the name resp onto the parameter stack. When the function print_greeting() runs, the command print_greeting(user_name) says "assign the value of the first parameter to the variable user_name". At this stage the function doesn't know whether the parameter is a string, int or real (the data in the parameter can only be reasonably treated as a string). You have to tell the function that the parameter is a string, which you do in the first instruction that uses user_name (the print statement).

In greeting_3.py

Scope is the range in which a variable can be referenced. In programming languages,

Let's look at this modified version of greeting_3.py. Call it greeting_4.py. It has extra print statements (to debug the code).

	Note
	Since adding extra print statements to code fills the screen with output, it's helpful to prepend the name of the routine/function to each output so you can decipher the extra output. When you're done debugging the code, you will first comment out and then later delete these extra lines.

#! /usr/bin/python 
""" greeting_4 
    greeting as a function
"""
#-----------------
#functions

def print_greeting(user_name):
	"""prints greeting
	"""
	print "hello %s, nice to meet you" % user_name
	print "print_greeting: the value of resp is " + resp

# print_greeting--------------

#main()

#get input
resp = raw_input ("What is your name? ")
#call the function
print_greeting(resp)
print "main: the value of user_name is " + user_name

# greeting_4.py ---------------------------------

Variables are only known/visible to functions at levels above in the calling hierachy. Variables are not visible to functions at the same level in the calling hierachy. (The calling level of functions is also described as nesting. Outside the function, variables are known only to nested layers above.)

Scope allows people to write functions independantly of each other knowing that only main() will be able to see the contents of their function. A variable user_name could be used in several functions, and each user_name could have different values.

Because of scope

While you don't have to know about variable names used by other function, you do have to know names in the code at levels above (which call your routine/function).

In greeting_4() there are two variables

The program has instructions to output the value of both variables in each of the two functions main() and print_greeting(). Predict the output of the running the program. Then run the code and explain the output from the two extra print statements ^[79] .

	Note
	End Lesson 10

The namespace of the code that starts execution is called "global namespace". When you pass a file to the python interpreter, python starts executing on the first instruction it finds in global namespace, i.e. the first code that isn't a function. In other languages, execution starts at a function called main(), which has the top level (global) name space and which calls all other functions. For these lessons, in each file which has function(s), I've put a comment #main() in the place where the python interpreter starts executing. In code thus far, both python and other languages would start executing at the same place. However in a later section on making a module, we will see how python's assumption, about where to start executing, gets us into trouble.

Variables are declared in a namespace: if you declare a variable user_name in a function print_greeting(), then user_name exists in the print_greeting() namespace.

Any function can change the values of a global variable (including functions that haven't been written yet). If none of your functions change any of the global variables, you can't guarantee that someone else, at some later time, won't write one that does. Or you could blunder and write one that changes a global variable without realising it.

Note

In some other languages, any function can read and write global variables; in others you declare a global variable to be read-only by declaring it to be const (constant). In python, global variables are read-only. If you want to write to a global variable (i.e. change it), you have to explicitely declare that you're going to do so inside the function. This "are you sure?" step, requires you to think twice about doing it, but doesn't stop you from doing it. (Thanks to Rick Zantow for his explanation at global variables http://www.velocityreviews.com/forums/t354870-unboundlocalerror-local-variable-colorindex-referenced.html)

Global variables are regarded as sign of poor (or unsafe) programming style. If you use a writeable global variable, expect someone to ask you to justify doing so and to justify not using some safer programming practice. In later demonstration code, you will see global variables being used: usually the safer programming practice would add complexity, obscuring the point of the demonstration code. If you're programming rocket guidance systems, air traffic control or heart monitors, you won't be using any global variables.

As well as requiring data from the calling routine (passed as parameters), functions are often written to generate a result that's returned to the calling routine.

Call this file volume_sphere.py.

#! /usr/bin/python
""" volume_sphere.py
    returns volume of sphere
"""

#---------------
#functions

def volume_sphere(r):
	pi=3.14159
	result=(4/3)*pi*r*r*r
	return result

# volume_sphere---------------

#main()

#get input
resp = raw_input ("What is the radius of your sphere? ")

#call the function
volume=volume_sphere(resp)	#assign the result of the function call to volume

#formatted output
print "the volume of your sphere is %f" % volume
#this gives the same output
#print "the volume of your sphere is " + repr(volume)

# volume_sphere.py ---------------------

Let's walk through this line by line

#get input
resp = raw_input ("What is the radius of your sphere? ")

You'll be prompted for a number, which you'll enter on the keyboard. Your entry will be assigned to resp

#call the function
volume=volume_sphere(resp)	#assign the result of the function call to volume

When there's several things to do on a line, the parser starts from the right hand side. The parser has a value for resp so next looks at volume_sphere() which, because of the (), is recognised as a call to a function. The interpreter finds volume_sphere and passes resp to volume_sphere(). Note: the interpreter doesn't look at volume=, at least yet. Execution now passes to volume_sphere()

	pi=3.14159

the variable π is assigned the value 3.14159.

	result=(4/3)*pi*r*r*r

the value of (4/3 pi*r³) is assigned to result. Why didn't we just splice the value of π into the formula for the volume of the sphere?

(We'll shortly replace the line pi=3.14159; stay tuned.)

Here's the new part for this lesson.

	return result

The function pushes result onto the stack and exits. (You can only return one item. If you want to return two numbers, you could return one list containing two numbers.) Execution returns to main() where the interpreter now sees

	volume=result

volume is assigned the value of result

Run the program. Did you get a run-time error about non-ints? An error at run-time means that the code syntax was correct and the program started to run, but you asked the program to do something it can't do. Writing code is full of situations like this; it's a coder's life.. You have to figure your way out of them all the time.

The error message you get here is misleading. The problem is with resp; is it a number? If not, what is it ^[80] ? You sent the volume_sphere() a string and asked it to use a string as a number, which it can't do, so the program crashed.

The problem with resp is fixed by turning the string contained in resp into a real (floating point) number. All languages have ways of interconverting numbers and strings and they all look much the same. Here's how it's done in python.

number=float(string)

Note

Strong Typing and Weak Typing. Strongly typed languages will not allow you to send incompatible data types to a function. They check data types before the code is ever allowed to run. Strongly typed languages require the declaration of the function to include the type of any result and the type of each parameter. The above code in a strongly typed language would fail at compile/interpret time and you'd get an error message about mismatched parameter types ("resp" is a string and can't be passed to "volume_sphere" which requires a real). The development time for a strongly typed language is longer, but there will be less chance of errors at run time. Python is a weakly typed language. Weak typing allows you to quickly write short pieces of code without the relatively large overhead of writing header files and explicitely typing (declaring) parameters. For small programs in a strongly typed language like C or C++, the overhead of writing header files can be daunting, Weak typing also allows you to write sloppy code, which passes syntax checking, but which crashes at run-time. So you don't use python for rocket guidance, air traffic control or heart monitors. You do use it for short programs, for programs that will only be run a few times and for programs where errors will not be catastrophic and where someone will be on hand to fix them when they do cause a crash. Both strongly typed languages with a long development time and low error rate, and weakly typed languages with a short development time and high error rate have their place. A bank chooses a bullet proof armoured wagon to transport its money and you choose a paper bag to transport your groceries.

You can fix the code however you want. Here's my fixed code.

#! /usr/bin/python
""" volume_sphere.py
    returns volume of sphere
"""

#---------------
#functions

def volume_sphere(r):
	r=float(r)      #allow function to accept both strings and numbers
	pi=3.14159
	result=(4/3)*pi*r*r*r
	return result

# volume_sphere---------------

#main()

#get input
resp = raw_input ("What is the radius of your sphere? ")
#resp=float(resp)       #don't turn resp from a string into a float here

#call the function, passing a string
volume=volume_sphere(resp)

#formatted output
print "the volume of your sphere is %f" % volume

# volume_sphere.py ---------------------

Run this program with a radius of 1 (which gives me the result 3.141590). Now you have to check that your result is OK. Try with bc.

echo "(4/3)*3.14159" | bc
3.14159

It's the same answer as the python program. Does the answer look OK? If you're not sure, do it by hand (give π a value of 3 just to get a rough estimate) ^[83] .

So what went wrong? In case you've forgotten: (4/3) is an integer operation which gives a result of 1. The volume of a sphere of radius 1 is 4.188 (and not 3.14159 as found by the program). To invoke real number operations with bc, you need the math libraries, invoked with bc -l. To stop yourself falling into a trap, always invoke bc with the -l option.

echo "(4/3)*3.14159" | bc -l
4.18878666666666666665
echo "(4.0/3.0)*3.14159" | bc -l
4.18878666666666666665

Here's the fix

	result=(4.0/3.0)*pi*r*r*r

It's not a good idea to type in the value of π you could make a mistake. Any language that does math, has the value for π (and other useful numbers) built in. You don't need to remember exactly how to include π. Instead look up the documentation or look in google for a piece of working code (google for "python math pi") - I found CGNS Python modules -- Introduction Python" (http://cgns-python.berlios.de/Ipython.html).

There's a couple of places you can put the import statement.

All of these will work. However functions are designed to be easy to move/copy to other code. What will happen to your function if you put the import math statement at the beginning of the file or at the beginning of main() and you copy the function to another python file, by swiping it with your mouse ^[84] ?

Modify your code to use the value of π from the python math module.

Here's what good documentation for a function looks like. If ever you hope to code with other people, or to remember in 6 months what your code does, you will need to do something like this for every function/piece of code you write.

#! /usr/bin/python
""" volume_sphere.py
    returns volume of sphere
"""

#---------------
#functions

def volume_sphere(r):
	"""name      - volume_sphere(r)
	   version   - 1.0 Feb 2008
	   method    - volume =(4.0/3.0)*pi*radius^3
	   parameter - radius as string, integer or real, valid all +/-ve numbers
	   returns   - volume of sphere, float
	   Author    - Homer Simpson, Homer@simpson.com (C) 2008
	   License   - GPL v3.
	"""

	import math     #import the whole math module (including lots of stuff you don't need)
	r=float(r)      #convert string input to float.
	                #this allows the function to accept either a string or a number

	                #pi in module math is called math.pi
	result=(4.0/3.0)*math.pi*r*r*r
	return result

# volume_sphere---------------

#main()

#get input
resp = raw_input ("What is the radius of your sphere? ")

#call the function, passing a string
volume=volume_sphere(resp)

#formatted output
print "the volume of your sphere is %f" % volume

# volume_sphere.py ---------------------

Add appropriate documentation to your version of this function showing

It's easy to disguise where a function returns. Functions can have many conditional statements, any of which (depending on the flow of the program) can give the return value. Code is allowed to return from anywhere in the function, but doing this makes it hard to read. Badly written code looks like this.

def my_function (param1, param2...)
	if (some complicated condition)
		.
		.
		return (some complicated expression)
	elif (some complicated condition)
		.
		.
		return (some complicated expression)
	elif (some complicated condition)
		.
		.
		return (some complicated expression)
	else
		return (some complicated expression)
	endif

It can be hard to debug this sort of code: it's hard to find all the return statements. Instead there must (should) only be one return statement; the last statement in the function. All output should be put in some obviously named variable (like "result") which will be the argument to the return statement. The code you should write is this

def my_function (param1, param2...)
	if (some complicated condition)
		.
		.
		result = (some complicated expression)
	elif (some complicated condition)
		.
		.
		result = (some complicated expression)
	elif (some complicated condition)
		.
		.
		result = (some complicated expression)
	else
		result = (some complicated expression)
	endif

	return result

This is still hard to read, but you've done your best (multiple conditional statements aren't easy to read). If you're looking for the return values, you can find the string "result" with your editor,

procedure/subroutine/functions have the following properties

	Note
	End Lesson 11

Let's say your function is one or both of

You make a module out of your code. From there, you can import the function into your other python code.

In its most basic form, a python module is just a file containing a function; making a module only a matter of copying the function into a file in a directory in PYTHONPATH (which includes your current directory). However you can do a much nicer job than that.

Before you unleash this code on the unsuspecting world: is it safe? Once people start copying it to their own machine, you've lost control of it and you can't fix it if it's broken. At that stage, you'll have to endure e-mails from irate users telling you that your function is crashing their machine, or their rockets are blowing up.

Before we let the code go, we want to to exercise all the features of volume_sphere.py (it can accept string, reals, +ve/-ve). Add these lines at the bottom of the code and rerun it to check your file.

#make a few calls to volume_sphere() to test it
print "the volume of a sphere of radius  1 is " + repr(volume_sphere(1))
print "the volume of a sphere of radius -1 is " + repr(volume_sphere("-1"))
print "the volume of a sphere of radius  4 is " + repr(volume_sphere("4"))

For the moment, we'll put the module in your class_files directory (so we don't have to change PYTHONPATH). Let's first make sure it can be called as a module. Make a copy of volume_sphere.py as volume.py. volume.py will be your module file (it will eventually hold all your functions that calculate volumes of solids; for the moment it only has the function volume_sphere). Then at the command line, do this:

# python -i volume.py
What is the radius of your sphere? 6
the volume of your sphere is 904.778684
the volume of a sphere of radius  1 is 4.1887902047863905
the volume of a sphere of radius -1 is -4.1887902047863905
the volume of a sphere of radius  4 is 268.08257310632899
>>> volume_sphere(0)	#run your own commands from the prompt
0.0

You told python to execute volume.py, which it did, asking you for a radius and returning a volume. Then python ran through the added tests and stopped (python had reached the end of volume.py) but python didn't exit; it stayed in immediate mode (the -i option; to see the difference run the same command without the -i option). Having loaded your file, python now knows about the new function volume_sphere() and you can run any extra tests e.g. volume_sphere(0).

Create and run this file (call_volume.py). It mimicks a piece of code that calls your module.

#! /usr/bin/python
from volume import volume_sphere

print
print "Testing code from call_volume.py"
print "the volume of a sphere of radius  2 is " + repr(volume_sphere(2))
print "the volume of a sphere of radius -2 is " + repr(volume_sphere("-2"))
print "the volume of a sphere of radius  0 is " + repr(volume_sphere("0"))

The new piece of code is the line

from volume import volume_sphere

which says "from the module file volume.py (or volume.pyc) import the function volume_sphere().

Note

The first time python imports your module file, it will make a bytecode version with a pyc extension, which it will subsequently use. Look for a new file volume.pyc in your class_files directory. bytecode: a platform independant binary version of your .py file. It loads slightly faster, but runs at the same speed. You hand around the pyc file, if you want to make it difficult for the user to figure out the source code. (This is antithetical to the idea of GNU computing and the idea of the internet being a place to freely exchange information. A pyc file enables someone to make money by depriving you of information.) How would you check that the py wasn't being used [85] ? How would you check that the pyc file was being used (in the absence of a py rather than the py file [86] ?

The output is

# ./call_volume.py
What is the radius of your sphere? 7
the volume of your sphere is 1436.755040
the volume of a sphere of radius  1 is 4.1887902047863905
the volume of a sphere of radius -1 is -4.1887902047863905
the volume of a sphere of radius  4 is 268.08257310632899

Testing code from call_volume.py
the volume of a sphere of radius  2 is 33.510321638291124
the volume of a sphere of radius -2 is -33.510321638291124
the volume of a sphere of radius  0 is 0.0

You asked python to load the function volume_sphere and then execute the code in the calling routine (print statements which call the function volume_sphere). You see the requested output in the 2nd block above (starting "Testing code ...").

The output starting "What is the..." is from code in the module file's global namespace. You did not ask python to execute this code, but it did anyway. Every module needs built-in testing code, but you don't want it executed every time you load the module. Handling this is the next step in building a module.

In a normal language, code starts and ends execution in main(). Python is different: it starts executing with the first code it finds in global namespace. You would like your program to start execution in the same place, no matter what order you load your files, but python will start executing in a different place if you reorder your files. This is not what you want.

You could solve this problem by commenting out the global namespace code in your module, but then you'd have to write a separate file to test your module. The chances of keeping your module file and your module testing files together for decades is small. Your module has to be able to test itself and not rely on external code for testing.

The Python 2.5 Documentation (http://docs.python.org/download.html) in section 6.1 (More on Modules) says that the global namespace code in executed once on loading to allow initialisation of the function(s).

	Note
	compiled languages have a linker/loader to handle much of this

Note

What/why you'd do initialisation at load time (from a posting to the TriLUG mailing list - the poster didn't want credit, he said it was all common knowledge): The most common case is a module will load more modules when loaded. That is, you import a module and the first lines in the module you load will be import statements which load yet more modules (and in turn, those modules may import still more modules). The module might selectively load certain other modules in a plug-in architecture based on things the module can introspect from the environment into which it was loaded. You might want to cache some data from somewhere. If you are doing metaprogramming in the module, you might want to create some new objects of the metatype you just loaded. If you've defined some factory methods in the module, you might want to fire out some objects. I think there's no end of uses for executing global namespace code when loading a module. I can think of no case where a module would be useful without executing code in its global namespace. Basically, all the statements in a Python module are executed in order when loaded, whether as a script or as an import. The execution of most statements in a module, though, just define code objects (i.e., function objects, class objects, etc.) which will be used later (in the script, or by another script or module importing the module). So you might think of that as function and class objects getting "initialized" with their byte code.

There is a fix for code you don't want run from global namespace (it makes python execute in the same way that other languages execute their code).

First here's an normal module with its test code. You'll execute ./module.py.

#/usr/bin/python
#this file is called module.py
def function():
	#code for the function

#global namespace
#main() for this executable
#module initialisation code (if you need it)
#	eg if your function generated random numbers, 
#	you would seed (initialise) the random number generator (with say the date)
#	do that the random number generator would give different random numbers 
#	each time your ran the module

function(1)	#this is your testing code. You only want this to run when you execute the module by itself.

Python loads the function and then coming to the global namespace starts executing. What there is to execute is the line of test code function(1). You'll see the output from calling function(1) (whatever that may be).

output from function(1)

Next here's a bit of code that calls the function.

#!/usr/bin/python
#this file is called call_module.py
from module import function

#global namespace
#main() for this executable
function(2)

When you run the command call_module.py, function() is imported (loaded into memory, but not executed). As well (the bit you don't want) python starts executing in the first global namespace it finds, which is the line function(1) in the module. Next python will start executing in the global namespace for call_module.py, which is the call function(2). What you'll see will be

output from function(1) 	#you didn't ask for this (it's your testing code)
output from function(2)		#this is what you wanted

However the global namespace for module.py is not main() for call_module.py (the piece of code that's making all the calls). Once the code is in memory and just before it starts to run, here's what it looks like.

#module
def function():
	#code

#global namespace
#module initialisation code (if you need it)
#this is NOT main() anymore
function(1)
.
.
.

from module import function

#global namespace
#main() for this executable
function(2)

Since python runs global namespace code (whether it's in main() or not) whereever it finds it, you'll see the output from calling function(1) then function(2).

You apply the fix, a conditional, to the module file. It says "if this isn't main(), don't run it". Here's the fix.

#module
def function():
	#code

#global namespace
#module initialisation code (if you need it)
if __name__ == '__main__':
	#this is NOT main() 
	function(1)
.
.
.

from module import function

#global namespace
#main() for this executable
function(2)

Here's the output, with python now behaving like other languages.

#module initialisation code (if you need it)
output from function(2)		#this is what you wanted

The fix is a conditional statement

if __name__ == '__main__':

which says "If code execution starts here, then execute these indented statements". For the moment, you can regard syntax as magic (I do). The testing code is now governed by the conditional statement. If you execute the module by itself, the first global namespace code python will see will be in your module, it will also be main() for this executable, python will call this code "__main__" and the conditional will be executed. If code execution starts in global namespace code in some other file, (e.g. some routine calling volume_sphere.py) then the testing code in the module will not be called "__main__" and will not be run.

Here's the fixed version of volume.py containing the function volume_sphere.py.

#---------------
#! /usr/bin/python
#functions
def volume_sphere(r):
	"""name      - volume_sphere(r)
	   version   - 1.0 Feb 2008
	   method    - volume =(4.0/3.0)*pi*radius^3
	   parameter - radius as string, integer or real, valid all +/-ve numbers
	   returns   - volume of sphere, float
	   Author    - Homer Simpson, Homer@simpson.com (C) 2008
	   License   - GPL v3.
	"""

	import math     #import the whole math module (including lots of stuff you don't need)
	r=float(r)      #convert string input to float.
			#this allows the function to accept either a string or a number

			#pi in module math is called math.pi
	result=(4.0/3.0)*math.pi*r*r*r
	return result
#volume_sphere

#---------------
#main()

print "initialising volume.py"
print

if __name__ == '__main__':
	print "self tests for volume.py"
	print

	print "self tests for function volume_sphere()"

	##turn off the interactive test(s)
	##get input
	#resp = raw_input ("What is the radius of your sphere? ")
	##call the function, passing a string
	#volume=volume_sphere(resp)
	##formatted output
	#print "the volume of your sphere is %f" % volume

	#make a few calls to volume_sphere() to test it
	print "the volume of a sphere of radius  1 is " + repr(volume_sphere(1))
	print "the volume of a sphere of radius -1 is " + repr(volume_sphere("-1"))
	print "the volume of a sphere of radius  4 is " + repr(volume_sphere("4"))

# volume.py ---------------------

I've added code in the place where you would initialise the module on loading ("self tests...).

Here's the module being run in self testing mode (note a string indicating initialisation, a string indicating that file volume.py is being run, and a string indicating that tests are being run for module volume_sphere).

# ./volume.py
initialising volume.py

self tests for volume.py

self tests for volume_sphere()
What is the radius of your sphere? 7
the volume of your sphere is 1436.755040
the volume of a sphere of radius  1 is 4.1887902047863905
the volume of a sphere of radius -1 is -4.1887902047863905
the volume of a sphere of radius  4 is 268.08257310632899

Here's the module being run in interpreted mode. You can run your own tests at the end (here we outputted the volume for a sphere of radius "0").

# python -i ./volume.py
python -i ./volume.py
initialising volume.py

self tests for volume.py

self tests for function volume_sphere()
What is the radius of your sphere? 7
the volume of your sphere is 1436.755040
the volume of a sphere of radius  1 is 4.1887902047863905
the volume of a sphere of radius -1 is -4.1887902047863905
the volume of a sphere of radius  4 is 268.08257310632899
>>> volume_sphere("0")
0.0

Here's the module being called from another piece of code. Note that the initialisation code (before the "if" statement) is still run, but that now the self testing code is not run.

# ./call_volume.py 
initialising volume.py


Testing code from call_volume_sphere.py
the volume of a sphere of radius  2 is 33.510321638291124
the volume of a sphere of radius -2 is -33.510321638291124
the volume of a sphere of radius  0 is 0.0

Summary: make the module self testing. That way a developer, unsure of a problem they're having can quickly re-run self tests on suspect modules. (People do accidentally change the wrong files, or "upgrade" them.)

By now I hope you realise

Much code is short programs that are used a few times and then thrown away. People don't document them well (or at all) and usually don't get into too much trouble for the lack of documentation.

Other code will be upgraded, or used unmodified for decades. The original author may be long gone, while the code is run hundreds of times a day, or once a year.

It turns out that the major cost in writing software, is not the initial cost of producing the first working version or even the production version, but maintaining the code (over decades). Someone who's never seen the code, doesn't know what it does or how it works, will be assigned the job of fixing/updating your code. It may take months before they can make the first changes. This may be longer than it took the original author to write the code. The prime aim of the coder should be to write code that a human can read.

Because code is unmaintainable, it's often simpler to write a new version from scratch. This is a terrible waste of resources. Managers, who often know nothing about writing code, are willing collaborators in this tradgedy. They will be sucked in by the arrival of a new language and say "this code was written in (some old language), it's time for us to be running (this year's fancy new language)". The effort to re-write the code from scratch in a new language is seen as a development cost. It's not, it's maintenance and is unneccessary maintenance at that. The old code was already debugged, an enormous investment; you can't throw away tested and debugged code. The computer doesn't care about the language the code was written in. Adding two numbers is the same to the computer, whether it was written in Fortran or Java. The new version will often be functionally worse than the original (and no more maintainable), but management will not realise this. To justify the money and time spent on the rewrite, management will tout the code rewritten in this year's new language, as an advance.

You code must be documented, and easy to read. Because reading other people's code is hard, anyone can write compact, fiendishlessly smart and cryptic code, that no-one else can read and that even the author won't understand 6 months later. Don't be too smart: make your code simple. Your first objective is to write maintainable code. Your prime objective in writing code, is not that the code works (which it must), or that it works correctly (which it must) or that it's fast (nice, but not required), but that it's maintainable.

Most GPL packages now have testing code which is run before installation.

Once your managers realise you have working code, they can claim their money from the customer, and they'll want you to move onto something else (which will make them money). You'll get little support for writing documentation and testing routines. When someone says to you "it works? then we don't need testing routines" you're talking to someone who doesn't know how to code and doesn't understand the enormous cost of maintaining bad code. You may have to do what they say, but you don't have to agree with them.

In the module volume, one of the functions does an import math to get the value of π. If you have multiple functions each import'ing math, you could move this instruction to the module global namespace initialisation area, so math only has to be imported once. Your module would work fine, but you shouldn't do this. Why ^[87] ?

	Note
	End Lesson 12

Functions are an essential part of a program. Functions look much the same in all languages and are invoked the same way (by passing parameters/arguments and accepting a returned variable). Most action in a program will be in functions. Any self contained logical block in a program is a candidate for being made into a function. Functions can be called from anywhere (i.e. from main() or from other functions).

The calling routine passes parameters; the called routine is passed parameters. The name of the parameters in the calling routine is usually different to those used in the called routine. The use of different names is for readability. The name of the parameter passed from the calling routine is usually specific to the calling routine. The name used for a parameter in the function is generic, a name that will work no matter what it's called by. If the same name is used for a parameter in the calling function and in the called function, then there are two variables, each in their own scope, with the same name: changes to the variable in the called function do not affect the variable of the same name in the calling function. The list of parameters is often called the interface - it's the exchange of information that must occur for a function to work.

#! /usr/bin/python

#my_program.py

#functions

def my_function(param1, param2..paramn):
	#calculate result
	return result

#main()
	#get data
	#do calculations
	#sample function call
	my_function(variable1, variable2..variablen)
	#output results

main() should be short, with as much work as possible handed off to functions. Functions can be called at any step in main(). Often the control of the flow of the program is complicated, in which case main() will be longer.

Code is turned into a function for

Functions must have documentation, enough for someone reading the docs to be able to recreate your code without having seen the code (for an example see documentation.

Continuing our exploration of solid objects, write a module to calculate the volume of a hexagonal prism.

A prism Prism (Geometry) (http://en.wikipedia.org/wiki/Prism_%28geometry%29) is an n-sided polygon translated through space. A cube (and a rectangular block) is a prism, because (at least one) of its faces is the same as you slice off pieces parallel to the face. A cylinder is not a prism because the face is a circle, not a polygon. An ice crystal is a prism: it has a regular hexagon base and extends at right angles to the base. Many crystals are prisms.

Hexagonal ice prisms in thin layers of cirrus clouds (in NC seen mainly in winter) are responsible for spectacular ice rainbows and halos e.g. Frequent Halos (http://www.atoptics.co.uk/halo/common.htm), The Cloud Appreciation Society (http://cloudappreciationsociety.org/version2/wp-content/uploads/2008/02/08-feb-high1.jpg) and sundogs seen in a circle of radius 22° from the sun (see 22° halo http://en.wikipedia.org/wiki/22%C2%B0_halo).

The formula for the area of a regular hexagon (in algebra) from Areas and Perimeters of Regular Polygons (http://www.algebralab.org/lessons/lesson.aspx?file=Geometry_AreaPerimeterRegularPolygons.xml) and Area of a Regular Polygon (http://www.mathwords.com/a/area_regular_polygon.htm) is

Area = ((3*sqrt(3)/8)*d^2	#d=vertex to vertex (diameter of circumcircle)
Area = ((3*sqrt(3)/2)*r^2	#r=center to vertex (radius of circumcircle) 
Area = 2*sqrt(3)*apothem^2	#apotherm=cente to face (radius of inscribed circle)
Area = (2*sqrt(3)/4)*face-to-face-diam^2

	Note
	In python the square root can be calculated a couple of ways >>> number=2 >>> number**(1.0/2.0) 1.4142135623730951 >>> number=9 >>> from math import sqrt >>> sqrt(number) 3.0

Write a self contained file volume_hexagonal_prism.py (later you will then turn the file into a module and add it, along with self tests, into volume.py). Here are the specifications:

	Note
	The simplest thing to do would be to copy volume_sphere.py to volume_hexagonal_prism.py and start hacking on that code. You never write code from scratch if you have something similar, that's debugged and tested.

Check your output using a few test inputs. Make the test code appropriate for the volume_hexagonal_prism() case. Finish by commenting out the tests requiring user input and check that the non-interactive tests run OK. Check that the output is equivalent to running the volume_sphere() code in volume.py.

Remember the formula I gave you for the volume of a sphere that contained "(4/3)" which the computer evaluated as an integer, and the formula gave the wrong result? Always check your answer by hand. You might think that you could plug the numbers into a 4 function calculator and get the right answer and in this case you likely will. However one day you'll make a mistake without knowing it. The people who built the Hubble Space telescope, did a complicated test on the mirror surface and got the wrong answer. An amateur telescope maker using a knife edge and a flashlight (torch) bulb (the Foucault Test, http://en.wikipedia.org/wiki/Amateur_telescope_making#Foucault_test) would have detected the problem in a matter of minutes. The simplest test wins.

For a complicated shape like a hexagon, you have to think a bit to get a simple calculation that gives a good enough answer. First check the area of a hexagon. Here's two checks

You now know that the area of a hexagon of vertex-to-vertex diameter=1 is between 0.5 and 0.78.

Check the volume of a hexagonal prism with height=1. Then try combinations of height/diameter of 0,1. Check the change in volume if you double the height, double the diameter (by what ratio should the volume change for a doubling of the base diameter, a doubling in the height?).

Here's my code for volume_hexagonal_prism.py together with the output from three pairs of checks run from the command line. ^[92]

The print statement in main() does not have repr() for diameter and height, but does have it for volume. Why ^[93] ?

	Note
	End Lesson 13

You're likely to write lots of functions. You put them in modules with similar functions. Here we're going to put the function volume_hexagonal_prism() into the module volume We've already put the function volume_sphere() into the module volume. Let's look at what happens when we add a 2nd function to the module.

The simplest scheme is to copy/paste volume_hexagonal_prism.py into volume.py this way.

	Note
	This code is called pseudo code - it's half English, half computer code; it's sufficiently generic that it could be implemented in any computer language.

shebang (if needed)
#scheme 1, for small modules

	"""
	collection of functions that calculate the volume of solids
	Authors attributed in each function.
	"""

volume_sphere()
	#code

volume_hexagonal_prism()
	#code

#----------
#main()

global namespace initialisation code

if __name__ == '__main__': 

	tests on volume_sphere		#multiple lines
	.
	.

	tests on volume_hexagonal_prism	#multiple lines
	.
	.

#--module volume

If you do something once, you can do it almost anyway you want. As soon as you do something twice, (like having two functions in a module) you have to consider what will happen if you do it 100 times. This is the problem of scaling mentioned in order of algorithm.

Inserting the functions is easy - you just stack them in order from the top of the file. It's what to do with the tests that's the problem. If you have 10 functions, each with 10 lines of tests (and comments) which you move into the module's global namespace (at the bottom of the file), you will have 100 lines of tests. This is too long a block of code to read and now the tests are a long way from the code they're testing. You'll have to do some thinking if you ever want to modify the file (like to turn off some of the tests).

On looking at your code with two functions, programmers will first ask "but does it scale?". Most often there aren't elegant solutions; mostly there are solutions that people will accept (or put up with). More often than you would like, the best solution anyone can think of is plain downright ugly (that's life).

You have some latitude on what to do here. Your main idea is for someone else to be glad to read your code.

One of the principles of code writing is that when a function (or main()) becomes too big, you split out a self contained logical block into a function. To do that here, put the tests for each function into their own test_function(), and call the test_function()s.

shebang (if needed)
#scheme 2, for medium sized modules

	"""
	collection of functions that calculate the volume of solids
	Authors attributed in each function.
	"""

volume_sphere()
	#code

test_volume_sphere()
	#tests

volume_hexagonal_prism()
	#code

test_volume_hexagonal_prism()
	#tests

#----------
#main()

global namespace initialisation code

if __name__ == '__main__': 

	test_volume_sphere()		#one line
	test_volume_hexagonal_prism()	#one line

This will work for 100 functions - you'd have 100 lines of test_function_name() at the bottom of the file. People would realise that this was a big module file and that the calls at the bottom are just a list of calls to test functions. Once the number of lines of calls in main() becomes unmanageable, you collect them into a module.

shebang (if needed)
#scheme 3, for huge modules

	"""
	collection of functions that calculate the volume of solids
	Authors attributed in each function.
	"""
volume_sphere()
	#code

test_volume_sphere()
	#tests

volume_hexagonal_prism()
	#code

test_volume_hexagonal_prism()
	#tests

run_all_tests():
	test_volume_sphere()		#one line
	test_volume_hexagonal_prism()	#one line
	.
	.

#----------
#main()

global namespace initialisation code

if __name__ == '__main__': 

	run_all_tests()

If you've got 100 functions, this 3rd way would be the best.

How do you know whether you have a small, medium or huge module? If you (or the users) can't understand the mess of code in main(); if they can't, you need to go to next scheme.

For the moment, let's make the module using Scheme 2.

Here's the specificiations

I've commented out the interactive tests, so the tests will run without keyboard input. Here's my version ^[94] and here's the output

# ./volume.py
initialising volume.py

self tests for volume.py

self tests for function volume_sphere()
the volume of a sphere of radius  1 is 4.1887902047863905
the volume of a sphere of radius -1 is -4.1887902047863905
the volume of a sphere of radius  4 is 268.08257310632899

self tests for function volume_hexagonal_prism()
the volume of your hexagonal_prism of diameter 1.0 and height 1.0 is 0.649519052838329
the volume of your hexagonal_prism of diameter 2.0 and height 1.0 is 2.598076211353316
the volume of your hexagonal_prism of diameter 1.0 and height 2.0 is 1.299038105676658

	Note
	End Lesson 14

A person installing your module and running the above tests has no idea if the output is correct. We have to compare the answer found with the expected answer, then output a pass/fail message. The user doesn't have to see the actual numbers (although they can).

Comparing reals is a problem (as we'll see in don't equate reals). The number on the screen (or in your code) may differ by a couple of low order bits from the number in memory. The number "0.1" will actually be represented by "0.10000000000000001". You have to compare the ratio or find the difference. The ratio (division) test will fail if the divisor is "0.0", so it's probably best to test the difference between the expected and found numbers. Reals in python are 64-bit giving a precision of 2^-52. The numbers output by the test will either be as correct as the computer can make them, or will be obviously wrong.

Let's say to be safe that we'd like the difference between the numbers to be 2^-50 (giving us 2 bits margin of error in detecting whether there is an error). What's 2^-50 in decimal ^[95] ?

Note

You all know what 232 is [96] . Here's some more useful numbers 2^10 = (approx) 10^3 = 1k 2^20 = (approx) 10^6 = 1M 2^30 = (approx) 10^9 = 1G 2^40 = (approx) 10^12 = 1T 2^50 = (approx) 10^15 = 1P 2^-10 = (approx) 10^-3 = 1m 2^-20 = (approx) 10^-6 = 1u 2^-30 = (approx) 10^-9 = 1n 2^-40 = (approx) 10^-12 = 1p 2^-50 = (approx) 10^-15 = 1f Hard disk manufacturers have been successfully sued for labelling disks as having (for example) 1TByte of storage, which people expect to mean 1099511627776 (1.1*1012) bytes, when they have only 1012 bytes (10%less). Attempts to diambiguate 1024 and 1000 as a kilobyte have lead to an alternate nomenclature, which is unfortunately is clumsy (Mebibyte http://en.wikipedia.org/wiki/Mebibyte).

	Note
	For the real way to compare numbers see Lahey Floating Point Arithmetic (see the section on "Safe Comparisons") (http://www.lahey.com/float.htm).

The test requires the numbers to be the same at an accuracy of 10^-15, which in turn requires the output from your code to display 16 significant figures after the decimal point. If your output has only 12 significant figures (the default for the python print statement which uses str() which in turn prints 12 significant digits, see Floating Point Arithmetic, http://docs.python.org/tut/node16.html), the test will erroneously return fail. We will explore this more in Real Numbers. For the moment note

>>> j=0.12345678901234567890	#j has 20 digits after the decimal point
>>> print j
0.123456789012			#standard python printing has 12 digits
>>> j
0.12345678901234568		#64 bit real has about 16 significant digits
>>> print "%5.20f" %j		
0.12345678901234567737		#64-bit real printing 20 digits and showing garbage after 16 digits
>>> 

Check that your tests in volume.py give output with 16 figures after the decimal point before proceeding.

Start a file test_compare_results.py with a function compare_results() which has these specifications

Here's my code ^[98] and here's my output ^[99] . Did one of the last pair pass? If so look at the note in the 2nd bullet above.

Why did I test the pairs of numbers, twice, the 2nd time in the reverse order ^[100] ? Why did the first pair pass, while the second pair failed ^[101] ?

The parameters passed to compare_results() are reals. How would you change the function to handle parameters than were string representations of numbers ^[102] ? (You don't need to handle this here. The numbers you're comparing are generated by the computer and will be numbers and not strings.)

	Note
	End Lesson 15

Shortly we will put compare_results() into volume.py, but first we can test it interactively (before we risk messing up volume.py by changing the code). Here's the interactive code

# python -i test_compare_results.py
pass
pass
fail
fail
>>> from volume import volume_sphere
initialising volume.py

>>> print "the volume of a sphere of radius  1 is " + repr(volume_sphere(1))
the volume of a sphere of radius  1 is 4.1887902047863905

In the interactive code we

We want to modify this test to compare the output with the expected output. At the end of this last line, add a call to compare_results() to compare the output with the expected output, so as to print a pass/fail ^[105] .

Notice how the symbol "1" and "4.1887902047863905" are used multiple times and/or they are the arguments to expressions? Constants (numbers that aren't ever going to be changed by the program) should be assigned to variables and the variables used in expressions. The reasons for this are

Assign the constants to variables (e.g. radius, expected_result) and use the new variable, rather than a number, in the expressions. Here's my code ^[106] .

You have two calls to (e.g. volume_sphere(radius)) in this print line. You never calculate anything twice: instead assign the result to a variable and use the value of the variable twice. Assign the value of the volume to volume_found and rerun the test. Here's my code ^[107]

Here's the final version, which you developed using python interactively

# python -i test_compare_results.py
python -i test_compare_results.py
pass
pass
fail
fail
>>> from volume import volume_sphere
initialising volume.py

>>> radius=1
>>> expected_result=4.1887902047863905
>>> volume_found=volume_sphere(radius)
>>> print "the volume of a sphere of radius " + repr(radius) + " is " + repr(volume_found) + " " + compare_results(expected_result, volume_found)
the volume of a sphere of radius 1 is 4.1887902047863905 pass
>>> 

volume.py consists of the functions we're really interested in, volume_sphere() and volume_hexagonal_prism(), plus a matching pair of functions, which test the functions with various test inputs. However as yet we don't have any way to test if the output is correct. You now are going to merge the code from test_compare_results.py and the interactive code above, into volume.py, so that you can compare the output of the tests with the known answer.

Now do the merge

Here's my module volume_2.py ^[110] and here's the output ^[111] .

You are now safe to let modules out into the world. Congratulations.

	Note
	End Lesson 16

Note

I didn't use this material in class, as I couldn't see any teaching value in it. It's left here to show that sometimes there isn't a way to write good code. The same line of code print "the volume of a sphere of radius " + repr(radius) + \ " is " + repr(volume_found) + " " + compare_results(expected_result,volume_found) is used in multiple places in the tests, no matter what values are passed to the test. It's not to difficult for a reader to guess that the lines all might be the same (the text lines up from line to line), however it would take a bit of checking to be sure. If you wanted to change the line, you would have to do it multiple times (this is not good, from the point of view of maintenance, you might miss one line). It would be better to have only a single copy of the line of code. You could call it in a function, but a function with only one line of code, whose only purpose is to output a string for the user, is a bit pointless. Instead I rewrote the code to run the print statement inside a loop, using the loop to feed in the parameters. After looking at the resulting code, I couldn't see that functionally is was any different to calling a function. Since the normal idiom is to call a function, anyone trying to maintain the code would have to puzzle as to why it was written in a loop. I decided that for the number of times the tests would be run, that it really didn't matter if it was left the original way. Sometimes bad code is OK enough, or at least not worth the trouble of fixing.

We have several parameters (radius, expected_result) to feed to the line. However loops feed parameters one at a time e.g.

for item in items:
	#do something to item

Let's group the parameters into a list, which as far as a loop is concerned is a single object. Each iteration of the loop will feed a list to the loop code. Since we're doing multiple tests, we'll need to feed many lists, so let's have a list of lists.

Before doing anything with lists, look at this info about lists and list of lists.

Here parameter_list[] is a list of reals. loop_list[] is a list of parameter_list[]s (i.e. a list of lists). The code for using lists to feed parameters to our tests is

#construct lists from input data
loop_list = []	#initialise list

#first test
parameter_list = [ radius_1, expected_number_1, volume_1 ]
loop_list.append(parameter_list)		#add parameter_list to the tail of loop_list 

#second test
parameter_list = [ radius_2, expected_number_2, volume_2 ]
loop_list.append(parameter_list)		#add parameter_list to the tail of loop_list 

#we're reusing the variable parameter_list here. It's reinitialised in this statement
for parameter_list in loop_list:		#recover parameter_list one at a time
	volume = parameter_list.pop() 	#pop() removes from the tail, so last in is first out
	expected_number = parameter_list.pop() 
	radius = parameter_list.pop()
	print "the volume of a sphere of radius " + repr(radius) + " is " + repr(volume_sphere(radius)) + " " + compare_results(expected_result,volume_sphere(radius))

Modify your code to run your tests in a loop (copy volume_2.py to volume_3.py). Here's my code ^[112] .

To add (or modify) a test, it's just a matter of adding (or modifying) a paragraph like this

	diameter = 1.0
	height = 2.0
	expected_result=1.299038105676658
	volume=volume_hexagonal_prism(diameter, height)
	parameter_list = [diameter, height, expected_result, volume ]
	loop_list.append(parameter_list)

	Note
	I wouldn't say that this version of the tests with a for loop is particularly more readable or modifiable than the previous version. However it's an alternate way of handling the scalability problem.

Code in compare_results() looked a little like the following. Notice any difference in readability between these two functionally identical pieces of code?

expected_result=1.0
volume=call_to_function(param1,param2...paramn)
compare_results(expected_result, volume)

expected_result=1.0
compare_results(expected_result,call_to_function(param1,param2...paramn))

The second line of the last piece of code is called a train wreck. It's not too big a train wreck as far as train wrecks go, but it's still a train wreck.

It's very easy, once you've set up all your functions, to have lines of code like this, that have calls to functions to any depth you can imagine (the params themselves could be calls to functions). The code is compact and demonstrates recursive parsing of code.

	Note
	Recursive parsing is the ability to recognise an expression, which will evaluate to a number, as being equivalent to a number. Early languages (e.g. Fortran) could not do recursive parsing. If instead of a number, an expression which evaluated to a number was found, the compiler/interpreter would crash. This was really stone age computing.

Compact code makes a section of code look neater (there's less lines of code). This is all very nice, except that compact code is hard to read and no-one (including you in 6 months), will have a clue what it's doing. The next person won't look forward to fixing it when it stops working.

When writing code, make sure in your mind that you separate what the code does from how you implement it.

You should be able ask someone to reimpliment the function (change the how), changing the code, without changing the what. Anyone using/testing the function should not be able to tell that the code in the function has been changed underneath them.

Here's an example piece of code from compare_results()

	difference = n_expected - n_found

	if (difference < 0.0):
	        difference *= -1.0

If you were explaining this code, you'd say that the function of this piece of code is to calculate the magnitude of (size of) the difference between n_expected,n_found. Depending on your audience, you may have to explain why you want the magnitude of the difference, or you may not. You would not explain what each line of code was doing - doing so is explaining the implementation and not the function. (In some circumstances, you might want to explain the implementation, but assume you're showing the code to a programmer - in this case you'll be explaining the function.) When explaining the function, you allow for the possibility that someone else will reimplement your code. Next time you see the code it might look like

	difference = mag_diff(n_expected,n_found)

and your explanation would be the same.

It should be possible for someone to reimplement your code in a computer language you've never seen, using a character set you don't know (Chinese, Cyrillic, Arabic) and for you to come into a room full of people who know this character set and computer language and for you to give the same explanation of how the code works, even though you can't read the displayed code.

People may ask you to explain your implementation of a function, in which case you'll have to justify it, or listen to a better one from the audience. Part of the reason you give a talk is to find out from people who know more that you. You should accept such help with humility: it's far better to find out from your friends or co-workers, than your rivals or the purchasing public, who decided to stay away from your product in droves.

Quite often there's a dozen ways of implementing something, all of which are equivalent. Instead of

	if (difference < 0.0):
	        difference *= -1.0

you could have written

	if (difference < 0.0):
	        difference = -difference

If you've thought about the other ways of implementing this (and before you give a talk, you should have) you'll know your way is as good as any other and you can just say "it was the first method I thought of" and leave it at that.

You should plan on describing the function of the code. Only explain the implementation if someone requests it, or you've done something really nifty. Be carefull though - people are wary of really nifty code - and may regard it as code that will be unmaintainable when you move to another project.

Seminars have a standardised format. You can go anywhere in the world, to a seminar on any topic and the format will be exactly the same.

There a couple of rules

I'll give a demo of how you should do it and then you'll follow. If it seems strange for you to stand up and do the same thing as I've just done, or to follow a student who's followed me, then just pretend you're on the soccer field and you're being asked to practice shooting goals, or you're in music class and you all are being asked to play "Fur Elise" one after another for the teacher. I want to see you shooting goals until I'm sure that when I walk away, you can shoot a goal by yourself. You may have to explain volume.py to multiple groups of people over a long period of time, so learn to enjoy doing it. Just think how many time Mick Jagger has sung "Satisfaction". Assume that you'll be explaining volume.py that many times. You may think that your skill (here coding) is all that counts. It's not till you sell it, that you collect your money.

wiki, real numbers (http://en.wikipedia.org/wiki/Real_number), wiki, IEEE 854 floating point representation (http://en.wikipedia.org/wiki/IEEE_854), wiki, floating point (http://en.wikipedia.org/wiki/Floating_point). IEEE 754 Converter (http://www.h-schmidt.net/FloatApplet/IEEE754.html) a java applet to convert floating point numbers to 32 bit IEEE 754 representation.

The Python docs on floating point Floating Point Arithmetic (see http://docs.python.org/tut/node16.html).

A compiler writer's view of floating point Lahey Floating Point (http://www.lahey.com/float.htm).

The reasoning behind the IEEE-754 spec: What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg, published in the March, 1991 issue of Computing Surveys. Copyright 1991, Association for Computing Machinery, Inc. The paper is available at ACM portal (http://portal.acm.org/citation.cfm?id=103163).

Scientists and engineers are big users of computing and do calculations on non integer numbers (e.g. 0.12). These numbers are called "real" to differentiate them from "imaginary" (now called "complex") numbers. Real numbers measure values that change continuously (smoothly) e.g. temperature, pressure, altitude, length, speed.

	Note
	You will never have to reproduce any of the following information in normal coding, since all the routines have been written for you, but you will have to understand the limitations of floating point representation, or get nonsense or invalid results.

One possible representation of real numbers is fixed point e.g., where the decimal point is in its usual position. Fixed point isn't useful for small or large numbers in fixed width fields (e.g. 8 columns wide) as significant digits can be lost. A small number like 0.000000000008 will be represented at 0.000000 and 123456789.0 will be represented as (possibly) 12345678 or 456789.0. So floating point notation is used. Here's some floating point numbers and their 8-column representation

real                fixed point    floating point
123456789.0         12345678       1.2345E8  (ie 1.2345*10^8)
0.0000000012345678  0.000000       1.234E-9  (ie 1.2345*10-9)

A decimal example of a floating point number is +1.203*10². There is one digit before the decimal point which has a value from 1-9 (the base here being 10). Subsequent digits after the decimal point are divided by 10,100... etc. The term 1.023 is called the mantissa, while the term 2 is called the exponent. A mantissa and exponent are the two components of scientific notation for real numbers. Here's how we calculate the decimal value of +1.203E2

sign     + positive
mantissa 1.203 = 1/1 + 2/10 + 0/100 + 3/1000
exponent 10^2 = 100
+1.203*10^2 = +120.3

Most computers use 64, 80 or 128-bit floating point numbers (32-bit is rarely used for optimisation calculations, because of its limited precision, see floating point precision). Due to the finite number of digits used, the computer representation of real numbers as floating point numbers is only an approximation. We'll explore the reason for this and the consequences here.

One of the early problems with computers was how to represent floating point numbers. Hardware manufacturers all used different schemes (which they touted as a "marketing advantage"), whose main effect was to make it difficult to move numbers between machines (unless you did it in ASCII, or did it between hardware from the same manufacturer). The purpose of these non-standard schemes was to lock the customer into buying hardware from the same manufacturer, when it came time to upgrade equipment.

You'd think that manufacturers would at least get it right, but these different schemes all gave different answers. A C.S. professor in the early days of programmable calculators (1970's) would stuff his jacket pockets with calculators and go around manufacturers and show them their calculators giving different answers. He'd say "what answer do you want?" and give it to them. Eventually a standard evolved (IEEE 854 and 754) which guaranteed the same result on any machine and allowed transfer of numbers between machines in standard format (without intermediate conversion to ASCII).

Hop over to IEEE_854 (http://en.wikipedia.org/wiki/IEEE_854) and look at the bar diagrams of 32- and 64-bit representation of real numbers.

Worked examples on 32- or 64-bit numbers take a while and because of the large number of digits involved, the results aren't always clear. Instead I'll make up a similar scheme for 8-bit reals and work through examples of 8-bit floating point numbers. If we were working on an 8-bit (rather than 32- or 64-bit) computer, here's what an 8-bit floating point number might look like.

Sseemmmm
S    = sign of the mantissa (ie the sign of the whole number)
see  = exponent (the first bit is the sign bit for the exponent)
mmmm = mantissa

Here's a binary example using this scheme.

00101010
Sseemmmm
S      = 0: +ve
see    = 010: s = 0, +ve sign. exponent = +10 binary = 2. multiply mantissa by 2^2 = 4 decimal
mmmm   = 1010: = 1.010 binary = 1/1 + 0/2 + 1/4 + 0/8 = 1.25 decimal
number = + 1.25 * 2^2 = +5.00

Here's another example

11011101
Sseemmmm
S      = 1: -ve number
see    = 101: s = 1, -ve sign. exponent = -01 binary = -1 decimal. multiply mantissa by 2^-1 = 0.5 decimal
mmmm   = 1101: = 1.101 binary = 1/1 + 1/2 + 0/4 + 1/8 = 1.625 decimal
number = - 1.625 * 2^-1 = -0.8125

Doing the same problem with bc -l

11011101
Sseemmmm
#sign is -ve

#checking that I can do the mantissa
# echo "obase=10; ibase=2; 1.101 " | bc -l
1.625

#checking that I can do the exponent. "see" is 101 so exponent is -01
# echo "obase=10; ibase=2; (10)^(-01)" | bc -l
0.5

#the whole number
# echo "obase=10; ibase=2; -1.101 * (10)^(-01)" | bc -l
-0.8125

How many floating point numbers is it possible to represent using a byte (you don't know the value of these numbers or their distribution in the continuum, but you do know how many there are) ^[114] ? (How many reals can you represent using 32-bits ^[115] ?) We could choose any 8-bit scheme we want (e.g. seemmmmm, seeeemmm, mmeeeees, mmmmseee...), but we will still be subject to the 8-bit limit (we will likely have different reals in each set).

In the the seeemmmm scheme used here, where do the numbers come from? The first bit (the sign) has 2 possibilities, we have 3 bits of exponent giving 8 exponents, and 4 bits of mantissa giving us 16 mantissas. Thus we have 2*8*16=256 8-bit floating point numbers.

	Note
	End Lesson 17

Student Example: What real number is represented by the 8-bit floating point number 01101001 ^[116] ? Use bc to convert this number to decimal ^[117]

You aren't allowed to have 0 as the first number of the mantissa. If you get a starting 0 in a decimal number, you move the decimal point till you get a single digit 1-9 before the decimal point (this format is also called "scientific notation"). The process of moving the decimal point is called normalisation.

0.1234*10^3 non-normalised floating point number
1.234*10^2  normalised floating point number.

Is 12.34*10¹ normalised ^[118] ?

Let's see what happens if we don't normalise 8-bit reals. For the moment, we'll ignore -ve numbers (represented by the first bit of the floating point number). Here are the 8 possible exponents and 16 possible mantissas, giving us 8*16=128 +ve numbers.

   eee   E
   000   2^0  =1
   001   2^1  =2
   010   2^2  =4
   011   2^3  =8
   100   2^-0 =1
   101   2^-1 =1/2
   110   2^-2 =1/4
   111   2^-3 =1/8

   mmmm   M
   0000   0
   0001  1/8
   0010  1/4
   0011  3/8
   0100  1/2
   0101  5/8
   0110  3/4
   0111  7/8
   1000  1
   1001  1 + 1/8
   1010  1 + 1/4
   1011  1 + 3/8
   1100  1 + 1/2
   1101  1 + 5/8
   1110  1 + 3/4
   1111  1 + 7/8

Because we have mantissas <1.0, there are multiple representations of numbers, e.g. the number 1/4 can be represented by (1/8)*2 or (1/4)*1; 1 + 1/4 can be represented by (5/8)*2 or 1 + 1/4. We want a scheme that has only one representation for any number. Find other numbers that can be represent in multiple ways - twice ^[119] - trice ^[120] - four times ^[121] - five times ^[122] .

Normalisation says that the single digit before the decimal point can't be 0; we have to start the number with any of the other digits. For decimal the other digits are 1..9. In binary we only have one choice: "any of the other digits" is only 1. All binary floating point numbers then must have a mantissa starting with 1 and only binary numbers from 1.0 to 1.11111.... (decimal 1.0 to 1.99999...) are possible in the mantissa.

Since there's no choice for the value of the first position in the mantissa, real world implementations of floating point representations do not store the leading 1. The algorithm, which reads the floating point format, when it decodes the stored number, knows to add an extra "1" at the beginning of the mantissa. This allows the storage of another bit of information (doubling the number of reals that can be represented), this extra bit is part of the algorithm and not stored in the number.

Here are all the possible normalised mantissas in our 8-bit scheme (the leading 1 in the mantissa is not stored). How many different mantissas are possible in our scheme ^[123] ?

stored				binary		decimal
shortened	true		value		value
mantissa	mantissa	mantissa	mantissa
0000		10000		1.0000		1
0001		10001		1.0001		1 + 1/16 
.
.
1111		11111		1.1111		1 + 15/16

Here's a table of the +ve 8-bit normalised real numbers. (How many +ve 8-bit real numbers are there ^[124] ?) There are only 112 numbers in this list. The 16 entries with see=000 are missing: these are sacrificed to represent 0.0 - we'll explain this shortly - see the non-existance of 0.0. Since a byte can represent an int just as easily as it can represent a real, the table includes the value of the byte, as if it were representing an int.

In my home-made 8-bit real number, the exponent is a signed int. For clarity I've separated out the sign bit of the exponent. If you use the signed int convention, then the first bit of an integer is the sign, with 1 being -ve and 0 being +ve. Remember that if you're using the signed int representation, and you have a 3 bit int, then 000-1=fff.

S see mmmm	int 	M	E	real
0 001 0000	16	1.0	2.0	2.0
0 001 0001	17	1.0625	2.0	2.125
0 001 0010	18	1.125	2.0	2.25
0 001 0011	19	1.1875	2.0	2.375
0 001 0100	20	1.25	2.0	2.5
0 001 0101	21	1.3125	2.0	2.625
0 001 0110	22	1.3750	2.0	2.75
0 001 0111	23	1.4375	2.0	2.875
0 001 1000	24	1.5	2.0	3.0
0 001 1001	25	1.5625	2.0	3.125
0 001 1010	26	1.625	2.0	3.25
0 001 1011	27	1.6875	2.0	3.375
0 001 1100	28	1.75	2.0	3.5
0 001 1101	29	1.8125	2.0	3.625
0 001 1110	30	1.875	2.0	3.750
0 001 1111	31	1.9375	2.0	3.875
0 010 0000	32	1.0	4.0	4.0
0 010 0001	33	1.0625	4.0	4.25
0 010 0010	34	1.125	4.0	4.5
0 010 0011	35	1.1875	4.0	4.75
0 010 0100	36	1.25	4.0	5.0
0 010 0101	37	1.3125	4.0	5.25
0 010 0110	38	1.375	4.0	5.5
0 010 0111	39	1.4375	4.0	5.75
0 010 1000	40	1.5	4.0	6.0
0 010 1001	41	1.5625	4.0	6.25
0 010 1010	42	1.625	4.0	6.5
0 010 1011	43	1.6875	4.0	6.75
0 010 1100	44	1.75	4.0	7.0
0 010 1101	45	1.8125	4.0	7.25
0 010 1110	46	1.875	4.0	7.5
0 010 1111	47	1.9375	4.0	7.750
0 011 0000	48	1.0	8.0	8.0
0 011 0001	49	1.0625	8.0	8.5
0 011 0010	50	1.125	8.0	9.0
0 011 0011	51	1.1875	8.0	9.5
0 011 0100	52	1.25	8.0	10.0
0 011 0101	53	1.3125	8.0	10.5
0 011 0110	54	1.375	8.0	11.0
0 011 0111	55	1.4375	8.0	11.5
0 011 1000	56	1.5	8.0	12.0
0 011 1001	57	1.5625	8.0	12.5
0 011 1010	58	1.625	8.0	13.0
0 011 1011	59	1.6875	8.0	13.5
0 011 1100	60	1.75	8.0	14.0
0 011 1101	61	1.8125	8.0	14.5
0 011 1110	62	1.875	8.0	15.0
0 011 1111	63	1.9375	8.0	15.5
0 100 0000	64	1.0	1.0	1.0
0 100 0001	65	1.0625	1.0	1.0625
0 100 0010	66	1.125	1.0	1.125
0 100 0011	67	1.1875	1.0	1.1875
0 100 0100	68	1.25	1.0	1.25
0 100 0101	69	1.3125	1.0	1.3125
0 100 0110	70	1.375	1.0	1.375
0 100 0111	71	1.4375	1.0	1.4375
0 100 1000	72	1.5	1.0	1.5
0 100 1001	73	1.5625	1.0	1.5625
0 100 1010	74	1.625	1.0	1.625
0 100 1011	75	1.6875	1.0	1.6875
0 100 1100	76	1.75	1.0	1.750
0 100 1101	77	1.8125	1.0	1.8125
0 100 1110	78	1.875	1.0	1.875
0 100 1111	79	1.9375	1.0	1.9375
0 101 0000	80	1.0	0.5	0.5
0 101 0001	81	1.0625	0.5	0.53125
0 101 0010	82	1.125	0.5	0.5625
0 101 0011	83	1.1875	0.5	0.59375
0 101 0100	84	1.25	0.5	0.625
0 101 0101	85	1.3125	0.5	0.65625
0 101 0110	86	1.375	0.5	0.6875
0 101 0111	87	1.4375	0.5	0.71875
0 101 1000	88	1.5	0.5	0.750
0 101 1001	89	1.5625	0.5	0.78125
0 101 1010	90	1.625	0.5	0.8125
0 101 1011	91	1.6875	0.5	0.84375
0 101 1100	92	1.75	0.5	0.875
0 101 1101	93	1.8125	0.5	0.90625
0 101 1110	94	1.875	0.5	0.9375
0 101 1111	95	1.9375	0.5	0.96875
0 110 0000	96	1.0	0.25	0.25
0 110 0001	97	1.0625	0.25	0.265625
0 110 0010	98	1.125	0.25	0.28125
0 110 0011	99	1.1875	0.25	0.296875
0 110 0100	100	1.25	0.25	0.3125
0 110 0101	101	1.3125	0.25	0.328125
0 110 0110	102	1.375	0.25	0.34375
0 110 0111	103	1.4375	0.25	0.359375
0 110 1000	104	1.5	0.25	0.375
0 110 1001	105	1.5625	0.25	0.390625
0 110 1010	106	1.625	0.25	0.40625
0 110 1011	107	1.6875	0.25	0.421875
0 110 1100	108	1.75	0.25	0.4375
0 110 1101	109	1.8125	0.25	0.453125
0 110 1110	110	1.875	0.25	0.46875
0 110 1111	111	1.9375	0.25	0.484375
0 111 0000	112	1.0	0.125	0.125
0 111 0001	113	1.0625	0.125	0.132812
0 111 0010	114	1.125	0.125	0.140625
0 111 0011	115	1.1875	0.125	0.148438
0 111 0100	116	1.25	0.125	0.15625
0 111 0101	117	1.3125	0.125	0.164062
0 111 0110	118	1.375	0.125	0.171875
0 111 0111	119	1.4375	0.125	0.179688
0 111 1000	120	1.5	0.125	0.1875
0 111 1001	121	1.5625	0.125	0.195312
0 111 1010	122	1.625	0.125	0.203125
0 111 1011	123	1.6875	0.125	0.210938
0 111 1100	124	1.75	0.125	0.21875
0 111 1101	125	1.8125	0.125	0.226562
0 111 1110	126	1.875	0.125	0.234375
0 111 1111	127	1.9375	0.125	0.242188

a couple of things to notice

reals have this format:

fixed number of evenly spaced numbers from 1..2 multiplied by 2<superscript>integer_exponent</superscript>.

In the 8-bit scheme here, there are 16 evenly numbers in each 2:1 range. In the range 1..2 the interval is (2-1)/16=1/16=0.0625; the numbers then are 1, 1 + 1/16.. 1 + 15/16. For the range 2..4, you multiply the base set of numbers by 2. There are still 16 evenly spaced intervals and the resulting spacing is (4-2)/16 = 1/8 = 0.125.

In our scheme, how many 8-bit reals are there between 8..16 and what is the interval ^[126] ? How many numbers are there between 0.25..0.5 and what is interval ^[127] ?

Here's a graphical form of the chart of 128 real numbers above. To show the numbers as a 2-D plot

Only the +ve numbers are shown. The first 16 of the possible 128 +ve numbers are not plotted (these are sacrificed to represent the 0.0 and a few other numbers, see the non-existance of 0.0). The absence of these 16 numbers is seen as the gap on the left side, where there are no red dots between 0 and the smallest number on the x-axis.

Note

This is about not being able to read your own code 6 months after you write it: I wrote this section (on floating point numbers) about 6 months before I delivered it in class. (I initially expected I would use it in primitive data types, but deferred it in the interest in getting on with programming.) On rereading, the explanation for the gap on the left hand side, made no sense at all. I couldn't even figure out the origin of the gap from my notes. After a couple of weekends of thinking about it, I put in a new explanation, and then a week later I figured out my original explanation. It was true, but not relevant to the content.

Figure 1. graph of the 7 rightmost bits in a floating point number represented as an integer (X axis) and real (Y axis)

The plot is a set of dots (rather than a line or set of lines), showing that you can only represent some of the real space; a few points on the infinity of points that is a line.

	Note
	End Lesson 18

The representation of the +ve reals is

1..2 * 2^exponent.

The 8-bit real numbers (from 0.25 to 4.0) on a straight line look like this (because of crowding, reals<0.25 aren't represented)

Figure 2. graph of the 8-bit real numbers from 0.25..4 plotted on a linear real number line.

There are 16 evenly spaced 8-bit real numbers in any 2:1 section of the real number line (e.g. from 0.25..0.50, 0.50..1.0, 1.0..2.0, 2.0..4.0). Because the length of each segment decreases by half each step going leftward (towards 0.0), the real numbers never get to 0.0 (0.0 would be 2^-∞).

Logarithms:

If you don't know about logarithms

Logarithms are the exponents when reals are represented as base^exponent.

Logarithms do for multiplication, what linear plots do for arithmetic.

Let's say you want a graph or plot to show the whole range of 8-bit reals from 0.125=2^-3 to 16=2⁴. On a linear plot most of the points would be scrunched up at the left end of the line (as happens above). If instead you converted each number to 2^x and plotted x, then the line would extend from -3 to +4, with all points being clearly separated.

If n=2^x, then x is log₂(n). So log₂(1)=0, log₂(2)=1, log₂(4)=2.

If n=2^0.5, then n²=2^0.5*2^0.5=2. Then n=sqrt(2) and 2^0.5 is sqrt(2).

If you multiply two numbers, you add their logarithms.

Logarithms commonly use base 10 (science and engineering; the log of 100 is 2), base 2 (for computing), and log_e (for mathematics, e=2.718281828459)

The 8-bit real numbers (from 0.25 to 4.0) on a logarithmic line look like this

Figure 3. graph of the 8-bit real numbers from 0.25..4 plotted on a log real number line.

Note there are 16 logarithmically spaced real numbers in any 2:1 section of the real number line (e.g. from log₂(number) =-2..-1, -1..0, 0..1, 1..2. Because 0.0 is 2^-∞, 0.0 on this graph is off the left end of the line (at -∞).

How many real numbers can be represented by a 32 bit number ^[128] ? How many floating point numbers exist in the range covered by 32-bit reals ^[129] ?

You can make real numbers as big as you have bits to represent them and as small as you have bits to represent them. However you can't represent 0.0 on the simple floating point scheme shown here.

However 0.0 is a useful real number. Time starts at 0.0; the displacement of a bar from it's equilibrium position is 0.0.

In the above images/graphs, there is no way to represent 0.0 - it's not part of the scheme. Since a normalised decimal number starts with 1..9 and normalised binary numbers start with 1, how do you represent 0.0? You can't. You need some ad hoc'ery (Ad hoc. http://www.wikipedia.org/wiki/Ad_hoc).

To represent 0.0, you have to make an exception like "all numbers must start with 1-9₁₀ or 1₂, unless it's 0.0, in which case it starts with 0". Designers of rational systems are dismayed when they need ad hoc rules to make their scheme work. ad hoc rules indicate that maybe the scheme isn't quite as good as they'd hoped. Exceptions slow down the computer, as rather than starting calculating straight away, the computer first has to look up a set of exceptions (using a conditional statement) and figure out which execution path to follow.

How do we represent 0.0? Notice that the bits for the exponent in the previous section can be 100 or 000; both represent 1 (2⁰ or 2^-0). 8 of our exponents are duplicates. We could handle the duplication, by having the algorithm add 1 to the exponent if it starts with 1, giving us an extra 8 numbers (numbers that are smaller by a factor of 2). However since we don't yet have the number 0.0 anywhere, the convention (in this case of a 1 byte floating point number) is to have two of the 8 numbers with the exponent bits=000 as +/- 0.0. As well we need to represent +/-∞ and NaN, so we use some of the remaining 6 numbers for those.

S  see  mmmm 
0  000  0000 +0.0
1  000  0000 -0.0 (yes there is -0.0)
0  111  0000 +infinity
1  111  0000 -infinity
*  111 !0000 NaN

NaN (not a number): In the IEEE-754 scheme, infinities and NaN are needed to handle overflow and division by 0. Before the IEEE-754 standard, there was no representation in real number space for ∞ or the result of division by zero or over/underflow. Any operation which produced these results had to be trapped by the operating system and the program would be forced to exit. Often these results came about because the algorithm was exploring a domain of the problem for which there was no solution, with division by zero only indicating that the algorithm needed to explore elsewhere. A user would quite reasonably be greatly aggrieved if their program, which had been running happily for a week, suddenly exited with an irrecoverable error, when all that had happened was that the algorithm had wandered into an area of the problem where their were no solutions. The IEEE-754 scheme allows the results of all real number operations to be represented in real number space, so all math operations produce a valid number (even if it is NaN or ∞) and the algorithm can continue. Flags are set showing whether the result came about through over/underflow or division by 0, allowing code to be written to let the algorithm recover.

NaN is also used when you must record a number, but have no data or the data is invalid: a machine may be down for service or it may not be working. You can't record 0.0 as this would indicate that you received a valid data point of 0.0. Invalid or absent data is a regular occurrence in the real world (unfortunately) and must be handled. You process the NaN value along with all the other valid data; NaN behaves like a number with all arithmetic operations (+-*/; the result of any operation with NaN always being NaN). On a graph, you plot NaN as a blank.

The number -0.0: is needed for consistency. It's all explained in "What every computer scientist should know about floating-point arithmetic" (link above). It's simplest to pretend you never heard of -0.0. Just in case you want to find out how -0.0 works...

>>> 0.0 * -1
-0.0
>>> -0.0 * 1
-0.0
>>> -0.0 * -1
0.0
>>> -(+0.0)
-0.0		#required by IEEE754
>>> -(+0)		
0		#might have expected this to be -0
>>> -0.0 + 0.0
0.0
>>> -0.0 - 0.0
-0.0
>>> -0.0 + -0.0
-0.0
>>> -0.0 - +0.0
-0.0
>>> -0.0 - -0.0
0.0
>>> 0.0 + -0.0
0.0
>>> 0.0 - +0.0
0.0
>>> 0.0 - -0.0
0.0
>>> i=-0.0
>>> if (i<0.0):
...     print "true"
... 
>>> if (i>0.0):
...     print "true"
... 
>>> if (i==0.0):
...     print "true"
... 
true		#(0.0==-0.0) is true in python
>>> if (i==-0.0):
...     print "true"
... 
true
>>> if (repr(i)=="0.0"):
...     print "true"
... 
>>> if (repr(i)=="-0.0"):
...     print "true"
... 
true

The result of the test (0.0==-0.0) is not defined in most languages. (As you will find out shortly, you can't compare reals, due to the limited precision of real numbers. You can't do it even for reals that have the exact value 0.0).

How many numbers are sacrificed in the 8bit float point representation used here, and in a 32 bit floating point representation (23bit mantissa, 8 bit exponent, 1 bit sign) to allow representation of 0.0, ∞ and NaN ^[130] ? What fraction of the numbers is sacrificed in each case ^[131] ?

	Note
	(with suggestions from the Python Tutorial: representation error).

We've found that only a limited number of real numbers (256 in the case of 8-bit reals, 4G in the case of 32-bit reals) can be represented in a computer. There's another problem: most (almost all) of the reals that have finite representation in decimal (e.g.0.2, 0.6) don't have a finite representation in binary.

Only numbers that evenly divide (i.e. no remainder) powers of the base number (10decimal) can be represented by a finite number of digits in the decimal system. The numbers which evenly divide 10 are 2 and 5. Any number represented by (1/2)ⁿ*(1/5)^m can be represented by a finite number of digits. Any multiple of these numbers can also be represented by a finite number of digits.

In decimal notation, fractions like 1/10, 3/5, 7/20 have finite represenations

1/10=0.1
3/5=0.6
7/20=0.35

In decimal notation, fractions like 1/3, 1/7, 1/11 can't be represented by a finite number of digits.

1/3=0.33333....
1/7=0.14285714285714285714...

Do these (decimal) fractions have finite representations as decimal numbers: 1/20, 1/43, 1/50 ^[132] ?

Here's some numbers with finite represenation in decimal.

#(1/2)^3
1/8=0.125

#3/(2^3)
3/8=0.375

#note 1/80 is the same as 1/8 right shifted by one place. Do you know why?
#1/((2^4)*5)
1/80=0.0125

#1/(2*(5^2))
1/50=0.02

#7/(2*(5^2))
7/50=0.14

Using the above information, give 1/800 in decimal ^[133]

Divisors which are multiples of other primes (i.e. not 2,5) cannot be represented by a finite number of decimals. Any multiples of these numbers also cannot be represented by a finite number of digits. Here's some numbers which don't have finite represenation in decimal.

#1/(3^2)
1/9=0.11111111111111111111...

#1/(3*5)
1/15=0.06666666666666666666...

#4/(3*5)
4/15=0.26666666666666666666...

(the Sexagesimal (base 60) http://en.wikipedia.org/wiki/Sexagesimal numbering system, from the Sumerians of 2000BC, has 3 primes as factors: 2,3,5 and has a finite representation of more numbers than does the decimal or binary system. We use base 60 for hours and minutes and can have a half, third, quarter, fifth, sixth, tenth, twelfth, fifteenth, twentieth... of an hour.)

In binary, only (1/2)ⁿ can be represented by a finite number of digits. Fractions like 1/10 can't be prepresented by a finite number of digits in binary: there is no common divider between 10 (or 5) and 2.

# echo "scale=10;obase=2;ibase=10; 1/10" | bc -l
.0001100110011001100110011001100110

In the binary system, the only number that evenly divides the base (2), is 2. All other divisors (and that's a lot of them) will have to be represented by an infinite series.

Here's some finite binary numbers:

echo "obase=2;ibase=10;1/2" | bc -l
.1
echo "obase=2;ibase=10;3/2" | bc -l
1.1

#3/8 is the same as 3/2, except that it's right shifted by two places.
echo "obase=2;ibase=10;3/8" | bc -l
.011

From the above information give the binary representation of 3/32 ^[134]

Knowing the binary representation of decimal 5, give the binary representation of 5/8 ^[135]

Here's some numbers that don't have a finite represenation in binary

"obase=2;ibase=10;1/5" | bc -l
.0011001100110011001100110011001100110011001100110011001100110011001

#note 1/10 is the same as 1/5 but right shifted by one place
#(and filling in the empty spot immediately after the decimal point with a zero).
echo "obase=2;ibase=10;1/10" | bc -l
.0001100110011001100110011001100110011001100110011001100110011001100

with the above information, give the binary representation of 1.6 ^[136]

Here is 1/25 in binary

echo "obase=2;ibase=10;1/25" | bc -l
.0000101000111101011100001010001111010111000010100011110101110000101

From this information, what is the binary representatin of 1/100 ^[137]

Just as 1/3 can't be represented by a finite number of decimal digits, 1/10 can't be represented by a finite number of binary digits. Any math done with 0.1 will be done with a 32-bit truncated representation of 0.1

>>> 0.1
0.10000000000000001
>>> 0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1+0.1
0.99999999999999989
>>> 0.1*10
1.0

The problem of not being able to represent finite decimal numbers in a finite number of binary digits is independant of the language. Here's an example with awk:

echo "0.1" | awk '{printf "%30.30lf\n", $1}'
0.100000000000000005551115123126

echo "0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1" | awk '{printf "%.30lf\n", $1+$2+$3+$4+$5+$6+$7+$8+$9+$10}'
0.999999999999999888977697537484

Do not test if two reals are equal. Even though the two numbers (to however many places you've displayed) appear the same, the computer may tell you that they're different. Instead you should test

abs(your_number - the_computers_number) < some_small_number

For the real way to do it see Lahey compiler.

Here's an example of how things don't work quite the way you'd expect.

>>> i=0.1
>>> if (i==0.1):	#is (i==0.1)?
...     print "true"
... 
true
>>> j=0
>>> for x in range(0,10):
...     j += i
... 
>>> j			#j should be 1.0, but you loose a 1 in the last place every time you add 0.1
0.99999999999999989
>>> print j		#outputs field width of 12 (and will round the number)
1.0
>>> print "%20.12f" %j	#print with explicit 12 digits after decimal point
      1.000000000000
>>> print "%20.20f" %j	#wider field
0.99999999999999988898
>>> print 1.0-j
1.11022302463e-16
>>> if (j==1.0):	#is (j==1.0)?
...     print "true"
... else:
...     print j
... 
1.0

Is (j==1.0) true? The answer is not yes or no; you aren't going to get a sensible answer. Knowing that, you should never ask such a question. In case you aren't convinced, let's see what might happen if you do.

(About 747s Airliner Takeoff Speeds, Boeing 747 http://www.aerospaceweb.org/aircraft/jetliner/b747/ Boeing 747-400 http://en.wikipedia.org/wiki/Boeing_747-400 http://www.aerospaceweb.org/question/performance/q0088.shtml). About diesel locomotives Electric Diesel Hauling II (http://www.cfr.ro/jf/engleza/0310/tehnic.htm) Steam Ranger Enthusiast Pages (http://www.steamranger.org.au/enthusiast/diesel.htm).

You've written the computer code to automatically load the Boeing 747-400ER cargo planes each night at the FedEx Memphis terminal. You have a terminal with 1500 palettes, each palette weighing 0.1Mg (about 200lbs). The max cargo payload is 139.7Mg.

It's 2am. All the packages have arrived by truck and are now sitting on palettes. Let's start loading cargo.

#! /usr/bin/python

palettes_in_terminal=[]	#make an empty list for the palettes in the terminal awaiting loading
palettes_in_747=[]	#make another empty list for the palettes loaded onto the plane
number_palettes=1500
weight_palette=0.1
weight_on_plane=0.0
max_weight_on_plane=139.7

#account for palettes in terminal
#make a list containing 1500 entries of 0.1
for i in range(0,number_palettes):
	palettes_in_terminal.append(weight_palette)
	

#load plane
#while plane isn't full and there are more palettes to load
#	take a palette out of the terminal and put it on the plane
while ((weight_on_plane != max_weight_on_plane) and (len(palettes_in_terminal) != 0)):
	#the next line is a train wreck version of this
	#moved_palette=palettes_in_terminal.pop()
	#palettes_in_747.append(moved_palette)
	palettes_in_747.append(palettes_in_terminal.pop())
	weight_on_plane += weight_palette
	#print "weight of cargo aboard %f" % weight_on_plane

print "weight on plane %f" % weight_on_plane
print "max weight on plane %f" % max_weight_on_plane
#-------------
	

What's this code doing? Swipe the code and run it. Is the plane safe to fly? How overloaded is it? How did the code allow the plane to become overloaded ^[138] ?

In an hour (at 3am), the plane, with a crew of 2, weighing 412,775 kg (about 400 tons, or 2-8 diesel locomotives), will be sitting at the end of the runway ready to take off, with its cargo and enough jet fuel (241,140l.) to make it across the country in time for the opening of the destination airport at 6am. (Why do the pilots have to wait till 6am, why couldn't they fly in earlier ^[139] ?) The pilot will spool up the engines, hold the plane for several seconds to check that the engines are running at full power, then let off the brakes, when the plane will start to roll. After thundering 3km down the runway, about a minute later alarms will alert the pilot that the apparently fully functional plane hasn't reached the lift off speed of 290 km/hr. The plane is too far down the runway to abort takeoff ^[140]. As far as the pilot knows, there is nothing wrong with the plane, only the speed is low. Is the indicated low speed an error and the the plane is really moving at the expected speed? Should the pilot rotate the plane and attempt lift off? He has only seconds to decide. Assuming the runway is long enough for the plane to eventually reach lift off speed, will it lift off ^[141] ? You'll read all about it in the papers next morning. In fact the pilot has no options and no decision to make. He was committed when he let the brakes off.

Why aren't there crash barriers at the end of the runway ^[142] ? There are better solutions for handling an overweight plane (see below).

The airline will blame the software company. The owners of the software business will protest that it was a problem with one of the computers that no-one could possibly have anticipated (there are people that will believe this). The american legal system will step in, the business will declare bankruptcy, the owners will have taken all their money out of the business long ago and have enough money to live very comfortably, the workers will be out of a job, and will find their pension fund unfinanced. The owners of the software business will be sued by the airline and the workers, but since the business is bankrupt, there is no money even if the suit is successful. Only the lawyers and the business owners will come out of this in good financial shape.

No-one would ever write such a piece of software. Even if you could equate reals, you wouldn't have written the code this way. Why not ^[143] ? Fix the code (you need to change the conditional statement). Here's my fixed version ^[144] .

While the algorithm I used here is completely stupid, the actual method that the FAA uses to calculate the take-off weight of a passenger plane is more stupid (and I'm not making this up): the airline just makes up a value. They guess how much the passengers weigh (and from watching how my luggage is handled at check in, they guess the weight of the baggage). The final decision is left to the pilots (at least in the case of US AirWays 5481 below), who have to decide by looking at the baggage and the passengers, whether it's safe to take off. This is the method that's used in the US. If you were a pilot, could you estimate the weight of passengers to 1% by looking out the cockpit door, down the aisles at the passengers? Pilots are trained to fly planes not to weigh passengers: scales are used to weigh planes.

Considering the consequences of attempting to take off overweight, it would be easy enough to weigh a plane before take-off, as is done for a truck or train: see Weigh Bridge (http://en.wikipedia.org/wiki/Weigh_bridge). Accuracy of current equipment is enough for the job: see Dynamass (www.pandrol.co.za/pdf/weigh_bridges.pdf) which can weigh a train moving at 35kph to 0.3% accuracy.

Planes attempt to take off overweight all the time.

The plane piloted by 7 yr old Jessica Dubroff, attempting to be the youngest girl to fly across the USA, crashed shortly after take off from Cheyenne (see Jessica Dubroff. http://en.wikipedia.org/wiki/Jessica_Dubroff). The plane was overloaded and attempted to take off in bad weather.

Most plane crashes due to overloading occur because the weight is not distributed correctly.

Air Midwest #5841 (http://en.wikipedia.org/wiki/US_Airways_Express_Flight_5481) out of Charlotte NC in Jan 2003 crashed on take off, because (among other things) the plane was overloaded by 600lbs, with the center of gravity 5% behind the allowable limit. The take-off weight was an estimate only, using the incorrect (but FAA approved) average passenger weight estimate, that hasn't been revised since 1936. The supersized passengers of the new millenium are porkers, weighing 20lbs more. The airline washed their hands of responsibility and left it to the pilots to decide whether to take off. Once in the air, the plane was tail heavy and the pilots couldn't get the nose down. The plane kept climbing till it stalled and crashed.

An Emirates plane taking off from Melbourne (http://www.smh.com.au/travel/takeoff-error-disaster-averted-20090430-aozu.html) with 257 passengers, weighed 100tonnes more than the pilots had been told. The Airbus A340-500 scraped its tail along the tarmac and grassland beyond the runway at Tullamarine on March 20, then hit airport landing lights and disabled a radio antenna before taking off. The pilots then dumped fuel over Port Phillip Bay for about 30 minutes to reduce the weight before making an emergency landing at Tullamarine.

One would think that 32-bits is enough precision for a real, and in many cases (just adding or multiplying a couple of numbers) it is. However most of the time spent on computers is optimising multi-parameter models. Much scientific and engineering computing is devoted to modelling (making a mathematical representation for testing - weather, manufactured parts, bridges, electrical networks). You are often optimising a result with respect to a whole lot of parameters; e.g.

Optimisation is similar an ant's method for finding the top of a hill. All the ant knows is its position and it's altitude. It moves in some direction and takes note of whether it has gone uphill or downhill. When it can no longer move and gain altitude, the ant concludes that it is at the top of the hill and has finished its work. The area around the optimum is relatively flat and the ant may travel a long way in the x,y direction (or whatever the parameters are), before arriving at the optimum. Successful optimisation relies on the ant being able to discern a slope in the hill.

Because of the coarseness of floating point numbers, the hill instead of being smooth, is a series of slabs (with consecutive slabs of different thickness). An ant standing next to the edge of a slab, sees no change in altitude if it moves away from the edge, but if it had moved towards the edge would have seen a large change in altitude. If the precision of your computer's floating point calculations is too coarse, the computer will not find the optimum. In computerese - your computation will not converge to the correct answer and you won't be able to tell that you've got the wrong solution, or it will not converge at all - the ant will wander around a flat slab on the side of the hill, never realising that it's not at the top.

If the computation changes from 32 to 64-bit, how many slabs will now be in the interval originally covered by a 32-bit slab ^[145] ?

Because you can't rely on a 32-bit computation converging or getting the correct answer, scientists and engineers have been using 64-bit (or greater) computers for modelling and optimisation for the past 30yrs. You should not use 32 bit floating point for any work that involves repetitive calculation of small differences.

How do you know that your answer is correct and not an artefact of the low precision computations you're using? The only sure way of doing this is to go to the next highest precision and see if you get the same answer. This solution may require recoding the problem, or using hardware with more precision (which is likely unavailable, you'll already be using the most expensive computer you can afford). Programmers test their algorithm with test data for which they have known answers. If you make small changes in the data, you should get small changes in the answers from the computer. If the answers leap around, you know that the algorithm did not converge.

Money (e.g. $10.50) in a bank account looks like a real, but there is no such thing as 0.1c. Money needs to be represented accurately in a computer in multiples of 1c. For this computers use binary coded decimal (http://en.wikipedia.org/wiki/Binary-coded_decimal). This is an integer representation, based on 1c. You can learn about on this your own if you need it.

	Note
	End Lesson 19

The smallest number representable by a 64-bit real is 2.2250738585072020*10^-308. Its value is mainly limited by the number of digits you can store in the exponent (multiplied by a normalised mantissa whose value has the range 1.0..1.999).

Another small number of interest to computer programmers is the value of the smallest change you can make to a number (let's say number==1.0). This is determined by the number of bits in the mantissa. If the 23 bit mantissa for 1.0 on a 32 bit machine (being normalised) is 000000000000000000000000, then the smallest change will be to 000000000000000000000001 (only the last bit changes). The number represented by this difference is called the machine's epsilon (http://en.wikipedia.org/wiki/Machine_epsilon). The number is called the "machine's epsilon" because before IEEE 754, every machine had a different way of representing reals, and each machine would have its own epsilon. Now with most machines (and except for specialised embedded machines) being IEEE 754 compliant, the machine's epsilon will be the IEEE 754 epsilon. Let's ask our machine what its epsilon is.

>>> epsilon = 1		#has a simple representation in binary
>>> while (epsilon+1>1):
...     epsilon = 0.5*epsilon	#still has a simple representation in binary 
...     print epsilon
... 
0.5
0.25
0.125
.
.
7.1054273576e-15
3.5527136788e-15
1.7763568394e-15
8.881784197e-16
4.4408920985e-16
2.22044604925e-16
1.11022302463e-16

If to the computer, 1+epsilon==1, the loop exits. The machine's epsilon for 64-bit IEEE 754 reals is 1.11022302463e-16 (for comparison, the smallest representable real is 2.2250738585072020*10^-308). (note: epsilon is a lot bigger than the smallest number representable by the machine's reals. Do you understand the difference between these two numbers?)

What this section is saying is that the computer can't differentiate reals that are more closely spaced than the spacing represented by the last bit in the mantissa.

After shifting x=1/5 one place to the right (in the previous section) to give the result x=1/10, note that the right hand "1" falls off the end and is lost. In the following sequence of calculations (right shift, losing the last 1, then left shift, which will insert a 0, in the place which previously had a 1)

x=0.2
y=x/2	#y=0.1
z=y*2	#z=0.2

would the resulting z have the same value as the original x? If you were using fixed point operations, you would get a different value. However a computer uses a mantissa and an exponent to represent a real, so dividing by a power of 2 just changes the exponent and no precision is lost. The mantissa will be unchanged following the above operations, and z would have the same value as x.

>>> x=0.2
>>> y=x/2
>>> z=y*2
>>> print z-x
0.0

Taking this example to the extreme...

>>> n_128=2**128
>>> print n_128
340282366920938463463374607431768211456
>>> x=0.1
>>> y=x/n_128                                  
>>> z=y*n_128
>>> print "%30.30g, %30.30g, %30.30g" %(x,y,z)
0.100000000000000005551115123126, 2.93873587705571893305445304302e-40, 0.100000000000000005551115123126

What if you did an operation that changed the mantissa (ignoring whether or not the exponent was changed)? We need a number that doesn't have a finite representation in binary.

>>> x=0.1		#x doesn't have a finite represenation in binary
>>> n_prime=149.0	#1/n_prime doesn't have finite representation in binary either
>>> y=x/n_prime
>>> z=y*n_prime
>>> print "%30.30g, %30.30g, %30.30g, %30.30g " %(x,y,z,z-x)
0.100000000000000005551115123126, 0.000671140939597315508424735241988, 0.100000000000000005551115123126, 0 

I would have thought that z and x would be different, but you get back the same number. I puzzled about this, till a co-worker Ed Anderson suggested the following as a method to find out what was going on: try a whole bunch of numbers. In particular, try some irrational numbers, like square roots (irrational numbers cannot be represented by a finite number of digits, in any base system). All square roots, except integers, are irrational (for the range 1..10, only 1,4 and 9 would have rational square roots). So use all the square roots from 1..some_large_number for both x and the divisor/multiplier. Ed's initial test was a 100x100 grid. Here's Ed's code for a 64*64 grid. There is a blank when the numbers are identical, and a letter/number if they were different.

#! /usr/bin/python

# fp_ed.py
# Joseph Mack (C) 2008, GPL v3.
# With suggestions from Ed Anderson

# takes a number, divides it by a 2nd, stores the result, then multiplies the result the 2nd number.
# in a perfect world, the result should be the original answer.
# due to limited precision of computer reals, sometimes a different answer will be returned.

from math import sqrt;
number_that_are_different=0;

end = 100;
for j in range(1,end):
	output_string = "" # can't output strings without \n in python and can't supress a blank either
	output_string += str('%3d ' %j);        #put blank after string
	for i in range(1,end):
	        x = sqrt(i);
	        y = x/j;
	        z = y*j;
	        if (x-z == 0.0):
	                output_string += " "

	        else:
	                #print x-z ;
	                output_string += "x";
	                number_that_are_different += 1;


	print output_string

print number_that_are_different
#- fp_ed.py -----------------------

Here's the output

	    1         2         3         4         5         6    
  1                                                                
  2                                                                
  3             x2                      33       3     w   3  3    
  4                                                                
  5  1     2     2               333     w                 3  3    
  6             x2                      33       3     w   3  3    
  7              x                                    w w3 w  3    
  8                                                                
  9                     w         w                                
 10  1     2     2               333     w                 3  3    
 11          x  2                     w 3w3w ww  w     3w 3  3     
 12             x2                      33       3     w   3  3    
 13             xx                               ww  w w   w w   3w
 14              x                                    w w3 w  3    
 15               x                                          w3w   
 16                                                                
 17                               w                                
 18                     w         w                                
 19          x                        w w ww                       
 20  1     2     2               333     w                 3  3    
 21                              w    3 w   3    3        3   3    
 22          x  2                     w 3w3w ww  w     3w 3  3     
 23               x                      3  3   3 w          www3  
 24             x2                      33       3     w   3  3    
 25          x   x                         w        w   w  w       
 26             xx                               ww  w w   w w   3w
 27              xx                              w   ww w  w3  w   
 28              x                                    w w3 w  3    
 29              2x                                        33w w3w 
 30               x                                          w3w   
 31                                                               w
 32                                                                
 33                                                                
 34                               w                                
 35       2                    3  3  3w                       3    
 36                     w         w                                
 37   y  2     x          33      w     w          w               
 38          x                        w w ww                       
 39       2                   33         3       33             w  
 40  1     2     2               333     w                 3  3    
 41              x              3 3 3               w      w    3  
 42                              w    3 w   3    3        3   3    
 43  1y    2 x2x  2              w 3  w    w   3   ww   w 3 w 33   
 44          x  2                     w 3w3w ww  w     3w 3  3     
 45                                 3                w  w    3   3 
 46               x                      3  3   3 w          www3  
 47         2x2 2                     w333 w   3     w33 w    w    
 48             x2                      33       3     w   3  3    
 49 0  1           2                        w w       3   3  3w    
 50          x   x                         w        w   w  w       
 51           2  x                          3ww33   w   3  w       
 52             xx                               ww  w w   w w   3w
 53              22                              3   3  3  3   3 3 
 54              xx                              w   ww w  w3  w   
 55   y        xx                                  w3 3w    3 3 w33
 56              x                                    w w3 w  3    
 57              2                                      3ww3     w 
 58              2x                                        33w w3w 
 59              2x                                       33w3 w   
 60               x                                          w3w   
 61                                                             w w
 62                                                               w
 63                                                                

What if we do the multiplication first

	    1         2         3         4         5         6    
  1                                                                
  2                                                                
  3   y        x  x                  w    w       3w  3  3     w w 
  4                                                                
  5           2  x                           3 3w       w  ww     3
  6   y        x  x                  w    w       3w  3  3     w w 
  7           x                              w w3           3      
  8                                                                
  9             x x                                    w 3w    w3 3
 10           2  x                           3 3w       w  ww     3
 11                                  3                w   3   w w3 
 12   y        x  x                  w    w       3w  3  3     w w 
 13  y     x 2  2               3 3w  w    3     3     3  w     3  
 14           x                              w w3           3      
 15     x              w                                  3     w  
 16                                                                
 17                                                          w3  w3
 18             x x                                    w 3w    w3 3
 19              22                               w   3 3w 3w3 3 w3
 20           2  x                           3 3w       w  ww     3
 21   y      2 x2                        3w3 w  3  w   3    3w    3
 22                                  3                w   3   w w3 
 23           2  2                  33         3        3  3 33    
 24   y        x  x                  w    w       3w  3  3     w w 
 25                                 ww  ww                   w     
 26  y     x 2  2               3 3w  w    3     3     3  w     3  
 27   1   2    2          w  w 3  w 3     3      w 3         3     
 28           x                              w w3           3      
 29                     w           w                        w3    
 30     x              w                                  3     w  
 31                                 3                        3     
 32                                                                
 33                                                             w3 
 34                                                          w3  w3
 35                                                      3  3      
 36             x x                                    w 3w    w3 3
 37             2                                      3w3    w  3 
 38              22                               w   3 3w 3w3 3 w3
 39   y       xx 2                             ww3 w  w  ww3  3 3  
 40           2  x                           3 3w       w  ww     3
 41             x x                       3   w 333    w  3    www 
 42   y      2 x2                        3w3 w  3  w   3    3w    3
 43               2                           3  ww   3 w      3   
 44                                  3                w   3   w w3 
 45   y      2xx                      w  w 3  ww w w               
 46           2  2                  33         3        3  3 33    
 47  y     x    2x                 w         33        3 3www   3  
 48   y        x  x                  w    w       3w  3  3     w w 
 49          x                   3    3    w ww         3w  3      
 50                                 ww  ww                   w     
 51   1   2   22               3        3     33   3          w    
 52  y     x 2  2               3 3w  w    3     3     3  w     3  
 53              2x              33           w   3        3   w   
 54   1   2    2          w  w 3  w 3     3      w 3         3     
 55      x               w w             3    3 w                3 
 56           x                              w w3           3      
 57                     3                     w       3 w          
 58                     w           w                        w3    
 59          x                        3    w  3          w         
 60     x              w                                  3     w  
 61                                           w                    
 62                                 3                        3     
 63                                           3                    

The figures are different. It seems that multiplication and division isn't always commutative.

Surprisingly (to me) there is a pattern, but then I didn't understand what was going on so any result would have been a surprise. "x"s appear in blocks on boundaries which are multiples of 8,16,32 and 64 (increasing the value of end shows that the block pattern continues).

With a slight modification of the code, the differences, if not 0.0, were all seen to be (1.11022302463e-16)*2ⁿ.

Is there anything in particular about 1.11022302463e-16? It turns out that it's the machine's epsilon (http://en.wikipedia.org/wiki/Machine_epsilon) (also see machine's epsilon).

What is 1.11022302463e-16 in binary?

dennis:~# echo "obase=2;1.11022302463*10^-16" | bc -l
.0000000000000000000000000000000000000000000000000000011111111111111

If you round this up by one bit, it's 2^-53. A 64 bit double has a 53 bit mantissa, indicating that the epsilon calculation is using a 64 bit real.

From the wiki entry, a single precision (32bit) real has a 24 bit mantissa. 1.0 then has an exponent of 0 and mantissa of 1.00000000000000000000000₂. The next biggest number has an exponent of 0 and a mantissa of 1.00000000000000000000001₂. The difference is 0.00000000000000000000001₂ or 2^-23.

Using IEEE 754 Converter (http://www.h-schmidt.net/FloatApplet/IEEE754.html), to do the conversion to a 32-bit real, I found

1.11022302463e-16=0x2500000=00100101000000000000000000000000. 

Breaking the binary representation into sign, exponent and mantissa (with the implied 1 for a normalised mantissa) gives

1.11022302463e-16 = 0 01001010 (1)00000000000000000000000
    decimal         s   exp             mantissa

Calculating this all gives: The sign bit: (bit 32, the leftmost bit) is 0 (positive).

The exponent: the IEEE-754 real math convention doesn't use the signed int representation. a 32 bit real, the exponent is a 9 bit integer (here 001001010). Floating point math is usually done on floating point hardware (or emulated on the the integer registers of the CPU, but the incremental cost of adding a floating point processor/unit (FPU) is small, so all x86 chips, starting with the 486, have a floating point processor built-in to the CPU chip). Since real math doesn't use the CPU's registers (each composed of 2,4,8 or 16 bytes), but instead uses the math co-processor's registers, then real math doesn't have to use the CPU's conventions (e.g. a signed int) to process a 9 bit int.

For 32 bit floating point, the convention is to subtract 127 from the value of the exponent. is is called "excess-127" (excess-1023 for 64 bit) and allows the exponent and the mantissa, when taken together, to vary monotonically (see IEEE-754 References http://babbage.cs.qc.edu/courses/cs341/IEEE-754references.html). (What this gets you, I don't know.) The 8-bit scheme 8-bit set of normalised reals I made up for this class is not monotonic in this respect (see the table and the figure, where there are jumps in value every 16 numbers).

In the current example, the exponent is -53 (from the next 9 bits). To get the actual value of the exponent, you subtract 127.

echo "obase=10; ibase=2; 001001010" | bc
74
then subtract 127
74-127=-53

The FPU int then is an unsigned int, from which is subtracted 127, so that -ve numbers are possible (it's called signed magnitude). Here are a few comparison numbers (assuming an 8 bit int):

bits		value, signed int	value, unsigned int	value, signed magnitude
00000000	0 			0 			-127
00000001	1			1			-126
01111111	127			127			0
10000000	-127			128			1
11111111	-1			255			128				

Since all of the evaluations can done with the same speed, with hardware designed to do these calculations, there is nothing to pick between them. Why IEEE-754 uses their particular representation of the exponent, rather than the signed int representation, I don't know.

The mantissa is 1.0 (after adding the implied 1 at the beginning of 23 zeroes).

Multiplying this all out gives the original number.

echo "1*2^-53" | bc -l
.00000000000000011102
#just checking that I can count 0s
 echo "1*2^-53/10^-16" | bc -l
1.11020000000000000000

The following table shows the number of times an entry occured for a grid size. The 2nd entry says that in an 16*16 grid, that 206 numbers were not changed and that -(1.11022302463e-16)*2³ was returned 9 times, +(1.11022302463e-16)*2*2 was returned 2 times and that +(1.11022302463e-16)*2*3 was returned 8 times. In the 16*16 grid 15*15=225 numbers are tested (=9+206+2+8).

grid size\n	-8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8

8*8		0 0 0 0 0 0 0 0 0 0 0 0 48 0 1 0 0 0 0 0 0 0 0 0 0
16*16		0 0 0 0 0 0 0 0 0 9 0 0 206 0 2 8 0 0 0 0 0 0 0 0 0
32*32		0 0 0 0 0 0 0 0 6 22 0 0 911 0 3 13 6 0 0 0 0 0 0 0 0
64*64		0 0 0 0 0 0 0 0 163 44 3 0 3566 1 6 33 153 0 0 0 0 0 0 0 0
128*128		0 0 0 0 0 0 0 99 346 82 10 0 15062 4 14 80 349 83 0 0 0 0 0 0 0
256*256		0 0 0 0 0 0 0 3085 723 163 24 0 56990 14 39 174 750 3063 0 0 0 0 0 0 0
512*512		0 0 0 0 0 0 1679 6215 1505 355 56 0 241228 35 87 381 1570 6360 1650 0 0 0 0 0 0
1024*1024	0 0 0 0 0 0 50447 12720 3101 752 121 0 911185 91 196 799 3209 12807 51101 0 0 0 0 0 0
2048*2048	0 0 0 0 0 26941 101676 25687 6251 1500 250 0 3863963 230 452 1669 6449 25610 102498 27033 0 0 0 0 0
4096*4096	0 0 0 0 0 822872 204267 51363 12576 3001 525 0 14575322 504 969 3366 12919 51254 205106 824981 0 0 0 0 0	
8192*8192	0 0 0 0 432519 1647991 409597 102641 25183 5958 1057 0 61830620 1105 2100 6885 25992 103109 411309 1652318 434097 0 0 0 0
16384*16384	0 0 0 0 13222697 3300406 821156 205169 50246 11827 2083 0 233145705 2345 4366 13903 52320 206826 824241 3306218 13233181 0 0 0 0 
32768*32768     0 0 0 6931293 26447787 6604323 1645180 409899 100085 23405 4099 0 989292832 4872 9022 28195 105323 414810 1651250 6614158 26462945 6926811 0 0 0
65536*65536	0 0 0 211603215 52895451 13214901 3295139 820314 200000 46646 8163 0 3730616481 9841 18209 56688 211726 830962 3306011 13229067 52912230 211561181 0 0 0

131072*131072
262144*262144
524288*524288
1048576*1048576	

8*8		            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.980 0.0 0.020 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
16*16		        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.040 0.0   0.0 0.916 0.0 0.009 0.036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
32*32		      0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.006 0.023 0.0   0.0 0.948 0.0 0.003 0.014 0.006 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
64*64		      0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.041 0.011 0.001 0.0 0.898 0.0 0.002 0.008 0.039 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
128*128		    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.006 0.021 0.005 0.001 0.0 0.934 0.0 0.001 0.005 0.022 0.005 0.0 0.0 0.0 0.0 0.0 0.0 0.0
512*512		  0.0 0.0 0.0 0.0 0.0 0.0 0.006 0.024 0.006 0.001 0.0   0.0 0.924 0.0 0.0   0.001 0.006 0.024 0.006 0.0 0.0 0.0 0.0 0.0 0.0
1024*1024	  0.0 0.0 0.0 0.0 0.0 0.0 0.048 0.012 0.003 0.001 0.0   0.0 0.871 0.0 0.0   0.001 0.003 0.012 0.049 0.0 0.0 0.0 0.0 0.0 0.0
2048*2048	0.0 0.0 0.0 0.0 0.0 0.006 0.024 0.006 0.001 0.0   0.0   0.0 0.922 0.0 0.0   0.0   0.002 0.006 0.024 0.006 0.0 0.0 0.0 0.0 0.0
4096*4096	0.0 0.0 0.0 0.0 0.0 0.049 0.012 0.003 0.001 0.0   0.0   0.0 0.869 0.0 0.0   0.0   0.001 0.003 0.012 0.049 0.0 0.0 0.0 0.0 0.0
8192*8192	0.0 0.0 0.0 0.0 0.006 0.025 0.006 0.002 0.0 0.0 0.0 0.0 0.922 0.0 0.0 0.0 0.0 0.002 0.006 0.025 0.006 0.0 0.0 0.0 0.0 
16384*16384	0.0 0.0 0.0 0.0 0.049 0.012 0.003 0.001 0.0 0.0 0.0 0.0 0.869 0.0 0.0 0.0 0.0 0.001 0.003 0.012 0.049 0.0 0.0 0.0 0.0
32768*32768     0.0 0.0 0.0 0.006 0.025 0.006 0.002 0.0 0.0 0.0 0.0 0.0 0.921 0.0 0.0 0.0 0.0 0.0 0.002 0.006 0.025 0.006 0.0 0.0 0.0
65536*65536	0.0 0.0 0.0 0.049 0.012 0.003 0.001 0.0 0.0 0.0 0.0 0.0 0.869 0.0 0.0 0.0 0.0 0.0 0.001 0.003 0.012 0.049 0.0 0.0 0.0

131072*131072
262144*262144
524288*524288
1048576*1048576	

After the multiplication/division, how many numbers are returned with a different value? For any particular j you have a 50% chance of losing a 0 off the right end, in which case you'll get back the original answer. So one possibility is 50% will be different. How about the divisor?; you only have a 50% chance of getting that right when you drop off extra digits. With two ways of getting it wrong, you only have 25% chance of getting it right, so 75% will be different. Here's the results.

grid_size	number_mul/div	number_changed	%changed	largest_difference
10*10		10^2		      2	 	  2		+    4 * 1.11e-16
100*100		10^4		    713		  7		+-  16 * 1.11e-16
1000*1000	10^6		 127980		 13		+-  32 * 1.11e-16	
10000*10000     10^8		8060303           8		+- 128 * 1.11e-16	
100000*100000	10^10	      927798399		  9		+- 512 * 1.11e-16

Let's look at this example with the IEEE-754 converter referred to above

>>> sqrt(14)
3.7416573867739413
>>> x=sqrt(14)
>>> y=x/3
>>> z=y*3
>>> z-x
-4.4408920985006262e-16

#from IEEE 754 converter
x = 3.7416573867739413 = 2 * 1.8708287 = 0 10000000 (1) 11011110111011101010001 
3.0 = 2 * 1.5                          = 0 10000000 (1) 10000000000000000000000
#from python
y = x/3.0 = 1.247219128924647          = 0 01111111 (1) 00111111010010011100000

#from bc
echo "obase=10;ibase=2; 1.11011110111011101010001/1.10000000000000000000000" | bc -l
1.24721916516621907552
echo  "obase=2;ibase=2; 1.11011110111011101010001/1.10000000000000000000000" | bc -l
1.0011111101001001110000010101010101010101010101010101010101010101010

#from python
z = 3.7416573867739409

from bc

echo "obase=10;ibase=2; 1.0011111101001001110000010101010101010101010101010101010101010101010 *11" | bc -l
3.7416574954986572265489474728439311945749068399891257286071777343750 

 
echo "obase=2;ibase=2; 1.0011111101001001110000010101010101010101010101010101010101010101010 *11" | bc -l
11.101111011101110101000011111111111111111111111111111111111111111111

The result after right shifting by one (to give an exponent of 1), and dropping the implied leading 1 from the mantissa give
a mantissa of 

11011110111011101010000	#new mantissa (after dropping extra digits)
11011110111011101010001 #original mantissa
	              1 #difference

since the exponent is 1 (a factor of 2) then the difference is

First we need to know the precision (number of significant figures) that python uses to represent real numbers. (What is the precision for integers ^[146] ?) Remember that a 32 bit computer can manipulate 2³² integers (you do remember the value of 2³²? - it's about 10 decimal digits). We'll learn more about Real Numbers later, but for the moment we won't go far wrong if we assume that a 32 bit computer is capable of representing 2³² real numbers. While there is only one integer between 0..2, there can be an infinite number of reals 0..2 and a computer can't hope to represent all of them. The computer finds the nearest one it can represent, whether it's correct or not, and say "it's near enough". If we go to 64-bits, we'll get a more accurate representation, but it still won't be exact.

Here's some code to show that floating point numbers longer than 12 digits are truncated i.e. any following numbers are garbage.

>>> for x in range (1,15):
...     print  x, 1+10**-x
... 
1 1.1
2 1.01
3 1.001
4 1.0001
5 1.00001
6 1.000001
7 1.0000001
8 1.00000001
9 1.000000001
10 1.0000000001
11 1.00000000001
12 1.0
13 1.0
14 1.0

10¹²12 is 2^{? [147]} This is not enough accuracy for a 64 bit number and is too much for a 32 bit number. We should suspect that we've goofed. The above truncation turns out to be due to limitations of the formatting. print uses str() for outputting, which only outputs 12 digits after the decimal point. If we increase the number of digits displayed, you can get about 15 significant figures before getting garbage. (Note: we are finding the machine's epsilon).

>>> for x in range (1,20):
...     print "%2d, %10.40f" %( x, 1+10**-x)
... 
 1, 1.1000000000000000888178419700125232338905
 2, 1.0100000000000000088817841970012523233891
 3, 1.0009999999999998898658759571844711899757
 4, 1.0000999999999999889865875957184471189976
 5, 1.0000100000000000655120402370812371373177
 6, 1.0000009999999999177333620536956004798412
 7, 1.0000001000000000583867176828789524734020
 8, 1.0000000099999999392252902907785028219223
 9, 1.0000000010000000827403709990903735160828
10, 1.0000000001000000082740370999090373516083
11, 1.0000000000100000008274037099909037351608
12, 1.0000000000010000889005823410116136074066
13, 1.0000000000000999200722162640886381268501
14, 1.0000000000000099920072216264088638126850
15, 1.0000000000000011102230246251565404236317
16, 1.0000000000000000000000000000000000000000
17, 1.0000000000000000000000000000000000000000
18, 1.0000000000000000000000000000000000000000
19, 1.0000000000000000000000000000000000000000

of course we hoping to get

 1, 1.1000000000000000000000000000000000000000
 2, 1.0100000000000000000000000000000000000000
 3, 1.0010000000000000000000000000000000000000
 4, 1.0001000000000000000000000000000000000000
 5, 1.0000100000000000000000000000000000000000
 6, 1.0000010000000000000000000000000000000000
 7, 1.0000001000000000000000000000000000000000
 8, 1.0000000100000000000000000000000000000000
 9, 1.0000000010000000000000000000000000000000
10, 1.0000000001000000000000000000000000000000
11, 1.0000000000100000000000000000000000000000
12, 1.0000000000010000000000000000000000000000
13, 1.0000000000001000000000000000000000000000
14, 1.0000000000000100000000000000000000000000
15, 1.0000000000000010000000000000000000000000
16, 1.0000000000000001000000000000000000000000
17, 1.0000000000000000100000000000000000000000
18, 1.0000000000000000010000000000000000000000
19, 1.0000000000000000001000000000000000000000

Any calculation of reals using python, will only be accurate to about the 15th place. 15 decimal places requires 50 bits (here "l" is log)

#  echo "15*l(10)/l(2)" | bc -l
49.82892142331043521840

The IEEE double precision (64 bit) representation of real numbers uses 52 bits (close enough to the 50 bits we see above) to represent the mantissa, 11 bits for the exponent and 1 bit for the sign - see IEEE 754-1985 (http://en.wikipedia.org/wiki/IEEE_floating_point_standard#Double-precision_64_bit). It looks like python is using double precision (64 bit) to represent real numbers.

You start with your number>1, and an estimate of the square root (Many algorithms require you to give it starting value(s). Often these can be almost anything.) The square root is going to be somewhere between the number and 1, so make the estimate the arithmetic mean (the average).

estimate = (number + 1.0)/2.0

You use this estimate to get a better estimate

new_estimate = (estimate + number/estimate)/2.0

How does this work? Let's find the square root of 9.

number = 9
estimate = (9+1)/2.0=5.0
new_estimate = (5 + 9/5))/2.0 = 3.4

We started with estimate=5.0 and get estimate=3.4, which is closer to the answer. The estimate gets closer to the known answer with each iteration (the algorithm is said to converge). When we eventually get to our answer (at least within the precision of a real), there will be no change when we plug in the estimate.

estimate = 3
new_estimate = (3 + 9/3))/2.0 = 3

We keep iterating till the difference between the estimate and then new_estimate is acceptable. (for the proper way to do real comparisons see Lahey compiler) We could use division to test if they are close enough; this will be slow, but will work no matter what size the number is. We could use substraction, which is fast, but the remainder will be large for large numbers and small for small numbers. Ideally we'd like to get the maximum precision possible (double precision), but for the exercise, we'll substract and use a relatively large difference to terminate the algorithm.

Before we write the code that iterates, we need to initialise some numbers

Write some code that does these 4 steps and prints out the value of number and the estimate. Here's my code ^[148] .

We are about to enter a loop. We don't know how many iterations are going to be needed, but we do have a test to show that we're done. Should we use a for loop or a while loop ^[149] ?

At the top of the while loop we test if we should enter the loop again. What's the test that we've found the square root and how do we use it in the while line? Here's my code ^[150]

Assume that the loop is entered (i.e. we don't have the square root yet): Before you entered the loop, when you initialised a few variables, you generated an estimate of the square root and from it a better estimate. You should use exactly the same steps to inside the while loop, to generate a better value of estimate. Heres what happens inside the loop:

Write code to do this, and inside the loop, print out the value for estimate. Here's my code ^[151] .

	Note
	In a for loop, the code for the calculations is all inside the loop. In a while loop, the code for the calculations must also be ahead of the loop (where it runs once), so that the first time the conditional test in the while statement is run, all values will be known.

Where to have the print() statement: For the final code the print() statement will be commented out, so this decision isn't all that important. You can print out values at any step in the while loop and everyone reading your code will know what you've done. The location given is slightly better, as you're printing out the value that caused the loop to be executed. After the loop exits, the last value of new_estimate will be printed by a statement in the following code. If you'd printed after the calculation of new_estimate, then the value for the first iteration of the loop would not be printed.

The only modification left is to print out the final value for the square root ^[152] .

The order of an algorithm e.g. O(n), O(nlogn) describes the rate at which algorithm's execution time increases with the amount of data that it processes. The calculation of the square root has no data (that can be varied), so we can't use this measure. Instead the measure of worth is the increase in running time needed to increase accuracy by one decimal digit (i.e. a factor of 10).

With a print() statement in the loop, you can see the number of iterations before the answer converges (arrives within the acceptable limit).

On Paul Hsieh's Square Root page (http://www.azillionmonkeys.com/qed/sqroot.html) you'll find that you only need 6 iterations to get 53₂ places (a 64 bit double precision real number) i.e. every step gets you 9 bits, a factor of 2⁹=512, closer to your answer. This is pretty fast.

	Note
	End Lesson 20

At this stage we have working code for the Babylonian sqrt(). The algorithm converges quickly (only 6 iterations are required to calculate to an accuracy of 53 bits). However saying that it "converges quickly" doesn't quantify anything (it's a phrase from marketing). You need to be able to say how quickly, using a measurement that is meaningful to just about anybody (e.g. your grandma), not just computer programmers. An obvious first test is to compare the Babylonian sqrt() to the best sqrt() you have (the math module sqrt()). If it's better, we'll do more tests, otherwise we'll just drop it.

To do speed tests on code, we need to be able to measure the execution time of a program. Python has a timer that measures wall clock time in seconds since 1 Jan 1970 (the start of unix time). (see Time access and conversions http://docs.python.org/lib/module-time.html). time() is a (apparently 32-bit) float. On some systems the resolution is only 1sec (programs that run for less than 1 sec will appear to run for 0 time). If the computer is otherwise idle, the wall clock time and the execution time of your program are nearly the same. Here's some code to measure time.

from time import time

n = 1000000

start = time()
#do something n times
finish = time()

print "execution time per iteration = %10.10f", (finish-start)/n

A single calculation of a square root is a little fast for most timers. We usually do a large number of iterations in a loop and then divide by the number of iterations. This brings its own problems: we have to figure out how much time is taken up by the looping code, but we can handle that too (see below). To time your square root code, first make your Babylonian sqrt() into a function. Copy your Babylonian sqrt code to a new file square_root_compare.py

You now have a function. From here on, this function is just a block of code, that you'll do timing runs on. You won't be changing anything inside of it; it's just a black box.

You need a main() to call the function. Write a two-line main() to test your function - assign the variable number the value 1000, and then call babylonian_sqrt() with the parameter number. Test that your code runs.

To do timing runs on the function, main() needs to generate a whole lot different numbers to feed to your the function. Why don't you just feed the same number over and over? Compilers and interpreters are supposed to be smart. If they see you do the same operation over and over on the same data, the compiler/interpreter will recognise what's going on and will just return the same answer, without going through the bother of calculating it over and over. You don't know if your python interpreter is this smart (how would you find out ^[153] ), but you should expect that one day you'll bump into a smart compiler/interpreter, and your timing runs will be meaningless. You may as well start writing your benchmarking code correctly from the start.

One way of producing a whole lot of different numbers is to use all the numbers produced by range(1,largest_number), where largest_number is say 10,000.

You now have a block of code in main() that calls the function largest_number of times. You now want to time that block

Here's my code ^[155] . This outputs the time required to calculate the square roots of all the numbers in the range 1..1000. This piece of timing code is a self contained block. You aren't going to change it (except for the value of largest_number)

The average time for the sqrt() calculation depends on the number of times you looped (it shouldn't, but it does). You'll find that the time/iteration depends on whether you did it 10,100 or 1,000,000 times. You need to find a range of largest_number for which the time/iteration is the fastest availalbe and relatively constant. To handle this, you need to call the timing loop with a range of values of largest_number. To do this, put the timing loop above, inside another loop, which feeds different values of largest_number into the timing loop. Putting one loop inside another is called "nesting loops".

To do this, you'll need to use a list. Go to the section explaining lists and then return here.

Create a list named largest_numbers[] holding the values 1,10,100....1000000. Read these values out, using a loop (for or while?) assigning the values one at a time to largest_number. The timing loop, being nested inside the new outer loop, will now have to be indented one step to allow it to be parsed correctly. Here's what the nested loops look like

largest_numbers=[1,10,100,1000,10000,100000,1000000]	#list of values for outer loop
for largest_number in largest_numbers:		#new outer loop
	for number in range(0,largest_number):	#original loop
	        babylonian_sqrt(number)

Here's my code ^[156] .

Here's my output

#  ./square_root_compare.py
largest number          1 time/iteration 0.0001909
largest number         10 time/iteration 0.0000348
largest number        100 time/iteration 0.0000295
largest number       1000 time/iteration 0.0000341
largest number      10000 time/iteration 0.0000398
largest number     100000 time/iteration 0.0000474
largest number    1000000 time/iteration 0.0000512

Look to see if there's a range of values for largest_number where time for doing a sqrt is independant of the loop size (it should be the fastest time). For largest_number having a value of 10 or more, the Babylonian square root function takes 30-50usec on my 1GHz machine. One one student's machine (MacOS), the time/interation kept dropping with in the range values for largest_number above, so I sent him off to do runs with increasing values of largest_number. Another student's machine (WinXP) reported 0 time for the first 3 runs, and presumably has a time() function with a resolution of only 1 sec.

	Note
	End Lesson 21

Next we want to compare the speed of the Babylonian sqrt() to the built-in sqrt() (in the math module).

Here's the new calling code in main() ^[157] .

Calling range() and setting up the looping consumes cpu cycles, so we need to subtract that time from our results. Code up a 3rd timing loop that doesn't do anything, so later you can subtract that time out. Here's a loop that does nothing (you can't have an empty line, the parser gets confused).

	for number in range(0,largest_number):
		pass

Here's my code ^[158] . Here's the output

# ./square_root_compare.py
babylonian_sqrt: largest number          1 time/iteration 0.0001969337
library_sqrt:    largest number          1 time/iteration 0.0000259876
empty_loop:      largest number          1 time/iteration 0.0000119209

babylonian_sqrt: largest number         10 time/iteration 0.0000346184
library_sqrt:    largest number         10 time/iteration 0.0000043154
empty_loop:      largest number         10 time/iteration 0.0000017881

babylonian_sqrt: largest number        100 time/iteration 0.0000295210
library_sqrt:    largest number        100 time/iteration 0.0000031996
empty_loop:      largest number        100 time/iteration 0.0000008512

babylonian_sqrt: largest number       1000 time/iteration 0.0000341501
library_sqrt:    largest number       1000 time/iteration 0.0000030880
empty_loop:      largest number       1000 time/iteration 0.0000007720

babylonian_sqrt: largest number      10000 time/iteration 0.0000398650
library_sqrt:    largest number      10000 time/iteration 0.0000031442
empty_loop:      largest number      10000 time/iteration 0.0000007799

babylonian_sqrt: largest number     100000 time/iteration 0.0000459829
library_sqrt:    largest number     100000 time/iteration 0.0000033149
empty_loop:      largest number     100000 time/iteration 0.0000009894

babylonian_sqrt: largest number    1000000 time/iteration 0.0000511363
library_sqrt:    largest number    1000000 time/iteration 0.0000033719
empty_loop:      largest number    1000000 time/iteration 0.0000010112

The time for the empty loop is 0.7-1.0usec. Here's the timing results, after subtracting the empty loop time.

The library sqrt() is 15-25 times faster than the Babylonian sqrt(). Even though the babylonian_sqrt() only needs 5 iterations to get the answer, the library routine gets the answer in the time of 1/3-1/5th of a loop.

Python being an interpretted language will be slower than a compiled language. As well, python is object oriented, making it slower yet. C is one of the languages of choice when speed is required. Here's the same program as above (comparing the Babylonian and math library sqrt()) written in C (it uses the low resolution timers). With your current programming knowledge, you should be able to read this code and know what it's doing. ^[159] . Here's the output

time/iteration: babylonian         10000000  0.00000122000000000000
time/iteration: library            10000000  0.00000010800000000000
time/iteration: empty              10000000  0.00000000700000000000

time/iteration: babylonian        100000000  0.00000135340000000000
time/iteration: library           100000000  0.00000010790000000000
time/iteration: empty             100000000  0.00000000680000000000

time/iteration: babylonian       1000000000  0.00000147792000000000
time/iteration: library          1000000000  0.00000010824000000000
time/iteration: empty            1000000000  0.00000000691000000000

Here's the comparison of the timing results, subtracting 7nsec for the empty loop (compared to 700nsec for the empty loop in python).

C is 20-40 times faster than python, at least for doing sqrt(). So why does anyone code in python? If speed is your main requirement you don't. However it takes longer to write C than it does to write python and people who don't program a lot, can write working python without too much trouble, but may not be able to get a C program to work. If you only require a few square roots and it takes 5 mins to write it in python and 10 mins to write it in C, then you do it in python. If you need to do 10¹² square roots, then you do it in C.

If your code only takes 1 min to run and you're only going to run it a couple of times, then you really don't care if it takes 10 times longer. However if you're going to be running the 1 minute program hundreds of times, or your program will take a week to run, then a factor of 10 in speed is important.

	Note
	End Lesson 22

You've just completed a piece of code and have a nice story to tell. For the rest of your life, whether you're a coder or something else, you're going to have to sell what you've done to someone else. Here you have a piece of code.

It's simple to tell other coders of your work. You just walk them through the code. They'll just go "yup, yup... yup. Nice story." and go back to their work. A little harder is a technically trained person, who isn't neccessarily a python coder.

Here's what you do.

This is usually phrased as

It's a conceptually simple story with a conclusion that everyone will agree on. A presentation on this piece of code shouldn't take more than 5-10mins.

	Note
	End Lesson 23

Let's look at the part of the circle in the first quadrant. This diagram shows details for the point (x,y)=(0.6,0.8).

	Note
	The python code for the diagrams used in this section is here [160] . Since this code will only be run once, there's no attempt to make it fast.

Figure 4. Pythagorean formula for circumference of a circle

Pythagorean formula for circumference of a circle

From Pythagorus' theorem, we know that the distance from the center at (0,0), to any point (x,y) on the circumference of a circle of radius=1, is 1. The square of the distance to the center of the circle (the hypoteneus) is equal to the sum of the squares of the two sides, the lengths of which are given by (x,y). Thus we know that the locus of the circumference of a circle of radius 1 is

x*x + y*y = 1*1

Note

locus: the path traced out by a point moving according to a law. Here's some examples. The circumference of a circle is the locus of a point which moves at constant distance from the center. A straight line is the locus of a point that moves along the path, which is the shortest distance between two points. An ellipse is the locus of a point which moves so that the sum of the distances to two points (the two focii) is constant. (The planets move around the sun in paths which are ellipses. The sun is one of the focii.) Note etymology: focus == fireplace

In the diagram above if x=0.6, then y=0.8. Let's confirm Pythagorus.

x^2=0.36, y^2=0.64; x^2+y^2=1.

Rearranging the formula for the (locus of the) circumference, gives y (the height or ordinate) for any x (the abscissa).

x^2 + y^2 = 1^2
y^2 = 1-x^2
y = sqrt(1-x*x)

i.e. for any x on the circumference, we can calculate y. e.g. x=0.6: what is y?

y = sqrt(1-x*x)
  = sqrt(1-0.36)
  = sqrt(0.64)
  = 0.8

If x=0.707, what is y ^[161] .

The point on the circle with y=0.5: what angle (by eyeball, or by marking off lengths around the circumference) does the line joining this point to the center, make with the x axis ^[162] .

Similarly we can calculate x knowing y.

x = sqrt(1-y*y)

The value of π is known to many more decimal places than we can ever use on a 32 bit computer, so there's no need to calculate it again. Instead let's assume that either we don't know its value, or that we want to do the numerical integration on a problem whose answer is known, to check that we understand numerical integration.

The area of a circle is A=pi*r², giving the area of the quadrant as pi/4. Knowing the locus of the circumference (we have a formula that gives y for any x), we will numerically calculate the area of the circle, thus giving us an estimate of the value of π.

If we divide the x axis into say 10 intervals, we can make these intervals the base of columns, whose tops intersect the highest part of the circle in that slice (which is on the left side of each slice).

Figure 5. Upper Bound of area under slices of a circle

Upper Bound of area under slices of a circle

When we add these slices together, we'll get an area that is greater than π; i.e. we will have calculated an upper bound for π.

Here's the set of columns that intersect the lowest part of the circle in each interval (here the lowest part of the circle is on the right side of the slice).

Figure 6. Lower Bounds of area under slices of a circle

Lower Bounds of area under slices of a circle

We can calculate the area of the two sets of columns. If we sum the sets of upper bound columns, we'll get an estimate which is guaranteed to be more than pi/4 and for the set of lower bound colums, a number guaranteed to be less than pi/4.

Note

We could have picked the point in the middle of the interval to calculate the area. The answer will be more accurate, but now we don't know how accurate (we don't even know if it's more or less than π). The advantage of the method used here is that we have an upper and lower bound for π, and so we know that the value of π is in this range. We could have used tighter bounds - lower bound by constructing a straight line joining the left and right end of each interval (giving a trapezoid), - upper bound by making a line tangent to the circle. (This is more bother than it's worth, and are left as an exercise for the student.)

If we progressively decrease the size of the interval (by changing from 10 intervals, to 100 intervals, to 1000 intervals..) the approximation to a circle by the columns will get better and better giving us progressively better estimates of pi/4. Here's the diagram with 100 slices.

Figure 7. Lower Bounds of area under 100 slices of a circle

Lower Bounds of area under 100 slices of a circle

We have a problem similar to the fencepost error problem: how many heights do we need to calculate if we want to calculate both the upper and lower bounds for 10 intervals. ^[163] ?

Here's code to generate the heights of the required number of columns (you'll need sqrt(), from the math module). I ran this code in immediate mode, but you can write a file called slice_heights.py if you like.

>>> from math import sqrt
>>> for x in range(0,11):
...     print x, sqrt(1-1.0*x*x/100)
... 
0 1.0
1 0.994987437107
2 0.979795897113
3 0.953939201417
4 0.916515138991
5 0.866025403784
6 0.8
7 0.714142842854
8 0.6
9 0.435889894354
10 0.0

Why did I use 11 as the 2nd argument to range() ^[164] ? What is the "1.0" doing in "1.0*x*x" ^[165] ?

	Note
	End Lesson 24

Note

Python is a scripting language, rather than a general purpose language. Python can only use integers for loop variables. To feed real values to a loop, python requires a construct like #num_intervals, interval are integers num_intervals=1000 for i in range(1,num_intervals): real_number=i*1.0/num_intervals function(real_number) where x,num_intervals make the real real_number. It's not immediately obvious that the calculations are ranging over values 0.0..1.0 (or more likely start..end). In most languages, real numbers can be used as loop variables, and can use the construct #interval, start, end, x are reals interval = 0.001 for (x = start; x < end; x += interval) function(x) Here it's clear that x is a real in the range start..end.

Write code pi_lower.py to calculate the lower bound for π using 10 intervals. Use variable names num_intervals, height, area for the number of intervals, height of each column, and for the cumulative area. Start with a loop just calculating the height of each slice, printing the loop index and the height each time. For the print statement use something like

I've been using the single letter variable 'x' as the loop parameter thus far, since it's the loop parameter conveying the position on the x-axis. However loop parameters (in python) are integers, while the x position is a real (from 0.0..1.0). Usually loop parameters (which are integers) are just counters and are given simple variable names 'i,j,k,l...', which are easy to spot in code as they are so short. In this case the loop parameter is the number of the interval (which goes from 0..num_intervals). I tried writing the code with the name "interval" mixed in with "num_intervals" and it was hard to read. Instead I will use 'i'.

print "%d, %3.5f" %(x, height)

When that's working, calculate the area of each slice and add it to the variable area. The area of each rectangle is height*base. You've just calculated height; what is the width of the base in terms of variables already in the code ^[166] ?

At the end, print out the lower bound for π with a line like

print "lower bound of pi %3.5f" %(area*4)

Here's my code for pi_lower.py ^[167] and here's my output.

0, 0.99499, 0.09950
1, 0.97980, 0.19748
2, 0.95394, 0.29287
3, 0.91652, 0.38452
4, 0.86603, 0.47113
5, 0.80000, 0.55113
6, 0.71414, 0.62254
7, 0.60000, 0.68254
8, 0.43589, 0.72613
9, 0.00000, 0.72613
lower bound of pi 2.90452

Do the same for the upper bound of π writing a file pi_upper.py and using variables num_intervals, Height, Area (note variable names for the upper bounds start with uppercase, to differentiate them from the lower bounds variables). Here's my code ^[168] and here's my output.

0, 1.00000, 0.10000
1, 0.99499, 0.19950
2, 0.97980, 0.29748
3, 0.95394, 0.39287
4, 0.91652, 0.48452
5, 0.86603, 0.57113
6, 0.80000, 0.65113
7, 0.71414, 0.72254
8, 0.60000, 0.78254
9, 0.43589, 0.82613
upper bound for pi 3.30452

The two bounds (upper and lower, from the output of the two programs) show 2.9<pi<3.3 which agrees with the known value of π.

The two pieces of code look quite similar. Also note some of the numbers in the outputs are the same (how many are the same ^[169] ?) We should check for code duplication (in case we only need one piece of code).

Figure 8. Upper and Lower Bounds of area under a circle

Upper and Lower Bounds of area under a circle

Looking at the diagram above which shows upper and lower bounds together, we see the following

Redrawing the diagram, shifting the lower bounds slices one interval to the right shows

Figure 9. Upper and shifted Lower Bounds of area under a circle

Upper and shifted Lower Bounds of area under a circle

This shows that the upper and lower bounds only differ by the area of the left most slice. This means only one loop of code is needed to calculate both bounds. Look in the output of pi_lower.py and pi_upper.py for the 9 numbers in common.

The duplication arises because the lower bound for one interval is the upper bound for the next interval and we only need to calculate it once. The first interval for the upper bound and the last interval for the lower bound are unique to each bound and will have to be calculated separately.

	Note
	This is a general phenomenon: When calculating a value by determining an upper and lower bound, if the curve is monotonic, you should expect to find values that are used for both the upper and lower bound.

Write code (call the file pi_2.py - there is no pi.py anymore) to calculate the area common to both bounds (i.e. except for the two end pieces) in one loop. Use x for the loop control variable (the slice number), h for the height of each slice and a to accumulate the common area. In the loop output the variables for each iteration with a line like

	print "%10d %10.20f %10.20f" %(x, h, a)

After exiting the loop, add the two end areas, one for the lower bound and one for the upper bound to give area and Area and output your estimates of the lower and upper bounds for π with a line like

print "pi upper bound %10.20f, lower bound %10.20f" %(Area*4, area*4)

Here's the code ^[170] and here's the output

pip:# ./pi_2.py
	 1 0.99498743710661996520 0.09949874371066200207
	 2 0.97979589711327119694 0.19747833342198911621
	 3 0.95393920141694565906 0.29287225356368368212
	 4 0.91651513899116798800 0.38452376746280048092
	 5 0.86602540378443859659 0.47112630784124431838
	 6 0.80000000000000004441 0.55112630784124427841
	 7 0.71414284285428497601 0.62254059212667278711
	 8 0.59999999999999997780 0.68254059212667272938
	 9 0.43588989435406727546 0.72612958156207940696
pi upper bound 3.304518, lower bound 2.904518

giving the same result: 2.9 < pi < 3.3

To get a more accurate result, we now increase the number of intervals, thus making the slices smaller and reducing the jaggyness of the approximation to the circle.

	Note
	End Lesson 25

Note

Engineers take pride in being able to do what is called a back of the envelope calculation. While it may take weeks of computer time to get an accurate answer, an engineer should be able to get an approximate answer using only the back of an envelope to do the calculation (extra points are added if you can do it in the elevator between floors). As part of the calculation, you also need to be able determine the accuracy of your answer (is it correct within a factor of 100,10,2, 10%?) or state whether it's an upper or lower bound. If the answer is within reasonable bounds, then it's worth spending a few months and a lot of money to find the real answer. e.g. how long does it take to get from Earth to Mars using the minimum amount of fuel. It's going to be somewhere between the half the time of Earth's orbit and half the time of Mar's orbit (both planets orbit using no fuel at all), i.e. between 6 months and a year - let's make it 9 months. The answer is 8.6 months (Earth-Mars Hohmann trajectory http://www.newmars.com/wiki/index.php/Hohmann_trajectory). Note The exact answers, the Hohmann transfer orbit time between planets (see Gravitational Slingshot http://en.wikipedia.org/wiki/Gravitational_slingshot) were in Clarke's tables (along with log and trig tables, and of course the properties of elements) which all high school students in my time (1960's) carried to science class. I remember looking for the element with the highest density (which I remember as being Osmium).

Before we change the number of intervals to some large number (like 10⁹), we need some idea of the time this will take. We could change intervals to 10,100... and measure the time for the runs (the quite valid experimental approach) and extrapolate to 10⁹ intervals. Another approach is to do a back of the envelope calculation of the amount of time we'll need. We don't need it to be accurate - we just need to know whether the run will take a second, an hour or a year, to see if it's feasible to run the code with intervals=10⁹. Lets say we use 1G intervals. The calculation goes like this

A 1GHz computer then will take about 100 secs for num_intervals=10⁹. If you have a 128MHz computer, expect the run to be 8 times longer (800secs=13mins). In this case, you may not be able to do a run with num_intervals=10⁹ in class time, but you should be able to do a run with num_intervals=10⁸.

Do multiple runs of pi_2.py increasing num_intervals by a factor of 10-100 each time, noticing the changes in the upper and lower bounds for π. (Comment out the print() statement inside the loop or it will take forever.)

	Note
	This is a problem with how I wrote the python code. It has nothing to do with numerical integration, but you're going to run into problems like this and you need to know how to get out of them. So we're going to take a little detour to figure this one out.

For a large enough value of num_intervals (try 10⁸ or 10⁹, depending on your computer), the program exits immediately with a MemoryError: the machine doesn't have enough memory to run the program. When you have a problem with code, that you've thought about for a while and can't solve, you don't just dump the whole program (with thousands of lines of code) in someone's lap and ask them what's wrong. Instead you pare down the code to the simplest piece of code that will demonstrate the problem and then ask them to look it over. Here's a simplified version of the problem showing the immediate exit:

>>> num_intervals=100000000
>>> height=0
>>> for interval in range(0,num_intervals):
...     height+=interval
... 
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
MemoryError

	Note
	If I'd gone one step further to >>> for interval in range (0,100000000): ... print interval ... Traceback (most recent call last): File "<stdin>", line 1, in ? MemoryError I probably would have asked myself why range failed for a large enough argument, and figured it out myself. But I didn't and I had to go ask for advice on a python mailing list.

With a smaller value for num_intervals (try a factor of 10 smaller), this code will now run, but take up much of the machine's memory (to see memory usage, run the program top). I thought that the code was only creating a handful of variables (x, height, num_intervals). In fact range() creates a list (which I should have known) with num_intervals number of values of x, using up all your memory. In all other languages, the code equivalent to

for i in range(0,num_intervals):

calculates one new number each iteration, and not a list at the start, so you don't have this problem. You don't need a list, you only need to create one value of x for each iteration.

While I was waiting for an answer on the python mailing list, I changed the for loop to a while loop. Here's the equivalent while loop

>>> num_intervals=100000000
>>> height=0
>>> interval=0
>>> while (interval<num_intervals):
...     height+=interval
...     interval+=1
... 

which works fine for large numbers of iterations (there doesn't need to be anything in particular inside the loop to demonstrate that a while loop can handle a large number of iterations).

Back to range(): Not only do you run into problems if you create a list longer than your memory can hold, but even if you have infinite memory, there is another upper limit to range(): you can only create a list of length=2³¹

>>> interval=range(0,10000000000)	#list 10^10
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
OverflowError: range() result has too many items

>>> interval=range(0,10000000000,4)	#list of every 4th number of 10^10 numbers = 2.5*10^9
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
OverflowError: range() result has too many items
>>> 

>>> interval=range(0,10000000000,5)	#list of every 5th number of 10^10 numbers = 2*10^9
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
MemoryError

Why did the error change from OverflowError with a list of length 2.5*10⁹ to MemoryError with a list of length 2.0*10⁹?

OverflowError indicates that the size of the primitive datatype has been exceeded. MemoryError indicates that the machine doesn't have enough memory to do that operation (i.e. your machine has less than 2G memory) - python will allow a list of this length, but your machine doesn't have enough memory to allocate a list of this size.

What primitive data type is being used to hold the length of lists in python ^[171] ?

If the transition from MemoryError to OverflowError had happened between 4.0*10⁹ and 4.5*10⁹, what primitive datatype would python have been using ^[172] ?

It turns out on a 64 bit system you can still only create a list of length 2³¹.

The python way of overcoming the memory problem of range() is to use xrange() which produces only one value of the range at a time (xrange() produces objects, one at a time, while range() produces a list). Change your pi_2.py to use xrange() and rerun it with a large value of num_intervals. Here's my fixed version of pi_2.py called pi_3.py. The only change to the code is changing range() to xrange() ^[173] .

	Note
	Use xrange() whenever you are creating a list of size comparable to the memory in your machine.

When coding your goals (in order) are

The current code runs correctly, but it runs slowly. Your first attempt at coding up a problem is always like this: you (and everyone else involved) wants to see if the problem can be handled at all. The process of changing code so that it still runs correctly, but now runs quickly, is called "optimising the code".

If you have a piece of software running on a $1M cluster of computers and you can speed it up by 10%, you just saved your business $100k. The business can afford to pay a programmer $100k for a year's work to speed up the program by only 10%. Optimising code is tedious and difficult and you can't always tell ahead of time how successful you'll be. People who optimise code can get paid a lot of money. On the other hand, the people with the money can't tell whether the code needs to be optimised (they think if it works at all, then it's finished) and can't tell a person who can optimise, from one who can't.

If you're selling the software to someone else, then the cost of the extra hardware, needed to run the unoptimised software, is borne by the purchaser and not you, so software companies (unless they have competition) have no incentive to optimise their code. Software companies, just like many businesses, will externalise their costs whenever they can (see Externality http://en.wikipedia.org/wiki/Externality, and Cost_externalizing http://en.wikipedia.org/wiki/Cost_externalizing). With little competition in the software world, there's a lot of unoptimised and badly written code in circulation.

The place to look for optimising, is the lines of code where most of the execution time is spent. If you double the speed of the code where only 1% of the execution time occurs, the code will now take 99.5% of the time, and no-one will notice the difference. It's better to tackle the part of the code where 99% of the time is spent. Doubling the speed there will halve the run time. Places in the code where a lot of time is spent are loops that are run many times (or millions of times). In pi_3.py the program is only a loop and it's running 10⁹ times.

Let's look at the loop in pi_3.py.

Here's my test_optimisation.py ^[178] and here's the output run on a 200MHz machine

pip:# ./test_optimisation.py
	                           iterations          sum   iterations/sec
unoptimised code:                     1000000 0.0000010000    38134.77894
precalculate num_intervals^2:         1000000 0.0000010000    63311.38723
multiple by reciprocal:               1000000 0.0000010000    88156.25065

showing a 2.5-fold increase in speed.

	Note
	End Lesson 26

Don't alter your pi_3.py; you aren't finished with test_optimisations.py yet.

Computers have a large but limited range of numbers they can manipulate. We're calculating π, whose value lies well in the range of numbers that computers can represent. However numerical methods, like we are using here, use values that are very small (and the squares of these numbers are even smaller). Other codes can divide small numbers by large numbers, large numbers by small numbers, square numbers that are smaller than the square root of the smallest representable numbers, or square numbers that the are bigger than the square root of the biggest representable numbers. You'll wind up with 0 instead of a small number or NaN (not a number, the computer result of dividing by 0). If you're close to the edge there, you drop significant figures and wind up with an inaccurate answer. Your rocket may not blow up for 10 launches, then one day when the breeze is coming from a different direction, the guidance system won't be able to correct and your rocket will blow up.

With the range of numbers we're using in the numerical integration, we're not going to have this problem. However you're writing code that for a different range of numbers will have this problem. Someone could swipe your code and use it for a new range of inputs, or you could re-use it in another piece of code. You want to write the code to minimise problems in the distant future. Here's where this problem shows up in our code

(1.0*x*x)/(num_intervals*num_intervals)

Depending on x, this is a small number squared, over a large number squared (x*x might evaluate to 0.0 and the division by num_intervals*num_intervals would be irrelevant; or num_intervals might be large, causing num_intervals*num_intervals to overflow i.e. NaN) - a disaster with your rocket's name on it. Instead you should write

(1.0*x/num_intervals)^2

Numbers and their scaling constants should not be separated. In computing, dividing a data value by some appropriate scaling constant is called normalisation. The word "normalisation" is used in different places to mean different things (unfortunately). Where else have we seen the term ^[179] ? You can take "normalised" to mean "anything that fixes up numbers so that something awful doesn't happen". What that may mean in any particular situation, you'll have to know.

The term used in science/engineering is "reduced" (rather than "normalised"), where the reduced value is dimensionless. Say a machine/system can do a maximum of 1000 operations/hr. If it's doing 500 operations/hr, its reduced speed is 0.5 (not 0.5 of its maximum speed, just 0.5 i.e. a number). Rather than refer to the speed of a plane (which has the dimensions of l/t), you can refer to its Mach number (the speed of the plane/speed of sound, this is dimensionless). The speed of sound depends on temperature. The temperature of air varies with altitude, so the speed of sound varies with altitude. Aerodynamics, especially near the speed of sound depends on the Mach number and not the speed (you'll need your speed to know when you'll arrive at your destination). If you're studying the aerodynamics of a plane at different altitudes, you want to know the Mach number.

In computations, often there will be dimensionless terms like x/y (which might be speed/speed of sound). In this case while each of the variables x,y might have a large range of values (e.g. 10**-160 to 10**160), the range of x/y might be from 10**-1 to 10**1.

Let's say we have a calculation which requires the value of (x/y)² when each of x,y have large ranges, but x/y has a small range. Here's how the ant climbing a hill, (floating point precision) looking for the highest point, gets stuck and starts wandering around aimlessly.

#small numbers (you have to guard against this possibility when doing optimisation problems)
>>> x=10**-161
>>> x*x
9.8813129168249309e-323
>>> x=10**-162
>>> x*x					#underflow: x*x is too small to be represented by a 64-bit real
0.0
>>> interval=10**-150
>>> (x*x)/(num_interval*num_interval)	#underflow: invalid answer
0.0
>>> (x/num_interval)*(x/num_interval)	#normalised number gives valid answer
9.9999999999999992e-25			#the correct answer is 10^-24, this answer is near enough

#big numbers (this doesn't happen a real lot)	
>>> num_interval=1.01*10**154
>>> num_interval*num_interval
1.0201e+308
>>> num_interval=1.01*10**155		#if I don't have the 1.01, the number will be an integer
>>> 					#python can have abitrarily large integers 
>>>					#and you'll get the correct answer
>>> num_interval*num_interval		#squaring a large number	
inf					#overflow
>>> x=1.01*10**150			#squaring a slightly smaller number
>>> x*x
1.0201000000000001e+300
>>> (x*x)/(num_interval*num_interval)	#unnormalised division
0.0					#wrong answer
>>> (x/num_interval)*(x/num_interval)	#normalised division
1.0000000000000002e-10			#correct answer

Whenever you have

y=(a*b*c)/(e*f*g)

figure out which denominator terms normalise which numerator term (you'll know from an understanding of the physics of the problem). You'll wind up changing the code to something like this

y=(a/e)*(b/f)*(c/g)

With the numbers we're using in the numerical integration, we aren't getting any benefit from using reduced numbers. None of the numbers under or overflow on squaring. However reduced numbers should always be used, just as seat belts should be used in cars. Here's the range of numbers involved

We're dealing with numbers 10^-18..1.0 whether we reduce or not.

Here's a case when reduced variables helps.

Here if we square x,y separately, we'll get an under or overflow. If we use reduced numbers, we'll get a correct answer.

Why did I get you to use reduced numbers when you don't need them? It's a good programming practice. You never know what someone else is going to do with your code. When another coder reads your code, he can tell that you know what you're doing and he won't have to have to look for stupid coding errors. This sort of code gives you good programming cred to people reading your code. You need to fix problems before they happen (not after). This country has blown up two Space Shuttles because people didn't fix problems when they came up to be fixed. You don't want to write code that blows up Space Shuttles.

In computers, division is done by long division, the same way you divide by hand. Long division is slow compared to multiplication, which on a computer (with special hardware) can be done in two steps. In a loop, if you divide by a constant, it's faster to multiply by the reciprocal.

#loop, slow way
	a=h/num_intervals

#faster way
interval=1.0/num_intervals	#outside loop
#loop
	a=h*interval

Add another section with a loop to test_optimisations.py using the reduced number x/num_intervals (since 1/num_intervals is used several times, you should multiply by the reciprocal, and so the code will actually use x*interval). Here's the code ^[180] and here's the ouput.

./test_optimisation.py
	                           iterations          sum   iterations/sec
unoptimised code:                     1000000 0.0000010000    37829.84246
precalculate num_intervals^2:         1000000 0.0000010000    63759.14714
multiply by reciprocal:               1000000 0.0000010000    99322.31288
normalised reciprocal:                1000000 0.0000010000    60113.31327

Here's my final version of test_optimisations.py ^[181] .

Because you're reducing numbers, you can no longer take advantage of the precalculated square, and your code is slower. In general, you don't care if normalisation slows down the code; you'd rather your rocket's navigation system worked correctly but slowly, than crashing the rocket quickly. If you're doing the runs yourself, and you know that the data is well conditioned and doesn't have to be reduced, you can do a run that takes a week instead of 10days. However, in 6 month's time, if your rocket blows up with someone's $1B satellite on board, no-one's going to be impressed when they find that your code didn't normalise. If your code will be released to a trusting and unsuspecting world, to people who will feed it any data they want, you must normalise.

The optimisations we tested were

If the loop had been longer (say 50 lines long), with lots of other calculations, then moving the generation of the constant out of the loop may not have made much difference. However the optimisations shown here should be done as a matter of course; they will assure anyone reading your code, that you know what you're doing.

The optimisations shown here are relatively routine and are regarded as normal programming practice, at least for a piece of code that takes a significant fraction of the runtime (e.g. a loop which is run many times). You wouldn't bother optimising if the lines of code were only run a few times. The unoptimised code that I've presented to you was deliberately written to be poor code so that you could make it better.

Note

End Lesson 27: At this stage I'd left the students with the results of running test_optimisation.py, expecting them to figure out on their own, the optimisations that would be incorporated into the numerical integration. By the next class, they couldn't even remember the purpose of reducing numbers, so I went over this lesson again, and had them disect out the optimisation that would be transferred to the numerical integration.

Which optimisations are you going to use in the numerical integration? We have tried these loops:

sum = 0

#unoptimised
for i in xrange(0,num_intervals):
     sum+=1.0/(num_intervals*num_intervals)

#faster
#precalculate constants (here num_intervals^2)
num_intervals_squared=num_intervals*num_intervals 
for i in xrange(0,num_intervals):
     sum+=1.0/num_intervals_squared
 
#fastest
#multiply by constants, don't divide
interval_squared=1.0/num_intervals_squared
for i in xrange(0,num_intervals):
     sum+=1.0*interval_squared

#safe, about the speed of "faster", not "fastest"
#use reduced numbers 
interval=1.0/num_intervals
for i in xrange(0,num_intervals):
     sum+=(1.0*interval)**2

The optimisations we can use for the final code are

We must use reduced numbers for safety. Once we've made this choice, the only other optimisation left is to multiply by the inverse of num_intervals. We can't use any of the other optimisations. We don't get the best speed, but we will get the correct answer for a greater range of x, num_intervals

As a programmer, if you're selling your code, you are in somewhat of a bind here. If your code is compared to a competitor's, which uses all the speed optimisations, and which doesn't reduce their numbers, your code will run at half the speed of the competition. It would take quite a sophisticated customer (one who can read the code) to understand why your code is better. In my experience, people who want the program, rarely know the back from the front of a computer, much less understand the trade-offs involved in coding for safety or speed. If you tell the customer that you can make your code run fast too, but it just won't be safe, and they buy your code and it does blow up a rocket or kill a patient in a hospital, they'll still take you to court for shoddy programming. The "poor unknowing customer", who bought the code knowing full well that wasn't safe, won't be blamed.

Solaris (the unix variant used as the OS on Sun's computers) is slow compared to their competitor's OSs and it's disparaged by calling it "Slolaris". It throws up a lot of errors in the logs (say when there's a network problem), while the logs of a competitor's machine, on the same network, seeing the same problems, are empty. People in the know defend Solaris saying it is slow because it's safer and that the other OSs are seeing errors too, but aren't reporting them. The people on the other machines say "my machine isn't having any problems, Slolaris is having the problem!". Expect this if you write good code.

Here's the way out of the bind.

We now need to incorporate the optimisations, and time our π calculation runs. Code takes a while to load from disk. If it's just been run, a subsequent run will load faster (the disk cache knows where the file is or will still have the file cached in memory). Speed then will be affected by how long since the code was last run (after about 15 secs, the cache has been flushed). When you do several runs, you may notice that the first run is slow.

Copy pi_3.py to pi_4.py and include the optimisations that we found to be useful for this code (hint ^[182] ). Here's my version of pi_4.py ^[183] .

	Note
	This section hasn't been presented in class. It was written while the kids are writing their presentations.

Once we have the upper and lower bounds, the value of π comes by comparing the two numbers and accepting from the left, the numbers that match. As soon as we find a non-match, we discard the remaining digits. Since we're calculating π only once, we could match by hand and write down the answer. However your high school wants you to know array operations, so we'll use array operations to do the matching. In python, arrays are implemented with lists. If you need to, first refresh your memory of list operations lists and list of lists. Then go to arrays. Here you'll learn about arrays (you don't need to know all the material in this section to do the following section here, but your school wants you to understand arrays, and it's a good opportunity to get aquainted with them). Then go to the following section turning reals into strings and return here.

Lets do the match by finding the char elements of a string. (I wish to thank members of the North Carolina Zope/Python User Group for help here. There is a neater solution using the python specific zip() which I don't want to use - I'm only using code that will work in any language.) Write code compare_reals.py that does the following

Should you output a string or a real? With the upper and lower bounds being reals, you might expect to write code returning π as a real. However if the string was "3.141", after conversion to a real it might be 3.140999999999995. With correct formatting in the print statement, it might be possible to output this as 3.141, but you would have to know ahead of time the number of places to output. It's pointless to do this when you already have the correct number of places in the string form of π.

	Note
	Optional: If you wanted to output a real for pi_result. If you were to output (or return) a real, you'd use the function float() on the string. Test your code when no digits match (what is float("")?). Will your rocket crash and burn [185] ?

Here's my code ^[186] . What is the order of the algorithm? i.e. if the length of the strings are increased by a factor of n, will the runtime increase by n times or n² times (or some other function of n) ^[187] ?

The value of π is truncated to the number of digits presented. It would be more accurate if we rounded the number up or down. Can we do that? If the first non-matching digits are both 5 or above, we can definitely round up. If the first non-matching digits are both 4 or below, we can definitely round down. What if one digit us 4 and the other 7? We can't say anything. We can't have a scheme where the last digit may be rounded or may be truncated, but we can only tell by looking up the original computer print out. Our only choice is to say that the value of π is truncated to the quoted number of places.

Copy the code to compare_reals_2.py. Make the code into a function string longest_common_subreal(real,real) (i.e. takes parameters of two reals, and returns a string). Change the names of variables to be more generic i.e. string_lower becomes string_1, so that the code will apply to any situation, not just to upper and lower bounds. In main() call the function with two reals (use typical values for upper and lower bounds) and print out the result. Here's my code ^[188] .

Copy pi_4.py to pi_5.py. Add longest_common_subreal() to the code, so that it outputs the truncated string for the estimate for π. Here's my code ^[189] .

Using pi_5.py pick a high enough value for num_iterations to get meaningful output (it looks like I used num_iterations=10⁹, you won't have time to run this in class). Do runs starting at num_intervals=1 (for Cygwin under Windows, the timer is coarse and will give 0 for num_intervals<10000), noting the upper and lower bounds for π. We want to see how fast the algorithm converges and get progressively tighter bounds for π. Here's my results as a function of num_intervals.

	Note
	End Lesson 28