How Can Throwing Rice Find π?

Use the square above. Press Start again whenever you like: your last runs stay on the chart as faint lines.

• a.Throw 10 grains at a time, five times. Write down the estimate after each throw. Does it get steadily better?
• b.Throw until you have 100,000 grains. How many decimals of 3.14159 does your estimate get right?
• c.Press Start again and do it once more. Do the two runs give the same estimate after 100 grains? After 100,000?

Why should rice know about circles? Because of area. The square runs from −1 to 1 across and from −1 to 1 up, so its side is 2 and its area is 2 × 2 = 4. The circle inside it has a radius of 1, so its area is π × 1 × 1 = π.

A grain that falls anywhere on the square, with every spot as likely as any other, lands inside the circle with a chance of π ÷ 4. So when many grains fall, about π ÷ 4 of them land inside:

grains inside ÷ all grains ≈ π ÷ 4, so π ≈ 4 × grains inside ÷ all grains

A computer “throws” a grain by picking two random numbers: x, how far across, and y, how far up, each anywhere from −1 to 1. Is the grain inside the circle? Its distance from the middle is √(x² + y²), by Pythagoras. It is inside when that distance is less than 1, which is the same as x × x + y × y < 1. No square root needed.

With blocks

The block program throws 100 grains on a square drawn on the stage. Grade 7’s random number block gives only whole numbers, so this square runs from −1,000 to 1,000 each way, and the circle’s radius is 1,000. (With whole numbers from −10 to 10, there would be only 21 places across for a grain to land, and the estimate would come out too low on average.) The variable d holds x × x + y × y, and the grain is inside when d < 1,000,000, which is 1,000 × 1,000.

In Python (from Grade 8)

Python’s random.uniform(-1, 1) gives a random decimal from −1 to 1, so the square can be the one in the figure. This program throws 400 grains and keeps a picture after every 25, as a film. pixels is this page’s own module for drawing squares of color (on a computer you would use a library such as Pillow).

import randomfrom pixels import show_frames random.seed(3)          # the same "random" grains every run (delete it for new ones)SIZE = 41               # the picture: 41 little squares each waypicture = []for row in range(SIZE):    line = []    for col in range(SIZE):        x = (col + 0.5) / SIZE * 2 - 1     # the middle of this little square        y = 1 - (row + 0.5) / SIZE * 2        if x * x + y * y < 1:            line.append(1)                 # inside the circle: lighter        else:            line.append(0)    picture.append(line) inside = 0frames = []for grain in range(1, 401):    x = random.uniform(-1, 1)    y = random.uniform(-1, 1)    col = min(int((x + 1) / 2 * SIZE), SIZE - 1)    row = min(int((1 - y) / 2 * SIZE), SIZE - 1)    if x * x + y * y < 1:        inside = inside + 1        picture[row][col] = 2              # a grain inside: yellow    else:        picture[row][col] = 3              # a grain outside: pale    if grain % 25 == 0:        frame = []        for line in picture:            frame.append(line[:])          # a copy of the picture now        frames.append(frame) print(inside, 'of 400 grains landed inside the circle')print('pi is about', 4 * inside / 400)show_frames(frames, ['#34403b', '#4d5c55', '#f4c542', '#ddd8c8'])

319 of 400 grains landed inside the circlepi is about 3.19

Picture

Press Play under the picture to watch the grains land. (line[:] makes a copy of a row; without the copies, every frame would show the last picture.)

Many grains

Now the throwing goes into a function, throw(grains), which gives back the estimate, so the program can try 100 grains, then 1,000, and so on up to a hundred thousand.

import random def throw(grains):    inside = 0    for i in range(grains):        x = random.uniform(-1, 1)        y = random.uniform(-1, 1)        if x * x + y * y < 1:            inside = inside + 1    return 4 * inside / grains random.seed(9)for grains in [100, 1000, 10000, 100000]:    print(grains, 'grains: pi is about', throw(grains))

100 grains: pi is about 3.21000 grains: pi is about 3.16810000 grains: pi is about 3.1468100000 grains: pi is about 3.14112

It finishes in about a second on a phone. Add 1000000 to the list for a million grains, and wait a few seconds: 3.142896. Change the seed to 10, 11, 12 and run it each time: the estimates change, but they always close in on 3.14159.

How far off, on average?

One run can be lucky or unlucky. To see what usually happens, this program makes 20 estimates with each number of grains and works out how far they were from the true π on average. math.pi is Python’s own value of π, 3.141592653589793, used here only to measure the miss. abs() turns a miss below π into a positive number, so misses on either side count the same.

import randomimport math def throw(grains):    inside = 0    for i in range(grains):        x = random.uniform(-1, 1)        y = random.uniform(-1, 1)        if x * x + y * y < 1:            inside = inside + 1    return 4 * inside / grains random.seed(6)for grains in [100, 1000, 10000]:    total_miss = 0    for run in range(20):        total_miss = total_miss + abs(throw(grains) - math.pi)    print(grains, 'grains: off by', round(total_miss / 20, 4), 'on average')

100 grains: off by 0.1408 on average1000 grains: off by 0.045 on average10000 grains: off by 0.0127 on average

Read the last program’s three lines together:

Ten times more grains made the miss only about three times smaller (0.14 to 0.045, then 0.045 to 0.013). A hundred times more grains made it about ten times smaller: 0.14 to 0.013. And ten times smaller is exactly one more right decimal place. So each new decimal of π costs about 100 times as many grains. (Why 100 and not 10 is a law of chance: the miss shrinks with the square root of the number of grains, and √100 = 10. Grade 11 explains it in More trials, a steadier answer.)

How far would that take you? This program follows the rule up from 10,000 grains, until the miss is small enough for five decimals, 3.14159:

grains = 10000miss = 0.0127            # the average miss with 10,000 grains, from the program abovewhile miss > 0.000005:   # five decimals right needs a miss under 0.000005    grains = grains * 100    miss = miss / 10    seconds = grains / 1000000    # a computer throwing a million grains a second    print(grains, 'grains: off by about', round(miss, 8))    print('   ', round(seconds), 'seconds, or', round(seconds / 60 / 60 / 24, 1), 'days')

1000000 grains: off by about 0.00127    1 seconds, or 0.0 days100000000 grains: off by about 0.000127    100 seconds, or 0.0 days10000000000 grains: off by about 1.27e-05    10000 seconds, or 0.1 days1000000000000 grains: off by about 1.27e-06    1000000 seconds, or 11.6 days

(1.27e-05 is Python’s short way of writing 0.0000127: move the point 5 places left.) A million million grains, and a computer throwing a million every second for eleven and a half days, to be all but sure of five decimals. (The program’s steps are 100 times each, so it overshoots: about 400,000 million grains, nearly five days, would already be within 0.000005 of π 19 times in 20.) Rice finds the first two decimals of π in a moment, and the fifth only with enormous patience.

That is the whole lesson of random sampling, and it is why scientists use it anyway. An estimate from random throws gets close quickly, and getting very close is very slow. For π there are far better ways: in China around 480 CE, Zu Chongzhi pinned π down to seven decimal places with polygons, and in Kerala around 1400, Madhava’s school found eleven decimals with sums of fractions (the Math, Going Further topic How Did Madhava Find π Without a Calculator? shows how). But for most real questions, how a cold spreads, how far a neutron travels in metal, what a game is worth, there is no formula at all, and throwing many random “grains” at the question is the only way to an answer. Grade 9 calls these Monte Carlo methods and tells how they got their name, in Algorithms that use chance.

The rule only works if every spot on the square is equally likely. Real hands don’t throw like that. Try throwing at the middle, the way a pinch of rice dropped from your fingers mostly lands near where you hold it.

Here is the same thing in Python. Each “aimed” number is the average of two even throws. Like the total of two dice, which comes out 7 more often than 2 or 12 (the Math, Going Further topic Which Total Comes Up Most with Two Dice?), an average of two crowds the middle.

import random def aimed():    # the average of two even throws: most land near the middle    return (random.uniform(-1, 1) + random.uniform(-1, 1)) / 2 random.seed(4)inside = 0for i in range(100000):    x = aimed()    y = aimed()    if x * x + y * y < 1:        inside = inside + 1print('Aimed at the middle: pi is about', 4 * inside / 100000)

Aimed at the middle: pi is about 3.89792

A hundred thousand grains, and the answer is off by 0.76. More grains make a random estimate steadier, but they can’t fix a biased one, where the throws favor some places. The same goes for a survey that only asks people in one town, or a test of a medicine given only to the healthiest patients. So a computer’s random numbers have to be good ones: spread evenly, with no pattern. They are worked out by a formula from a starting number, the seed, as Grade 9 explains (the same seed gives the same numbers, which is how the programs on this page give the same output every time).

Real rice. Draw a square of side 20 cm and the circle inside it. Drop 100 grains, one at a time, from about 50 cm up, without aiming, and count. How close to π do you get? Decide before you count whether a grain on the line is in or out. Do your grains spread evenly, or crowd under your hand?

A quarter is enough. Grade 9’s version uses only the top-right quarter: x and y from 0 to 1 (random.random()), and the quarter circle of radius 1. Change the throw program to do that. Does it still give π, and why?

Other shapes. Throwing grains measures the area of any shape, as long as a program can tell whether a point is inside it. Change the test to abs(x) + abs(y) < 1: what shape is that, and what area do the grains give it? Can you check the answer without a program?

★ Inside or outside? (blocks)

Finish the program so that it asks for a grain’s x and y (on the blocks’ square, from −1,000 to 1,000) and says inside if the grain is inside the circle of radius 1,000, and outside if it isn’t. A grain exactly on the line counts as outside. Then press Test.

★★ Rice in neat rows (Python, from Grade 8)

What if the grains weren’t thrown at random, but laid in neat rows, one on every point with whole-number x and y? Write grid_count(r) that counts the points (x, y), with x and y whole numbers from −r to r, that lie inside the circle of radius r or on it (x × x + y × y ≤ r × r). Then π is about grid_count(r) / (r * r).

def grid_count(r):    # count whole-number points (x, y) with x*x + y*y <= r*r    return 0 print(grid_count(5))

★★★ The luckiest throw (Python, from Grade 8)

With n grains, the estimate 4 × inside ÷ n can only take a few values, because inside is a whole number from 0 to n. If you were as lucky as possible, how few grains could give π right to a given number of decimal places? Write fewest_grains(decimals): the smallest n for which some number of grains inside, k, makes 4 × k ÷ n close enough, by the test in the starting program.

import math def fewest_grains(decimals):    # the smallest n for which some k from 0 to n gives    # abs(4 * k / n - math.pi) < 0.5 / 10 ** decimals    return 1 print(fewest_grains(2))

• Grade 7, Variables and Choices: random numbers in blocks, in a guessing game.
• Grade 8, Decisions and Loops in Python: counting inside a loop, as inside = inside + 1 does.
• Grade 9, Algorithms that use chance: Monte Carlo and its story, and random quicksort.
• Grade 11, Random simulations: Monte Carlo: how many runs are enough, and the 1 ÷ √(trials) law.
• How far does a wandering ant get? Another way chance grows like a square root: a random walk.
• How does a cold spread through a school? Chance, many times over, decides whether a cold dies out.
• Math, Going Further: How Did Madhava Find π Without a Calculator? π from a sum of fractions, eleven decimals by hand.