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

The number system in common daily use is called decimal. The base of the decimal system is 10. The decimal system is used because we have 10 fingers. The hardware in a computer only has two states (0/1), so a computer uses the binary system. There are no decimal numbers inside a computer. Binary numbers are transformed by software to a decimal representation on the screen for us. Number systems used in computing are

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

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.

Let's figure out formally what the decimal number 102 represents. This exercise is trivial, but it will be used next for the binary system. We assemble a list of the powers of the base of the system (the base here is 10, and powers of the base are 1=10^0, 10=10^1, 100=10^2...).Then starting with the smallest number in our list 1 (then 10, then 100...), we look for the first power of 10 that's bigger than our number 102. This is 1000. Then we go back by one member in our list arriving at 100. Next we divide our number (102) by 100.

102/100   = 1 
remainder = 02

The answer is 1 with a remainder of 02. We repeat the process on the carry, with the next smallest power of 10, in this case 10.

02/10     = 0 
remainder = 2

We repeat the process, dividing by 1, to get a carry of 0 and we're done

2/1 i     = 2 
remainder = 0

From this we find

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

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. We see that 0 and 1 are the same in the binary and the decimal system. (You may wonder if this is obvious or trivial or why I'm even saying it. Some of the numbers will be different in the decimal and the hexadecimal system, so wait till we get to the hexadecimal system.)

binary decimal
   0  =  0
   1  =  1

Here's a binary number. Since the base is 2, here's what it represents.

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

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

What is the decimal representation of 1101₂?

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

leading zeroes behave the same as for decimal numbers so

00=0
01=1
01101=1101

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, converting 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]

You ocassionally need to convert between decimal and binary. Many languages have built-in facilities for doing this.

Here's python code to convert decimal to binary (http://www.daniweb.com/code/snippet285.html). You'll need to know more about programming to use this.

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

Python can do the conversion directly for you (example at convert decimal to binary in python http://stackoverflow.com/questions/3528146/convert-decimal-to-binary-in-python). Fire up the python interpreter

python >>>

Then enter the following string followed by a carriage return. This converts 16₁₀ to binary.

python >>> bin(16)
'0b10000'

The prefix "0b" is python's way of saying the following number is binary. Convert from binary to decimal

>>> 0b1001
9

There's often several ways of doing the same thing in a language. Here's another in python from easy binary-decimal-denary conversions in python (http://movingtofreedom.org/2010/10/28/easy-binary-decimal-denary-conversions-in-python-2-6/)

>>> int ('1001', 2)
9
>>> int ('1001', 10)
1001

The second number is the base. The output is decimal.

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 in groups of 2ⁿ (e.g. 2,4,8,16,32,64,128) bits. For manufacturing, if you wanted to make something bigger, you make another copy of the piece of hardware next to the old piece, leading to increases in hardware in steps of 2. 8 bits is called a byte(B). Initially computer displays were text only and to display numbers, punctuation and letters (both lc and UC) you needed to be able to describe about 80 characters. An 8 bit (1 Byte) counter can handle numbers from 0..255 and can hold enough (80) numbers to represent standard text. A 4 bit counter (0..15) is too small to hold text. Because of the requirement for computers to manipulate characters, the byte became the standard size grouping of bits to manipulate at one time.

some 1Byte numbers expressed in decimal

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

The convention for the integer value of a byte that represents a character is called ASCII (http://en.wikipedia.org/wiki/ASCII) (the American Standard Code for Information Interchange). Here are some byte values for ASCII

byte		hexadecimal	decimal ASCII
00110000	0x30		48	'0'
00110001	0x31		49	'1'
.
.
00111001	0x39		57	'9'
.
.
01000001	0x41		65	'A'
01000010	0x42		66	'B'
.
.
01011010	0x5A		90	'Z'

01100001	0x61		97	'a'
.
.
01111010	0x7A		122	'z'		

Single quotes are used to denote a character e.g. 'A'. If you see an 'A' on your screen, the value stored in memory will be 01000001₂ == 65₁₀.

There is a difference between a character and a number. A character is a glyph (arbitary shape) on a screen. The character '0' is a glyph used to represent the value 0. If the computer wants to store the value 0, it stores 00000000₂. If the computer wants to write the character '0' on the screen, it's stored as 00110000₂. The program accessing the memory knows whether the storage represents a number or a character (or something else) and interprets it appropriately. The character 'Z' is stored in the machine as 01011010₂. 'Z' is not the value 90₁₀.

Let's say the computer is storing the number 10100000₂=160₁₀=A0₁₆. If the user asks the computer to display the number in decimal format on the screen, the frame buffer will contain '1','6','0' i.e. 00110001,00110110,00110000. If the user asks the computer to display the number in hexadecimal format, i.e.'A','0', the frame buffer will contain 01000001,00110000.

The only difference between the binary representation of UC and LC characters is that the 3rd most significant bit (00x00000) is flipped. This simplifies changing between UC and LC.

Convert the following decimal numbers to binary

3
14
27
72

Convert the following binary numbers in decimal

11
101
1100
10110110

Bytes are always 8 bits. However data is shifted around the computer according to the bus width of the computer (the number of wires/connections between the CPU and the rest of the computer). Most PCs (in 2008) were 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. 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. In 2012, many PCs are now 64bit machines. The term word comes from the HPC (high performance computing or supercomputers) world, where 64 bit computing has been the standard for 30 yrs, 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. A 64-bit computer has a 64-bit (8 byte) word size.

On a 32-bit machine, what is the largest integer that you can represent? ^[4] You can use python to calculate this for you.

>>> 2**32-1
4294967295L

A 32-bit machine can generate 4G different numbers (0..4G-1). If the computer uses these numbers as addresses, the computer than can address 4G different objects.

The way computers are built, memory is addressable as bytes. The machine reads or writes the whole byte at a time.

	Note
	While the reads and writes to memory occur in bytes, once the byte is fetched into the CPU, the bits can be changed individually.

What is the largest amount of memory that a 32-bit computer can address? ^[5]

How much memory is on your computer? In Windows do control panel->system.

A computer reads and writes to the hard disk. The computer accesses the hard disk in blocks. On PCs the block size is fixed for a particular machine and is 512-4096 bytes. (For Windows, the default block sizes are at Default Cluster Size http://support.microsoft.com/kb/314878. In Windows you can find the block size for your computer by running chkdsk C:. This will take a few minutes, and at the end the computer will display the block/sector size.) When a computer reads or writes to the hard disk, it reads or writes a block as one unit. A computer does not read or write parts of a block.

What is the maximum amount of storage a hard disk can have on a 32-bit computer if the block size is 512 bytes? ^[6]

^[7]

How much storage do you have on your hard disk? To find this in Windows go to explorer, then the properties of your disk.

The information displayed on the computer display is stored in separate memory on the video card. The display is made up of pixels (pixel = picture element i.e. the dots on the screen). Early PCs did not have a lot of video memory and the displays were monochrome. The eye can differentiate about 200 different intensities of light. In grayscale (http://en.wikipedia.org/wiki/Grayscale), on the right hand side at the top, look at the 16 level grayscale. How many bits are needed to represent a 16 level grayscale? ^[8] If the shading is more than 200 different intensities, the grayscale will appear to be continuous to the human eye. How many bits will you need to make a continuous grayscale? ^[9]

In the following, give the answer in bytes, using appropriate prefixes.

The early computer displays were monochrome and text only. Being text only, the computer hardware could only render letters of a certain type and size. Eventually color displays capable of displaying graphics became available. In a graphics display, every pixel on the screen can be independantly adressed. On a graphics screen you can draw arbitary shapes, e.g. a circle, while on a text only display, you can only display letters that are built into the hardware.

On a computer, color is composed of 3 channels; red, green, blue (r,g,b) allowing the display to present most of the colors in the rainbow and their combination. If you look at a white pixel on a computer display with a 5X magnifying glass, you will see the (r,g,b) subpixels. What is their physical arrangement (i.e. what do you see) ^[10] ? How many pixels are on your display (in Windows, Control Panel->Display. You'll get the w*h in pixels.)?

Digital and analogue cameras and TV sets all use the same (r,g,b) method of generating colors. This is because the human eye has 3 color receptors, which are also (r,g,b). The color receptors are in the retina and they're called cones. The sensitivity of the three types of cones, as a function of wavelength, is shown in the graph at Mechanism_of_trichromic_color_vision (http://en.wikipedia.org/wiki/Trichromacy#Mechanism_of_trichromic_color_vision). Note that the red and green cones overlap in wavelength, while the blue has little overlap. Light of 500nm will activate all three cones (S,M,L). The unit on the x-axis is what ^[11] What is the shortest and longest wavelength that humans can see ^[12] ? You can call these cones the r,g and b cones, if you like, and many people do, but people who study vision call them the S,M and L cones (short, medium and long wavelength). Red has the longest wavelength, while blue is the shortest.

To get an idea of the color at each wavelength, look at the diagram labelled "Relative brightness sensitivity of the human visual system as a function of wavelength" at Color_vision (http://en.wikipedia.org/wiki/Color_vision).

Here's an explanation of the rgb color model (http://en.wikipedia.org/wiki/RGB_color_model). To explore the combinations of colors, using the rgb model, look at color cube (http://www.morecrayons.com/palettes/webSmart/colorcube.php). (The color cube is an idea of Maxwell.)

In the rgb model the color is described by 3 numbers running from 0 (no color) to 1 (full color). Thus (r,g,b)=(1,0,0) is red (abbreviated 1,0,0). White is (1,1,1). (r,g,b)=(1,1,0) has full red and full green but no blue. What color is (1,1,0)? Try a few other combinations. The pairwise colors (1,1,0), (1,0,1), (0,1,1) have common names. What are they? ^[13] What is (1,1,1)?

If you wanted to present to the user all colors in continuous tones, how much memory would a computer need to store a pixel ^[14] ? It's convenient hardware-wize to manipulate bytes in groups of 2^n. In this case computer manufacturers use 4 bytes for each pixel, the 4th byte being used for transparency. A typical computer display is 1024 pixels wide * 768 pixels high. How much memory is needed to display full color on this screen ^[15] ? Modern video cards have lots more memory than this. It's used to store intermediate results when calculating the intensity of pixels.

If you want to look at how colors used to look on computers before full color displays became available, look at monochrome and rgb palettes (http://en.wikipedia.org/wiki/List_of_monochrome_and_RGB_palettes).

When you looked at your computer display, you might have noticed that the blue was darker than the r,g. This is because only about 5% of the sun's (white) light is blue. Thus the color we see as white is (r,g,b)=(45%,50%,5%). If there's any more than 5% b, the color looks blue, rather than white. When we recreate white, we don't need much b, just 5% of the total energy.

If you look at a color (e.g. blue, for long enough (30-60secs), the cones will become bleached (respond less to that color). If you then move your focus to a white area, you will see the complementary color (in this case you'll see white without the blue, e.g. only the red, and green = yellow). As an example follow the instructions on after image (http://www.yorku.ca/eye/afterima.htm).

The cones are relatively insensitive to light and only function in daylight. We have another light receptor, called a rod, which sees black/white/gray only. Rods are most sensitive at 498nm (green). Here is the light sensitivity of rods and cones (http://www.yorku.ca/eye/specsens.htm). The rods are more sensitive to light than are cones. At night only the rods are being use (this is why night scenes have little color). Despite what this curve shows, due to the much greater sensitivity of the rods, the rods detect light beyond the red and blue ends detected by the cones. Thus the color seen at the two ends of the rainbow is grey.

It's important to differentiate a colored light from a colored pigment. A red colored light emits red light (about 650nm). A green colored light emits green light (about 550nm). If you mix both lights together (with a prism, or by shining both on a piece of white paper), you will see both colors together, which the brain interprets as y. Thus lights are additive (r+g=y). For an image of adding lights see RBG illumination (http://en.wikipedia.org/wiki/File:RGB_illumination.jpg) or additive color (http://en.wikipedia.org/wiki/File:AdditiveColor.svg).

On the other hand, a red colored pigment (e.g. paint, printed dye on paper) reflects red. It reflects red, because when illuminated by white light, it absorbs the green and blue. Thus pigments are subtractive (http://en.wikipedia.org/wiki/Subtractive_color). If you mix red and green pigments, the red absorbs g,b, while the green absorbs r,b. All colors will be absorbed and none will be reflected. Thus the mixture of green and red pigments, when seen with white light, will be black (in practice they'll be a dirty dark color, because the pigments don't absorb all light).

What will happen if you mix c and y pigments ^[16] ?

Color printers use (c,m,y) dies, which are called the primary subtractive colors. You can't make a very good black by adding (c+m+y), so printers use have black (k) too (the symbol b being taken by blue). (As well b ink is cheaper than using up your c+m+y inks.) A printer then is (c,m,y,k).

The primary lights then are (r,g,b). The primary pigments are (c,m,y). When someone says "what are the primary colors?" you have to say "Do you mean primary lights or primary pigments?". (c,m,y are also called secondar lights, but you rarely hear the term used.)

The original experiment of splitting white light into its components with a prism and recombining them, was done by Newton. Newton found that the eye cannot distinquish spectrally pure yellow (produced from a prism) from a combination of red and green (each produced separately by a prism and then combined), even though physically these are distinct ways of creating the color yellow (a prism will show that yellow is still yellow, while the yellow which comes from a combination of red and green, will be split again by a prism into red and green).

To show trichromic vision, I'm now going to do a demonstration with colored light (white light filtered through r,g,b,c,m,y cellophane). If you illuminate a white object, with red light what color will the object appear to be ^[17] ? If you illuminate a red object with red light, what color will the object be ^[18] ? If you illuminate a green object with red light, what color will the object be ^[19] ? If you illuminate a blue object with red light, what color will the object be ^[20] ? If you illuminate a black object with red light, what color will the object be ^[21] ?

Now watch the demonstration.

Assume that white is equal amounts of r,g and with 5% b and that the eye isn't terribly sensitive to small changes in brightness. Group r,g,b,c,m,y according to their perceived brightness ^[22] . Note that c is also perceived as being bright like y (c is used on fire trucks in Chapel Hill, although more likely for political than practical reasons). This is because you can make a color (0.5,1,1) which looks the same to the eye as (0,1,1). Thus you can make a bright c.

The table below shows the perceived color of an object when illuminated by various colors. The r column is filled for you. Fill out the empty columns in this table.

The S and M pigments diverged about 500Mya (Trichromic vision in Primates http://physiologyonline.physiology.org/content/17/3/93.full). What geologic age was this, what sort of life forms existed on earth then ^[23] ?

The big event of the Cambrian, as far as life forms go, was the Cambrian Explosion (http://en.wikipedia.org/wiki/Cambrian_explosion), a period of rapid diversification of life forms and the transition animals from soft bodied to hard bodied. At about 540Mya, trilobites acquired segmented compound eyes (as found in modern arthropods, like insects). Previously trilobites has light sensitive patches, but over a period of about 20My, they developed eyes which could form images, and presumably allow them to locate prey.

Following the aquisition of eyes at 540Mya, it would appear then that by about 500Mya, separate S and M receptors has evolved.

While humans are trichromic, most mammals are dichromic. From Trichromic vision in Primates (http://physiologyonline.physiology.org/content/17/3/93.full).

Exceptions to dichromacy are rare.

Many marine mammals and a few nocturnal rodents, carnivores, and primates have secondarily lost the S cone pigment and become monochromats.

Many diurnal primates, on the other hand, have acquired a third cone pigment, the L cone pigment, which is maximally sensitive to the longer visible wavelengths (red).

The aquisition of the L (red) receptor allowed the primate to differentiate fruit (red, orange and yellow), which constitutes 90% of the diet and young leaves (slightly red), which are eaten as a last resort in the absence of fruit, from the bulk of the green older leaves. At the same time a significant loss of olefactory capacity occured. Presumably it was less risky to spot a ripe fruit with your eyes than to climb a tree and smell the fruit only to discover that it wasn't ripe yet.

Although vertebrates are primitively tetrachromic (many teleost fish, reptiles and birds have four different types of cone cells, each sensitive to distinct light spectra), the mammals are generally dichromic (having only two cone types in the eye). A popular notion is that this reflects a transition to a nocturnal way of life. In the primates, however, trichromic vision has re-evolved - there are three types of cone cell, differentiating colour in the red, green and blue regions of the spectrum. Trichromacy is linked not only to largely diurnal activity, but also to the adaptive advantages of being able to detect coloured fruit and leaves. An important corollary of this argument is that with the shift to acute vision the role of olfaction has become less important. There is probably a lot of truth in this idea, and certainly relative to other mammals (notably the rodents) human olfactory capacity is seriously impaired with many of the relevant genes turned off by being transformed into pseudogenes. (see trichromic vision in mammals http://www.mapoflife.org/topics/topic_328_Trichromic-vision-in-mammals/).

The dichromic vision of mammals may account for their less colorful coats (e.g. cow, moose) compared to birds. From the color of the skins of cartilaginous fish (e.g. sharks), would you expect they have tetrachromic vision as do the telost fish?

Then about 100Mya (lecture) primates acquired trichromic vision (http://mitworld.mit.edu/video/669). Also see Evolution of color vision in primates (http://en.wikipedia.org/wiki/Evolution_of_color_vision_in_primates).

Some human females have tetrachromic vision (http://en.wikipedia.org/wiki/Tetrachromic_vision). Also see tetrachromacy (http://en.wikipedia.org/wiki/Tetrachromacy).

Although beyond the scope of this class, it's worth noting that the brain does not directly receive the (r,g,b) information. The image is subject to neural processing in the retina. (see Trichromic vision in Primates http://physiologyonline.physiology.org/content/17/3/93.full).

The dynamic range over which the eye works is about 10⁷ (day to night). A black road in daylight is emitting more light than does a white sheet of paper indoors, but the brain perceives the road as black, and the paper as white. One mechanism for achieving high dynamic range is encoding intensity of light as the difference between the intensity at the receptor and those immediately surrounding it.

Color is encoded as the difference between b-y and r-g i.e. as a blue-yellow channel and as a red-green channel.

see http://www.yorku.ca/eye/trichrom.htm http://www.yorku.ca/eye/opponent.htm http://www.yorku.ca/eye/huecan.htm http://www.yorku.ca/eye/handj.htm

A computer display is too dim to use in bright light (broad daylight). Let's see what contributes to this.

The computer display is backlit i.e. a fluorescent light at the side(s) of the display illuminates a white sheet behind the display. The white light shines through the pixels, to deliver the light you see. Assume that white light is equal amounts of (r,g,b) (it's not; what are the approximate amounts of energy in each color r,g,b?). Also assume that each (r,g,b) subpixel is the same size.

If a pixel is (1,0,0), what color are we seeing ^[24] ? How much light is getting through to the viewer? Here's the pixel

----------------------
|      |      |      |
|      |      |      |
|      |      |      |
|      |      |      |
|  r   |      |      |
|      |      |      |
|      |      |      |
|      |      |      |
|      |      |      |
----------------------

Only r is on. What fraction of the white backlighting is being presented at the back of the pixel is illuminating the r subpixel ^[25] ? What fraction of the white light being presented at the back of the r subpixel, is being seen by the user ^[26] ? What fraction of the white light being presented by the backlight to the whole pixel is being seen by the viewer ^[27] ? If the pixel was (1,1,1) instead of (1,0,0), what fraction of the backlight would be seen by the viewer ^[28] ?

At most, only 1/3 of the backlighting gets through, even when the screen is fully white. If you were dichromic (or tetrachromic) with an appropriate TFT display, what fraction of the backlight would get through to you on a white screen?

The simple fix for having a computer display usable in daylight, is to cut out the metal case at the back of the display, allowing sunlight to shine through the TFT pixels. (You will have to be careful not to punch through various protective and functional layers behind the TFT pixels.)

Let's do some binary addition: The rules are the same as for 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. Note: 1+1!=10 (where "!=" is the symbol for "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₂ ^[29] ? What decimal numbers are represented?

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.the algorithm of 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? ^[30]

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, as a function of problem size, is called the Order of the algorithm.

Speed (Order) of the addition operation:

Addition of a pair of 1 byte numbers, one digit at a time (as we did above), takes 8 steps. The Order of addition, doing it 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 (note the first carry of 0 is put into the 2nd column from the right)
 ---
 001 add
1    carry
----
1001 add
0    no carry - we're done 

	Note
	The last carry/add is done in computer hardware as part of the adding the carry and add so there's really only two steps.

Addition done in parallel takes 2 steps, no matter how long the numbers being added. The Order of parallel addition algorithm 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 O(1) parallel addition algorithm scales better than the O(n) stepwise algorithm.

To do parallel addition, 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 the billions of addition operations that will be done during the lifetime of the computer.

What is amortization (http://en.wikipedia.org/wiki/Amortization)?

Amortization is the process of saving money, in a large project, by paying a large amount at the beginning.

An example: Let's say you want to dig a small hole in the ground. The hole can be dug by one person using a shovel in 1hr at a cost of $20. Now let's say you want to dig a hole 100 times bigger. You could hire 100 people with shovels at a cost of $2000. Or you could rent a BobCat ($500) to do the same job. You have amortised the larger cost of the BobCat over the bigger job, saving you money.

Let's say in the design of a computer, you have the choice of two pieces of hardware; a fast expensive piece and a slow cheap piece and the piece you choose only has to be paid for once, and you can use that process as many times as you want. In this case you choose the expensive piece of hardware, and you are said to be amortising the extra cost over the life of the machine.

e.g. cost of stepwise adder=$10. cost of parallel adder=$20. If you're going to be doing 10¹⁵ 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 hiring staff and buying machines 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. As well we don't care if the constant of proportionality is different for each process. 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+
 ----

^[31] ?

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. If both numbers have a 1 in the most significant bit, then you'll have 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

7.4.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 by specially designed hardware.
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.

7.5.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 computer, 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? ^[35]

Overflow is required to make this process work. 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).

In a 1 dit computer there is no such thing as 10 (you need two digits to make a 10). So the process of subtracting the number from 10 doesn't work. Your computer doesn't have subtraction either. We have to do something else. Instead if you want the ten's complement of 7, you start with your list of numbers in inverse order (9,8..0) and you ask the computer to come in 3 places from the biggest number (here 9), i.e.8,7,6 giving 6. 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 (automatically 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 9+ 7 = 6 the answer the computer gave us 6

7.5.2. Two's complement

We've done subtraction in decimal using the ten's complement. Now we're going to do subtraction in binary using the 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 ^[36] ?

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 ^[37] ?
• 1011 ^[38] ?
• 1001 ^[39] ?
• 1100 ^[40] ?

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 ^[41]

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- --------

  ^[42] ?

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

  ^[43] ?

	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, for the rightmost digit, there is never a carry (there's no column to its right to produce a carry bit). However since the computer has a carry bit for the other 7 digits, it's simplest to have a carry bit for the right most column, and it's set to 0 at the start of addition. This carry bit is useful for the ADC operation (see below).

RH
10
11+
----
?? sum
?0 carry for right colum set to 0

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

step 2: 2nd least significant digit, add two digits + carry digit. There is overflow
RH
10
11+
--
01 sum
00 carry (with overflow)

The computer did 10+11₂ and got 01₂ (i.e. 2+3=1, with an overflow of 1 in the 3rd least significant digit, 100₂=4).

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 ^[56] ? They are the (-ve of the number + the base) (in a 1 byte computer the base is 256).

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

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 ^[58] . 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 ^[59] ? 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 ^[60] ?

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 ^[61] ?

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 ^[62] ?

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

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 ^[64] ?

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 $ ^[65] ?

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 ^[66] ? 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? ^[67] ) 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 ^[68] ?

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 ^[69] ? How many prime numbers could we address ^[70] ?

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 ^[71]

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) ^[72] ? How many letters down the alphabet is the character represented by hex '51' ^[73] ? Knowing that the hex for 'A' is 41h, figure out the hex for 'Z' and then try it ^[74]

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 ^[75] ?

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 ^[79] ?

	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

(see http://stackoverflow.com/questions/19365508/when-is-semicolon-use-in-python-considered-good-or-acceptable)

Many lanuguages use ';' as a statement separator and ignore white space. In these languages both of these pieces of code are identical

foo; bar;

foo;
bar;

Python's statements are separated by tabs at the start of the line, and not an end of statement character. In python the ';' separates statements on the same line. (However the practice of having multiple statements on one line is discouraged.) Thus "foo; bar" is OK, but "foo bar" is not OK. "foo; bar;" and "foo;" is OK - the null statement at the end of the line is an OK statement.

#OK
	foo 
	bar

#OK
	foo; bar


#error
	foo bar

#OK
	foo; bar;

#OK
	foo;
	bar;

Thus programmers coming from languages with statements terminated by ';' can continue to use the ';'. If you're going to write a lot of python, you should consider using the python convention, i.e. without the ';'

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 ^[85] .

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 ^[86] .

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 ^[87] .

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. ^[88] .

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 ^[89] .

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 ^[90] . 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 ^[91] ?

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 ^[92] .

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 ^[93] . 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 ^[94] .

	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 ^[95] .

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

loop_over (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.

Note

There are 3 parameters for range(), start, stop and step. Every language has an equivalent function for iterating through loops (everyone has to choose a different one). The default for step is 1 in all languages, so you can safely leave it out. You may remember the default for start for any language (it will be 0, 1 or undefined), or you may not. I haven't coded in python for a couple of years. It's easier and safer to put in values, even if they are the default values.

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 ^[96] .

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 ^[97] .

	Note
	End Lesson 9

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

When writing loops, often it is hard to figure out the number of times you have to run a loop. The problem is called the fencepost problem or the fencepost error.

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 ^[98] and how many fenceposts do you need ^[99] ?

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

When you're programming, it's not always obvious whether a counter is a fence section or a fencepost. The fencepost error 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 fence sections. You test every time.

If you've got it wrong and you only have one iteration to do, the computer will think it only has 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.

Here's another iteration problem suggested by a friend, Steve. You're at the State Faire and you toss balls to knock over milk bottles, to win a teddy bear for your date. The usual setup is a row of milk bottles; if you miss the one you're aiming for, you'll possibly hit one next to it. But this stall is run by a fiendish mathematician; he sets up one milk bottle at a time for you to hit. Let's say you're a good shot and have a 0.3 chance of hitting a milk bottle. You're given 10 milk bottles to hit. You need to hit 5 milk bottles to win the teddy bear for your date. What's your chance of winning? Is throwing the balls at the milk bottles a sensible thing to do?

This is a standard perms and coms probability coin toss problem; you have a coin which comes up heads with a probability of p=0.3. What's the probability of getting at least 5 heads in 10 tosses?

	Note
	For scientific programming, you're not programming the computer hardware. All the math will be known and you're just coding up a problem that you can work through on a piece of paper. After a while, the coding is relatively routine.

Since the exercise is about computing and not probability, here's the math involved for calculating probability.

Note

You can find the coefficients of the bionomial expansion with Pascal's Triangle (http://en.wikipedia.org/wiki/Pascal's_triangle). This works well for binomial expansions of small powers (say to the 5th power), but it's tedious for expansions of higher powers. What is the order of finding a single term in the bionomial expansion (say the probability of getting 3 heads and 7 tails)? For the first order (a single coin flip), you need to write town 1+2 terms, for two coin flips you need 1+2+3 terms. For 5 coin flips you need to write 1..5 terms. The number of terms increases as the square of the order or O(n^2). What is the order of finding terms using the equations with factorials? No matter what the order you only have to do one calculation. (Sure it's a complicated equation). The order is O(1), i.e. the cost is the same no matter how many coin flips you do. Pascal's Triangle scales much worse than the factorial method. It doesn't matter that for any particular number in the bionomial expansion, Pascal's Triangle is easier than the factorial method. It only matters that for a sufficiently large problem, the factorial method will always be faster.

First you need to know how to do factorials. To get factorials in python do this

#!/usr/bin/python
"""


"""

import math

#----

print math.factorial(3)

#--factorial.py----------------

	Note
	You should separate importing of libraries from your own code. For readability, I've separated out the line "import math" from the functional code with a comment line.

Next you need to format the output rather than using the free format built into python

	Note
	free formatting: you just want something out and you don't care about the format. You let python do it any way it wants. Free formatting is what you did with print "hello world"

Use this model program to print some output. It demonstrates output in decimal (integer), exponent, and floating point formats. The formatting token "%3.5f" indicates 3 figures before the decimal point and 5 figures after the decimal point.

#! /usr/bin/python

"""


"""

var_1 = 100
var_2 = 3.14159

print '%5d %3.5e %3.5f' %(var_1, var_2*100, var_2)

# formatted_print.py ---------

Now you're going to start writing bottles.py. Start with a section (delineated by comments) showing the import function, (you can copy factorial.py), then add another section for user defined variables (with the label "user defined variables" or "user input"), with num_bottles = 10 and probability = 0.3.

Note

For short programs that you'll only use a couple of times, you usually code input values into the program. These variables will be in a separate section (delineated by comments), so that the user knows this is the area he's allowed to mess with. Sometimes with editing, you can fat finger the code and loose a line or mess up code. You don't want the person editing to have their cursor sitting in any part of the code. If the program is going to be used by a lot of people, who don't know how to edit, you arrange for the user defined variables to be supplied on the command line.

The program will look like this so far.

#! /usr/bin/python

"""
        bottles.py
        Joe
        simulate Steve throwing balls at milk bottles
"""

import math

#------
#user entered input

num_bottles = 10
probability = 0.3

# bottles.py------

Note

I use the long variable name "probability" rather than say "p". In the early days of computing, you would have chosen "p". This is because you were restricted in the number of characters you could enter on a line (usually 80), entering your program was tedious (card punches or later, even worse editors than vi), keys were mechanical and requires strength. Now none of those factors apply. You have on-screen editors, and you'll spend a lot of time debugging. You're likely to want to look for all the instances where the value of the probability is used or changed. You'll never find them if you have to search for "p", but there's some hope if you use the string "probability". As well as the program grows larger, you're likely to want to change a variable name, and you'll want to do a bulk search and replace. You can't do that if you've used "p" for your variable. Loop counters only exist within the scope of the loop. Once the loop has exited, the loop counter doesn't exist any more. Usually loops are short (in terms of space on a screen), so while variables usually are given long explicit names, loop counter are often (but not always) given single letter names e.g. "i","j","k". (This is a convention from Fortran, where many of us learned to program. Fortran automatically declared "i","j" and "k" to be integers. As well "x", "y" and "z" were automatically declared reals.)

As you enter new lines of code, you will go from a program that works, to a program that doesn't work for a while. You don't want to get too far from working code. How far is too far is up to you. You might get a phone call or have to go somewhere at short notice, or realise that you're out of time for today. You don't want to come back to a mess of code that doesn't work. You'll never be able to figure out what the hell it was supposed to be doing. It's better that your code does something trivial and works, that doesn't do something that's amazing. No-one is impressed by things that don't work. You build your code a line or two at a time. Rarely do you ever have to write more than a couple of lines to get working code again.

You're going to need to keep a tally of the number of times a combination occurs. How many different combinations are there? This number may not be obvious. If you're unsure, flip a coin once. How many different combinations are there ^[102] ? You can have a head or a tail. How many possibilities are there if you toss a coin twice ^[103] ? You can have two heads, one head and tail or two tails. How many combinations are possible if you toss a coin n n times, or toss n balls at bottles ^[104] ? We're going to calculate the combinations in a loop. First time through we'll count the number of times it's all heads (or all hits), next time we'll count 1 tail and the rest heads (or only one miss).

Add a loop that just prints out the loop counter, for the number of different combinations of heads and tails (or bottle hits and bottle misses). Here's my code ^[105]

Next print out the number of permutations of that combination (there is only one way of getting all heads). Here's my output (yours doesn't have to look exactly like this, but the numbers have to be the same. The text strings are optional.)

dennis:class_code$ ./bottles.py
dennis:/src/www.austintek.com/python_class/class_code# ./bottles.py
number of hits 0 number of permutations 1
number of hits 1 number of permutations 10
number of hits 2 number of permutations 45
number of hits 3 number of permutations 120
number of hits 4 number of permutations 210
number of hits 5 number of permutations 252
number of hits 6 number of permutations 210
number of hits 7 number of permutations 120
number of hits 8 number of permutations 45
number of hits 9 number of permutations 10
number of hits 10 number of permutations 1

This is an unreadable mess. Put the strings at the top of the output like this.

dennis: class_code$ ./bottles.py
number        number       
bottles       permutations 
hit                        
    0               1
    1              10
    2              45
    3             120
    4             210
    5             252
    6             210
    7             120
    8              45
    9              10
   10               1

Nex print out the probability of each permutation as well. Here's my output. (you can stop python printing a carriage return at the end of a line, by appending ',' to the end of the print statement).

dennis:/class_code$ ./bottles.py
number        number          probability 
bottles       permutations    each        
hit                           permutation 
    0               1         2.82475e-02 
    1              10         1.21061e-02 
    2              45         5.18832e-03 
    3             120         2.22357e-03 
    4             210         9.52957e-04 
    5             252         4.08410e-04 
    6             210         1.75033e-04 
    7             120         7.50141e-05 
    8              45         3.21489e-05 
    9              10         1.37781e-05 
   10               1         5.90490e-06 

The probability of a certain number of hits is the product of the number of possible permutations and the probability of those permutations. Print out the probability too. What is the most likely number of hits?

dennis:class_code# ./bottles.py
number        number          probability     probability  
bottles       permutations    each            of number of 
hit                           permutation     bottles hit
    0               1         2.82475e-02     0.02825
    1              10         1.21061e-02     0.12106
    2              45         5.18832e-03     0.23347
    3             120         2.22357e-03     0.26683
    4             210         9.52957e-04     0.20012
    5             252         4.08410e-04     0.10292
    6             210         1.75033e-04     0.03676
    7             120         7.50141e-05     0.00900
    8              45         3.21489e-05     0.00145
    9              10         1.37781e-05     0.00014
   10               1         5.90490e-06     0.00001

It looks like we're close to home, except we don't have any evidence that we have a sensible answer. The output says that the most likely number of bottles hit is 3; this is reassuring. How about you add up the probabilities line by line as you go and output them in another column. The output for the last line of the iteration should be 1.0, with each line giving the probability of n or less hits.

Declare another variable cummulative_probability and output that as well.

Note

Initialising variables: Before you use a variable, you must give it a value. If you accumulate a value like this a = a + b or a += b then the variable a must have a value before it can be used on the right hand side of an equation. If you don't, depending on the language, you can get nonsense or the program will crash. Many languages, python included, automatically give values to uninitialised variables. In the top equation above, the variable a would likely be assigned int(0) and your program would give the expected result. However once you've learned one programming language, subsequent languages are easier to pick up, and it won't take long before you're using a language that behaves badly with unitialised variables. For cummulative_probability make a new section in your code, under the import function, called something like initialise variables and add the line cummulative_probability = 0.0

dennis:class_code$ ./bottles.py
number        number          probability     probability  cummulative
bottles       permutations    each            of number of probability
hit                           permutation     bottles hit
    0               1         2.82475e-02     0.02825      0.02825
    1              10         1.21061e-02     0.12106      0.14931
    2              45         5.18832e-03     0.23347      0.38278
    3             120         2.22357e-03     0.26683      0.64961
    4             210         9.52957e-04     0.20012      0.84973
    5             252         4.08410e-04     0.10292      0.95265
    6             210         1.75033e-04     0.03676      0.98941
    7             120         7.50141e-05     0.00900      0.99841
    8              45         3.21489e-05     0.00145      0.99986
    9              10         1.37781e-05     0.00014      0.99999
   10               1         5.90490e-06     0.00001      1.00000

Let's see what cummulative_probability is telling us. It's the sum of the probabilities of all the events listed so far. It tells us something like "at least" or "at most". The last line has a value of 1.0. It's telling us that the sum of the probabilities of 0->10 hits is 1.0. We want to know the probability of 5 or more hits. The sum of the probabilities from 0->4 hits is 0.85. This means we have a 15% chance of getting 5 or more hits.

Note

That "the sum of probabilities is the probability that any of those events can occur" was not obvious to me. As I found out, the statement in this last sentence is so imprecise, that I couldn't figure out what it meant or think clearly about it. As a first attempt to understand what was going on I looked at formal explanations, like Ask Dr Math: Ten Coin Flips, Four Heads and Probability of at Least 45 Heads in 100 Tosses of Fair Coin . The formulae from the Dr Math pages all looked simple enough. However being a scientist, I know that mathematicians can come up with formulae that will do anything at all, even prove that 1==2. Instead I go from the experimental results and then look for formulae. Here I write out all possibilities for a 4 coin toss and calculate by hand p_exactly(0) and p_at_least(0) and show that I get the same result as the Dr Math webpage. Here's the worked example; I toss a coin 4 times. What are the probabilities of exactly a certain number of heads, and the probabilities of at least a certain number of heads. The tosses and partial sums are grouped into blocks of four lines for readability. The columns are the configurations of the coins (i.e.H or T), the number of heads in the block, then two blocks, the exact number of tosses in the block with exactly that number of heads followed by the number of tosses in the block with at least that number of heads. coins H number_heads exact at_least 0 1 2 3 4 0 1 2 3 4 HHHH 4 HHHT 3 HHTH 3 HHTT 2 0 0 1 2 1 4 4 4 3 1 HTHH 3 HTHT 2 HTTH 2 HTTT 1 0 1 2 1 0 4 4 3 1 0 THHH 3 THHT 2 THTH 2 THTT 1 0 1 2 1 0 4 4 3 1 0 TTHH 2 TTHT 1 TTTH 1 TTTT 0 1 2 1 0 0 4 3 1 0 0 Total 1 4 6 4 1 16 15 11 5 1 p_at_least(0) = p_exact(0) + p_exact(1) + p_exact(2) + p_exact(3) + p_exact(4) 16 = 1 4 6 4 1 p_at_least(1) = p_exact(1) + p_exact(2) + p_exact(3) + p_exact(4) 15 = 4 6 4 1 p_at_least(2) = p_exact(2) + p_exact(3) + p_exact(4) 11 = 6 4 1 p_at_least(3) = p_exact(3) + p_exact(4) 5 = 4 1 p_at_least(4) = p_exact(4) 1 = 1 Explanation: there are 4 possibilities for exactly 1 H. The number of possible ways of getting at least 3 H (5) is the same as the number of ways of getting 3 H (4) and 4 H (1). Note: the binomial factors 1,4,6,4,1 for the totals in the exact column. If you add probabilities, then you are OR'ing the events. What this means Since you add the probabilities, I saw that to get "at_least(3)" in a 4 coin toss, you were OR'ing "exactly(3)" OR "exactly(4)". Thus the event where you get 3 heads had to be independant of the event when you get 4 heads. I thought I was doing one 4 coin toss and calculating all the probabilities. What I was really doing was observing a 4 coin toss where I got 3 heads, then calculating the probability that this occured, then observing a 4 coin toss where I got 4 heads and then calculating that probability. Since these two 4 coin tosses were independant events, the probability than one OR the other happened, was the sum of their probabilities.

Of course you quickly programmed up this problem before you decided to step forward to put down your money to throw the balls. Let's say the side show operator is selling the teddy bear for $25 and the ten balls cost $10. Is stepping forward an economically sensible thing to do? Let's say you have 100% probability of knocking over 5 bottles. Your cost for the teddy bear by throwing balls is $10. Clearly you chose to throw the balls. If you have a 50% probability of knockng over 5 bottles, the cost of the teddy bear is $20. Still you throw the balls. The breakeven point is 40% probability of knocking over 5 bottles, when the cost of the teddy bear by direct purchase and by throwing balls are both $25. What's the cost for you with a probability of 15%?

$10/0.15 = $66.67

If you want your date to have a teddy bear, the sensible thing is to pay the side show operator $25 for the teddy bear, and not throw the balls. That is if you're a total nerd at least.

It turns out that if you'd won the teddy bear by tossing balls, you would have got $10 worth of being cool in front of your date. (There's no cool from programming in python to figure out whether you should buy directly or toss balls.) So your real cost for the teddy bear is $35.

The end parameter for range is num_bottles + 1. This loop looks generally useful. It's possible (likely?) that sometime in the future, you'll want to mouse-swipe this piece of code and use it somewhere else. It would be better if there were minimal relics specific to the original program in the loop. Next time you aren't likely to be calculating anything about bottles. To fix your code make the end parameter num_combinations and calculate the value of num_combinations at the bottom of the block holding the user defined variables.

What you want is readabilility; you are looping over (number of combinations), not (number of bottles + 1). Sure here they are the same, but if you want to mouse swipe and reuse the loop elsewhere, the parameter for range will always be (number of combinations); it won't always be (number of bottles + 1).

Here's my code ^[106]

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)" ^[107]

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 ^[110] .

	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 ^[111] ? 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) ^[114] .

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 ^[115] ?

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 [116] ? How would you check that the pyc file was being used (in the absence of a py rather than the py file [117] ?

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 ^[118] ?

	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. ^[123]

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

	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 ^[125] 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 ^[126] ?

Note

You all know what 232 is [127] . 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 ^[129] and here's my output ^[130] . 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 ^[131] ? Why did the first pair pass, while the second pair failed ^[132] ?

The parameters passed to compare_results() are reals. How would you change the function to handle parameters than were string representations of numbers ^[133] ? (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 ^[136] .

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 ^[137] .

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 ^[138]

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 ^[141] and here's the output ^[142] .

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 ^[143] .

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) ^[145] ? (How many reals can you represent using 32-bits ^[146] ?) 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 ^[147] ? Use bc to convert this number to decimal ^[148]

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 ^[149] ?

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 ^[150] - trice ^[151] - four times ^[152] - five times ^[153] .

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 ^[154] ?

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 ^[155] ?) 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 ^[157] ? How many numbers are there between 0.25..0.5 and what is interval ^[158] ?

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 ^[159] ? How many floating point numbers exist in the range covered by 32-bit reals ^[160] ?

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 ^[161] ? What fraction of the numbers is sacrificed in each case ^[162] ?

	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 ^[163] ?

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 ^[164]

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 ^[165]

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

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 ^[167]

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 ^[168]

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 ^[169] ?

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 ^[170] ?) 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 ^[171]. 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 ^[172] ? 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 ^[173] ? 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 ^[174] ? Fix the code (you need to change the conditional statement). Here's my fixed version ^[175] .

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 ^[176] ?

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.111).

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 ^[177] ?) 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¹² is 2^{? [178]} 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 ^[179] .

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 ^[180] ?

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 ^[181]

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 ^[182] .

	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 ^[183] .

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 in a 64 bit double precision real number, i.e. every iteration 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 ^[184] ), 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 ^[186] . 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 ^[187] .

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() ^[188] .

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 ^[189] . 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. ^[190] . 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 [191] . 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 ^[192] .

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 ^[193] .

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. ^[194] ?

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() ^[195] ? What is the "1.0" doing in "1.0*x*x" ^[196] ?

	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 ^[197] ?

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 ^[198] 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 ^[199] 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 ^[200] ?) 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 ^[201] 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 ^[202] ?

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 ^[203] ?

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() ^[204] .

	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 ^[209] 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 ^[210] ? 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 ^[211] 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 ^[212] .

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 ^[213] ). Here's my version of pi_4.py ^[214] .

	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 [216] ?

Here's my code ^[217] . 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) ^[218] ?

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 ^[219] .

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 ^[220] .

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

We're doing runs that take multiples of 1000sec. Another useful number: how long is 1000 secs ^[221] ?

For large enough number of intervals, the running time is proportional to the number of intervals. It seems that the setup time for the program is about 0.03 sec and the setup time for the loop stops proportionality at about 10 intervals.

The difference between the upper and lower bound of π is 1/num_intervals (the area of the left most slice). If we want to decrease the difference by a factor of 10, we have to increase num_intervals (and the number of iterations) by the same factor (10), increasing the running time by a factor of 10. If it takes an hour to calculate π to 9 places, it will take 10hrs to calculate 10 places. With about 10,000 hours in a year, we could push the value of π out to all of 13 places. We'd need 100yrs to get a double precision value (15 decimal places) for π.

Here's the limits we've found on our algorithm:

If someone wanted us to calculate π to 100 places, what would we do ^[222] .

Calculating π to 100 places by numerical integration would require 10¹⁰⁰ iterations. Doing 10⁹ calculations/sec, the result would take 10^(100-9)=10⁹¹ secs. To get an idea of how long this is, the Age of the Universe (http://en.wikipedia.org/wiki/Age_of_the_universe)=13Gyr (another useful number) =13*10⁹*365*24*3660=409*10¹⁵ secs. If we wanted π to 100 places by numerical integration, we would need to run through the creation of 10^(100-9-17)=10⁷⁴ universes before getting our answer.

As Beckmann points out, quoting the calculations of Schubert, we don't need π to even this accuracy. The radius of the (observable) universe is 13*10⁹ lightyears = 13*10⁹*300*10⁶ * 365*24*3600=10²⁶m=10³⁶Å. Assuming the radius of a hydrogen atom is 1Å, then knowing π to 10² places, would allow us to calculate the circumference of the (observable) universe to 10²⁶/10¹⁰⁰=10^-64 times the radius of a hydrogen atom.

Despite the lack of any practical value in knowing π to even 10² places, π is currently known to 10¹² places (by Kanada).

We calculated π by numerical integration, to illustrate a general method, which can be used when you don't have a special algorithm to calculate a value. Fortunately there are other algorithms for finding π. If it turned out there weren't any fast algorithms for π some supercomputer would have been assigned to calculate 32-,64- and 128-bit values for π long ago; and they'd now be stored in libraries in our software.

While there are fast algorithms to calculate π, there aren't general methods for calculating the area an arbitary curve, and numerical integration is the method of choice. It's slow, but often it's the only way.

It can get even worse than that: sometimes you have to optimize your surface. You numerically integrate under your surface, then you make a small change in the parameters hoping for a better surface, and then you do your numerical integration again to check for improvement. You may have to do your numerical integration thousands if not millions of times, before finding your optimal surface. This is what supercomputers do: lots of brute force calculations using high order (O(n² or worse) algorithms, which are programmed by armies of people skilled at optimisation. People think that supercomputers must be doing something really neat, for people to be spending all that money on them. While they are doing calculations for which people are prepared to spend a lot of money, the computing underneath is just brute force. Whether you think this is neat or not is another thing. Supercomputers are cheaper than blowing up billion dollar rockets, or making a million cars with a design fault. Sometimes there's no choice: to predict weather, you need a supercomputer, there's no other way anymore. Unfortunately much of supercomputing has been to design better nuclear weapons.

I said we need a faster algorithm (or computer). What sort of speed up are we looking for? Let's say we want a 64-bit (52-bit mantissa) value for π in 1 sec. If using numerical integration, we'd need to do 2⁵²iterations=4*10¹⁵=10^15.6 iterations. A 1GHz computer can go 10⁹ operations/sec (not iterations/sec, but close enough for the moment; we're out by a factor of 100 or so, but that's close enough for the moment). We'd need a speed up of 10^15.6-9.0=6.6. If we instead wanted π to 100 places (and not 15 places) in 1sec, we'd need a further speed up of 10⁸⁵. There's no likelihood of achieving that sort of speed-up in computer hardware anytime real soon. Below let's see what sort of speed up we can get from better algorithms.

	Note
	Again for convenience let's assume that there is one real operation/iteration (or per interval) in the loop calculating π (in fact there are about 10).

We now have a value for π. As with all answers, we aren't home yet; we first have to know how good our answer is (i.e. we need an estimate of the errors too). Usually finding/calculating the errors in an answer takes longer than calculating the answer.

Due to the limited precision (32,64 bit) of registers used in computer real operations, rounding errors can (will) occur in every real operation and will affect the result. In the worst case, every floating point operation will have a rounding error and the errors will add, giving us 1 error for every operation.

	Note
	(In real life, not all calculations will result in a rounding error, and many roundings will cancel. But without a detailed analysis to find what's really going on, we have to assume the worst case.)

These errors are seen in the rightmost (least significant) bits of the number/mantissa. For every doubling of the number of intervals (which doubles the number of iterations), an extra bit of error is added at the rightmost end of the upper and lower bounds.

Here's a diagram of the lower bounds calculation after 2⁸ iterations on a 32 bit machine (23 bit mantissa). The least significant 8 bits have errors due to rounding.

	Note
	Joe: FIXME: do you get the same answer if you have 2^8 intervals and do 2^8 iterations, and in the case where you do 2^8 iterations of a 2^32 interval problem?

lower bounds: llllllllllllllleeeeeeee
                |         |         |
               20        10         0

l - lower bounds value
e - bit in error due to rounding

We've found that increasing the number of intervals by a factor of 10, gives us one more place where the upper and lower bounds match i.e. for every doubling of the number of intervals, an extra bit is added to our value for π. Thus as our calculation proceeds, for every doubling of the number of intervals, an extra bit on the left hand end of the register matches in the upper and lower bounds estimates, allowing us to add an extra bit to the value of π. Thus whenever an error bit is added to the right hand end of the register (the least significant end). an extra bit of precision is added on the left of the register (the most significant end) to our value of π.

Here's a diagram showing the bits in the 23 bit mantissa (32 bit machine) for the upper and lower bounds after 2⁸ iterations. Note that 8 bits match on the left end, giving us the value for π to 8 bits, while 8 bits of rounding error have accumulated on the right end.

upper bounds: mmmmmmmmuuuuuuueeeeeeee
lower bounds: mmmmmmmmllllllleeeeeeee
                |         |         |
               20        10         0

m   - matching (correct value for pi)
l/u - lower/upper bounds value
e   - bit in error due to rounding

With a sufficient number of intervals, the error bits and the good bits will meet in the middle. When this occurs you've calculated π to the maximum accuracy for that mantissa size. The number of intervals required to do this is the integer that fits into half the mantissa width.

Here's the maximum possible number of places for π that can be calculated by numerical integration, for various machine bit widths.

The maximum accuracy of π that can be calculated on 32,64 or 128 bit computers by numerical integration, is 11,26,56 bits (≅3,8,19 decimal places), taking 9msec, 16mins and 10³²yr respectively.

For the 32 and 64 bit machines, we need to make sure that we don't run the calculations any longer than needed, otherwise correct bits will be overwritten by errors. (For the 128 bit machine, there won't be enough time, to get to the maximum number of bits.) If the number if intervals is increased past this point, then errors coming from the right end of the register, will overwrite the already calculated good values of π sitting at the left end of the register.

Our calculation of π ran for 10⁹ (≅ 2³⁰) iterations/multiplication/additions. This table shows the combined effect of the register bit width (size of reals) and the rounding errors on the result.

In the case of the 23 bit mantissa, all 23 bits are overwritten by 30 bits of rounding errors. In the case of the 52 bit mantissa, the 30 error bits overwrite 8 of the good bits of π that were being calculated in the left most bits and we are left with only 22/3=7 correct decimal places. In the 128 bit case, 30 error bits were produced, 30 good bits were produced and the 52 bits in the middle hold values for the upper and lower bounds.

The error is relative to the number, rather than an absolute number. For slice heights which are less than half of the radius of the quadrant, the error in the heights will be reduced to one half. In this case approximately half the errors will not be seen (they will affect digits to the right of the last one recorded). We have two factors which each can reduce the number of rounding errors by a factor of 2 (half of the errors underflow on addition, and half of the additions don't lead to a rounding error). It's possible then that the worst case errors are only 29 bits rather than 30 bits.

Still it's not particularly encouraging to wait for 10⁹ iterations, only to have to do a detailed analysis to determine whether 29 or 30 of the 52 bits of your mantissa that are invalid. If we did decide to calculate π to 15 decimal places (64 bit double precision, taking 100yrs), doing 10¹⁵ operations, what sort of error would we have in our answer ^[223] ? One of the reasons to go to 64 (or 128) bit registers, is to have a large enough representation of the number that there is room to store all the errors that accumulate, while getting enough correct digits to get a useful answer.

In most cases with numerical integration, we don't know the answer (otherwise we wouldn't be doing the calculation). In this case our answer is the known π so we can use the known value to check our result. Our answer agrees with the known value of π to 9 rather than 7 places. This tells us that we didn't have the worst case errors. In fact we had about 100 times less errors than the worst case. Why was this? Was it due to skillful programming on our part? No; it was just dumb luck - most of the errors cancelled. We can expect as many rounding errors to round down as round up, so we should expect most of them not to be seen. However to know the number of errors we did get, we would need a more complicated error analysis than we've done here.

When we present our value for π from the numerical integration, what accuracy do we give; 7 places or 9 places? We found that the upper and lower bounds agreed to 9 places, so you might think that we can give 9 places. However our error analysis shows that in the worst case only 7 of these places are correct. We have to give our answer to 7 places. (Remember if we were determing the value of π by numerical integration, we wouldn't know its real value.)

We've gone through this numerical integration algorithm in great depth. Back when computers were first being used for modelling and numerical integration, all algorithms were new and they were explored in great depth to analyse speed and errors. Nowadays numerical integration is well understood and no-one would go through their algorithm in the detail we have just done here. They just know that if it's brute force and many iterations, then they'll have to run the code on a 128 bit computer. If you're working in that environment, the 128 bit computers will already be there and you will likely not even think about the number of bits you need (or have a choice of the bit width). (Supercomputing has been using 128bit reals since the 1970s.)

This section was originally designed to illustrate numerical integration. Numerical integration isn't practical for π: it's too slow, and has too many errors. We need another (quicker) algorithm, where most of the 52 bits of the mantissa (in a 64 bit machine) are valid. Since we've got numerical integration covered, let's look at other algorithms for calculating π.

The ancient Egyptians in 2500BC had a value for π accurate to 0.04% (see History of π http://en.wikipedia.org/wiki/Pi#History). The first attempts to calculate a value (by Archimedes, 287-212 BC) used numerical integration. Archimedes drew two sets of polygons, inscribing and circumscribing a circle (to get a lower and upper bounds) (see Archimedes' Polygons http://en.wikipedia.org/wiki/Pi#Geometrical_period) and by increasing the number of sides of the polygon, was able to get a value of π accurate to 4 places. It wasn't till the invention of calculus almost 2000yrs later, that faster algorithms were found to calculate π

There are many mathematical identities involving π and many of them are fast enough to be useful for calculating π. Here's one discovered independantly by many people, but usually attributed to Leibnitz and Gregory. It helped start of the rush of the digit hunters (people who calculate some interesting constant to large numbers of places).

pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...

Swipe this code with your mouse, calling it leibnitz_gregory.py.

#! /usr/bin/python
#
#leibnitz_gregory.py
#calculates pi using the Leibnitz-Gregory series
#pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...
#
#Coded up from 
#An Imaginery Tale, The Story of sqrt(-1), Paul J. Nahin 1998, p173, Princeton University Press, ISBN 0-691-02795-1.

large_number=10000000
pi_over_4=0.0

for x in range(1, large_number):
	pi_over_4 += 1.0*((x%2)*2-1)/(2*x-1)
	if (x%100000 == 0):
		print "x %d pi %2.10f" %(x, pi_over_4*4)

# leibnitz_gregory.py-------------------------

In the line which updates pi_over_4, what are the numerator and denominator doing ^[224] ? Run this code to see how long it takes for each successive digit to stop changing. How many iterations do you need to get 6 places (3.14159) using Leibnitz-Gregory and by numerical integration ^[225] ?

Newton (http://en.wikipedia.org/wiki/Isaac_Newton), along with many rich people, escaped London during the plague (otherwise known as the Black Death http://en.wikipedia.org/wiki/Black_Death, which killed 30-50% of Europe and resulted in the oppression of minorities thought to be responsible for the disease) and worked for 2 years on calculus, at his family home at Woolsthorpe. There he derived a series for π which he used to calculate π to 16 decimal places. Newton later wrote "I am ashamed to tell you to how many figures I carried these computations, having no other business at the time". We'll skip Newton's algorithm for π as it isn't terribly fast.

Probably the most famous series to calculate π is by Machin - see Computing Pi (http://en.wikipedia.org/wiki/Computing_%CF%80), also see Machin's formula (http://milan.milanovic.org/math/english/pi/machin.html).

pi/4 = 4tan^-1(1/5)-tan^-1(1/239)

Machin's formula gives (approximately) 1 significant figure (a factor of 10) for each iteration and allowed Machin, in 1706, to calculate π by hand to 100 places (how many iterations did he need?). The derivation of Machin's formula requires an understanding of calculus, which we won't be going into here. For computation, Machin's formula can be expressed as the Spellbach/Machin series

pi/4 = 4[1/5 - 1/(3*5^3) + 1/(5*5^5) - 1/(7*5^7) + ...] - [1/239 - 1/(3*239^3) + 1/(5*239^5) - 1/(7*239^7) + ...]

250 years after Machin, one of the first electronic computers, ENIAC, used this series to calculate π to 2000 places.

	Note
	End Lesson 29

This series has some similarities to Gregory-Leibnitz. There are two grouping of terms (each within [...]). In the first grouping note the following (which you will need to know before you can code it up)

In each term in the series consider

The second set of terms is similar to the first with "239" replacing "5". The second set of terms converge faster than the first set of terms (involving "5"), so you don't need to calculate as many of the terms in 239 as you do for the 5's to reach a certain precision (do you see why?). However for simplicity of coding, it's easier to compute the same number of terms from each set (the terms in 239 will just quickly go to zero).

Unlike ENIAC, today's laptops have only 64-bit reals, allowing you to calculate π to only about 16 decimal places. (There is software to write reals and integers to arbitary precision; e.g. bc has it. Python already has built-in arbitary precision integer math.)

The Spellbach/Machin's series has similarities with the Leibnitz-Gregory series. Copy leibnitz-gregory.py to machin.py and modify the code to calculate π using the Spellbach/Machin series. Here's my code for machin.py ^[226] and here's my output for 20 iterations

dennis: class_code# ./machin.py
x  1 pi 3.18326359832636018865
x  2 pi 3.14059702932606032988
x  3 pi 3.14162102932503461972
x  4 pi 3.14159177218217733341
x  5 pi 3.14159268240439937259
x  6 pi 3.14159265261530862290
x  7 pi 3.14159265362355499818
x  8 pi 3.14159265358860251283
x  9 pi 3.14159265358983619265
x 10 pi 3.14159265358979222782
x 11 pi 3.14159265358979400418
x 12 pi 3.14159265358979400418
x 13 pi 3.14159265358979400418
x 14 pi 3.14159265358979400418
x 15 pi 3.14159265358979400418
x 16 pi 3.14159265358979400418
x 17 pi 3.14159265358979400418
x 18 pi 3.14159265358979400418
x 19 pi 3.14159265358979400418

How many iterations do you need, before running into the 64-bit precision barrier of your machine's reals ^[227] ? How many iterations does it take to get 6 significant figures? Compare this with the Leibnitz-Gregory series ^[228] . What is the worst case estimate for rounding errors for the 64-bit value of π as calculated by the Machin series? Look at the number of mathematical operations in each iteration of the loop, then multiply by the number of iterations. Convert this number to bits, and then to decimal places in the answer ^[229] .

Two different formulae by Ramanujan (derived about 200 yrs after Machin) gives 14 significant figures/iteration. Ramahujan's series are the basis for all current computer assaults on the calculation of π by the digit hunters. You can read about them Ramanujan's series Pi (http://en.wikipedia.org/wiki/Pi). Coding these up won't add to your coding skills any more than the examples you've already done, so we won't do them here. How many iterations would Machin have needed to calculate π to 100 places if he'd used Ramanujan's series, rather than his own ^[230] ? How many iterations of Ramanujan's formula would you have needed to calculate π to 6 places ^[231] ? Ramanujan also produced a formula which gives the value for any particular digit in the value of π (i.e. if you feed 13 to the formula, it will give you the 13th digit of π).

π is a number of great interest to mathematicians and its value is required for many calculations. π is a constant: its value, as a 32- and 64-bit number, is stored in all math libraries and is called by routines needing its value. Since the value of π is not calculated on the fly by the math libraries, we have to look elsewhere to see the speed of standard routines that calculate π.

bc can calculate π to an arbitrary number of places. Remember that it would take the age of 10⁷⁴ universes to calculate π to 100 places by numerical integration. Here's the bc command to calculate π to 100 places (a() is the bc abbrieviation for the trig function arctan()). Below the value from bc, I've entered the value from Machin's formula. and the value calculated by numerical integration (with lower and upper bounds), using 10⁹ iterations.

# echo "scale=100; 4*a(1) " | bc -l
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170676
3.14159265358979 Machin, 16 places, 11 iterations
3.14159265(2..5) lower and upper bounds from numerical integration, 10 places, 10^9 iterations

You might suspect that the quick result from bc is a table look up of a built-in value for π. To check, rerun the example, looking up the value of a(1.01). If the times are very different, then the fast one will be by lookup, and the slow one by calculation. (If both are fast, you may not be able to tell if one is 10 times faster than the other.) Assuming the bc run took 1sec, it's 10^100-9=91 times faster than numerical integration.

Note that numerical integration gets the correct result (at least to the precision we had time to calculate). In another section we will attempt to calculate e and get an answer swamped by errors.

Most numbers of interest to computer programmers, whether available from fast algorithms or slow algorithms, are already known and are stored in libraries. Unfortunately fast algorithms aren't available for most calculations.

Quickly converging series have been found for integration of most analytic formulae (anything you can write as an algebraic function, e.g. the area under a parabola, a circle) and you use these series when you can. Most real world problems aren't describable by algebraic formulae. Numerical integration belongs to a group of methods called numerical methods. Numerical methods are used when algebraic methods don't work. Numerical methods are slow, brute force methods.

e.g.

Numerical methods handle discontinuities in data without a blink, while algebraic methods fail.

Numerical integration and other similar brute force numerical methods are the backbone of supercomputing. If you want to model the economy, calculate the flow of a fluid over a surface (will this plane fly?), model the weather, the oceans or pollution in the air, you'll be using numerical methods (and brute force). If you want 16 bits of accuracy (and 16 bits of error) each time step of your simulation will take 2³² iterations - about 1000secs (16mins). The overall calculation could take days or weeks.

That said, don't forget how much was accomplished before computer enabled numerical methods.

The code to calculate the area in a quadrant of a circle, can be extended to calculate the area under any curve, as long as its surface can be described by a formula, or a list of (x,y) coordinates.

How would we use numerical integration to calculate the length of the circumference of a circle ^[232] ?