Rolling dice: getting "33" vs "34"

bobbym · 2010-11-24 08:02:35

Hi qmech;

This is not working for me.

http://www.dartmouth.edu/~chance/teachi … pter11.pdf

It looks like you are attempting to diagonalize that matrix but I am not sure. Usually I see the formula as P^-1 X M X P. I am running it with the eigenvectors and the eigenvalues and I am not getting your sum, if that is what you want.

If you are looking for this sum:

As you can see it is a big mess. I did it with a CAS.

To compute it using C++ you are going to have to break it down into programming statements. Is that what you want, to program the sum?

To begin to understand roundoff error due to smearing we start with a question. You have single precision variables. They can hold values like 123456 or .123456 or 34.3265, only six digits at a time. If I ask you to sum this column of numbers with your six digit machine:

128632.816245
276152.7615
123. 23434
87231.91823
90852.83645

How would you do it? Just explain your method. 8 out of 10 people come up with a nice solution but never spot the problem when it occurs in real life.

qmech · 2010-11-24 10:24:44

Take the rightmost 4 digits of each number, with all numbers having the same number of decimal digits (it looks like 6 here). Add them up. If there's overflow and the rightmost (least significant digit) isn't zero, round so you still only have 4 digits to worry about. Repeat until the entire set of numbers is added.

Here I'm only using 4 instead of (maybe) 5 digits, and I'm only shifting left 1 digit at a time when I could possibly shift more, but I think this conveys the idea.

Are you pointing out that the intermediate calculations may need more than 6 digits to store the result?

bobbym · 2010-11-24 10:49:48

Hi;

Nope, what I am saying is that the only way to prevent smearing is to add the list in 2 pieces. The integer parts and the fractional parts. You would add them as 2 lists and then print the 2 lists.

128632 .816245
276152 .7615
123 .23434
87231 .91823
90852 .83645
--------- ------------

This way you keep the most digits. Otherwise if you add 128632.816245 + 276152.7615 = 404784, remember you can only hold 6. You lost a bunch of significant digits with with only one addition. All the digits after the decimal point got smeared out. Smearing occurs anytime a large number and a much smaller one are added. Anytime a multiplication occurs 101 * .123456 = 12.4690, actual answer 12.469006. We lost 2 digits, they were smeared out.

That is why that recurrence is not robust. It loses 1 digit per multiplication. Pretty soon the accuracy is totally gone. We say it multiples the error in the forward direction.

There are other types of error, truncation, subtractive cancellation, rounding.

This is a property of floating point arithmetic. Did you know that computers cannot subtract correctly?
That they cannot divide correctly? Everyone thinks computers are infallible. The only thing a computer can do is make a lot of errors really quick.

How do you get around this? Well, a long time ago Capt. Goldstine, John Von Neumann and Wilkinson invented Numerical Analysis.
The greatest math field there is! Basically it is the design of accurate models, the control of roundoff, truncation, smearing and cancellation error. How to coax your machine into giving accurate answers in calculations.
Designing algorithms that are robust. The field is huge and full of amazing examples.

1) Ever heard of the quadratic formula? Taught to every kid in the world. It is total garbage.

2) Difficult to demonstrate but if you add up a big list of numbers in one direction and then add them up in the opposite direction you may not get the same answer.
a1 + a2 + a3 + a4 +...+ an ≠ an +...+ a4 + a3 + a2 + a1

3) In mathematics the number line is a smooth line. With computers there are holes in it. 1 / 3 . 1 / 10 and infinite amount of others does not exist to them.

4) Many math identitites are not true for a computer. (a + b)(a - b) ≠ (a^2 - b^2)

5) If I start at 10 and count down by 1 / 10 it will not generally hit, 9,8,7,6,5,4,3,2,1,0

6) Computing a recurrence in the forward direction might produce incorrect answers due to smearing. Running it backwards produces extremely accurate answers.

qmech · 2010-11-24 11:35:39

How do you run a recurrence in the backwards direction?

bobbym · 2010-11-24 14:29:52

Hi;

Take this difference equation and get for me a[10] using single precision arithmetic.

Print out your a[0] to a[10], all eleven numbers. Bear with me it is important and it will show you what I mean.

Let's see how close Borland's compiler gets.

qmech · 2010-11-25 09:20:19

For float (single precision variables):

a[ 0] = 0.117783
a[ 1] = 0.057736
a[ 2] = 0.038114
a[ 3] = 0.028421
a[ 4] = 0.022634
a[ 5] = 0.018929
a[ 6] = 0.015234
a[ 7] = 0.020985
a[ 8] = -0.042877
a[ 9] = 0.454128
a[10] = -3.533023

For double (double precision) variables:

a[ 0] = 0.117783
a[ 1] = 0.057736
a[ 2] = 0.038114
a[ 3] = 0.028419
a[ 4] = 0.022647
a[ 5] = 0.018821
a[ 6] = 0.016099
a[ 7] = 0.014064
a[ 8] = 0.012486
a[ 9] = 0.011225
a[10] = 0.010196

Try this link:
http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/Chapter11.pdf

bobbym · 2010-11-25 09:58:29

Hi qmech;

Very good!

The actual answer is a[10] = .01019617597387411983158761... to 26 places. As you can see the single precision answer is gibberish not even having the correct sign.

To demonstrate that double precision or any precision is not the answer get for me a[25].
Please print out your answers in double precision, also what you get for ln(9 / 8).

qmech · 2010-11-25 12:04:43

These are results from a 16-bit machine, double precision.

a[ 0] = 1.17783e-01 (0.117783)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 9] = 1.12255e-02 (0.0112255)
a[10] = 1.01962e-02 (0.0101962)
a[11] = 9.33967e-03 (0.00933967)
a[12] = 8.61595e-03 (0.00861595)
a[13] = 7.99547e-03 (0.00799547)
a[14] = 7.46481e-03 (0.00746481)
a[15] = 6.94820e-03 (0.0069482)
a[16] = 6.91438e-03 (0.00691438)
a[17] = 3.50851e-03 (0.00350851)
a[18] = 2.74874e-02 (0.0274874)
a[19] = -1.67268e-01 (-0.167268)
a[20] = 1.38814e+00 (1.38814)
a[21] = -1.10575e+01 (-11.0575)
a[22] = 8.85057e+01 (88.5057)
a[23] = -7.08002e+02 (-708.002)
a[24] = 5.66406e+03 (5664.06)
a[25] = -4.53124e+04 (-45312.4)
a[26] = 3.62499e+05 (362499)

bobbym · 2010-11-25 14:44:06

Hi;

Okay, notice that the numbers are oscillating plus and minus. The actual answer is:

a[25]=.004 291 231 914 255 080 933 09 7532 748... to 30 digits.

The double precision answer is - 45312.4, 10 million times too large and having the wrong sign. What went wrong. Nothing, no matter what the precision the round off error will eventually catch up to you. The problem is that the difference equation, is ill conditioned or pathogenic.

This can be added to the above 6 in post #28. 7) Computers cannot even input numbers correctly.

Remember a[0] = ln(9 / 8) you put it in, in double precision ( around 17 digits ). But that was not really ln(9/8), it had a small error. Each time time you iterated throught the recurrence you magnified that error by 8 until the error swamped out the correct digits, that is what smearing is.
That is why there are no correct digits in your answer. In order to get 1 correct digit you would have needed to work with 25 digits of precision.

Do you understand why the 8 as a multiplier is causing the problem. If you do not understand tell me now or otherwise the fix will have no meaning to you.

qmech · 2010-11-26 05:00:06

I'm interpreting your remarks to mean ln(9/8) is just some number that is irrational (so it can't be represented exactly in a finite number of bits), and the 8 coefficient is just some number greater than 1. With any coefficient greater than 1, continued iteration with lead that to grow and grow and grow.

bobbym · 2010-11-26 08:56:39

Hi qmech;

Yes , that is almost right. Not only transcendentals and irrational numbers but remember post #28? Good old fractions like 1 / 3 and 1 / 10 ... are not representable in a finite number of bits.

bobbym wrote:

Anytime a multiplication occurs 101 * .123456 = 12.4690, actual answer 12.469006. We lost 2 digits, they were smeared out.

Anytime you have a number n, unless you are dealing with integers you have a small error called little e. You do not really have n but instead n + e. Everytime you do a calculation the little e is either no trouble because it remains small or...

Here is what happened.

1st iteration) 8 (n + e) = 8n + 8e. For the sake of brevity I will just chart the e.
2nd iteration) 8 (n + 8 e) = 64e
3rd iteration) 8 (n + 64 e) = 512e
3th iteration) 8 (n + 512 e) = 4096e etc.

The small e is now 4000 times larger than it was. It started at e = .0000005 and now it is approximately .002. In other words it is affecting the second and third digits after the decimal point. You already have lost 4 - 5 significant digits!

That is the reason the forward direction is impossible. Since you lost 1 digit for every iteration you would have had to start with ln(9/8) to 31-32 digits and all subsequent calculations to 31-32 digits to have 5 or 6 correct at the end.

Do you follow me so far?

qmech · 2010-11-27 05:29:44

So far so good. Yes, the e*8^n will grow huge after many iterations and swamp out your signal. So how do you avoid this without resorting to doing things analytically?

bobbym · 2010-11-27 07:42:56

Hi qmech;

If multiplying by some constant is unstable when that constant is greater than 1 than dividing by it must be stable.

We can run the recurrence backwards, first you solve for a[n-1].

Now we have the 8 on the bottom as a divisor. The error is going to work for us! What does this difference equation say. Well it runs from the higher numbers downwards ( backwards direction ). We need to start from say a[50] then we would get a[49], a[48], a[47]... all the way down to a[25].

Only one problem, how do we know what a[50] is? We need a starting value as the initial condition. How do we get that when we could not even get a[25]?
What would you do?

qmech · 2010-11-27 10:27:40

I agree that if the coefficient is less than 1, the error will get smaller and smaller. So rewriting the recurrence as you have does that wonderfully.

However we've lost a[0] - the starting point. It's not obvious to me where you start. If you assume the existence of a fixed point (for n going to infinity), you can solve for the fixed point, but I don't see how to get back down to finite n.

So what do you do?

bobbym · 2010-11-27 11:11:53

Start at a[50] and pick any number you want! Say a[50]=1241 and work your way down to a[25]

Run it in double precision and you will get all the digits!

a[25]=.004 291 231 914 255 080 933 09 7532 748...

qmech · 2010-11-27 12:03:26

Amazing! I confirm this for 2 starting values and 2 iteration lengths (listed below).

Although I can believe that errors will grow smaller and smaller when the coefficient is less than 1, it isn't obvious to me why the sequence values are asserting themselves so strongly. But clearly they do...

a[ 51] = 1.00000e+03 (1000)
a[50] = -1.24998e+02 (-124.998)
a[49] = 1.56272e+01 (15.6272)
a[48] = -1.95085e+00 (-1.95085)
a[47] = 2.46460e-01 (0.24646)
a[46] = -2.81479e-02 (-0.0281479)
a[45] = 6.23588e-03 (0.00623588)
a[44] = 1.99829e-03 (0.00199829)
a[43] = 2.59112e-03 (0.00259112)
a[42] = 2.58309e-03 (0.00258309)
a[41] = 2.65330e-03 (0.0026533)
a[40] = 2.71712e-03 (0.00271712)
a[39] = 2.78536e-03 (0.00278536)
a[38] = 2.85696e-03 (0.00285696)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)

a[ 51] = 5.55500e+03 (5555)
a[50] = -6.94373e+02 (-694.373)
a[49] = 8.67991e+01 (86.7991)
a[48] = -1.08473e+01 (-10.8473)
a[47] = 1.35852e+00 (1.35852)
a[46] = -1.67156e-01 (-0.167156)
a[45] = 2.36118e-02 (0.0236118)
a[44] = -1.73701e-04 (-0.000173701)
a[43] = 2.86262e-03 (0.00286262)
a[42] = 2.54915e-03 (0.00254915)
a[41] = 2.65755e-03 (0.00265755)
a[40] = 2.71659e-03 (0.00271659)
a[39] = 2.78543e-03 (0.00278543)
a[38] = 2.85695e-03 (0.00285695)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)

a[ 71] = 5.55500e+03 (5555)
a[70] = -6.94373e+02 (-694.373)
a[69] = 8.67984e+01 (86.7984)
a[68] = -1.08480e+01 (-10.848)
a[67] = 1.35784e+00 (1.35784)
a[66] = -1.67864e-01 (-0.167864)
a[65] = 2.28769e-02 (0.0228769)
a[64] = -9.36541e-04 (-0.000936541)
a[63] = 2.07019e-03 (0.00207019)
a[62] = 1.72535e-03 (0.00172535)
a[61] = 1.80046e-03 (0.00180046)
a[60] = 1.82412e-03 (0.00182412)
a[59] = 1.85532e-03 (0.00185532)
a[58] = 1.88673e-03 (0.00188673)
a[57] = 1.91933e-03 (0.00191933)
a[56] = 1.95307e-03 (0.00195307)
a[55] = 1.98801e-03 (0.00198801)
a[54] = 2.02423e-03 (0.00202423)
a[53] = 2.06179e-03 (0.00206179)
a[52] = 2.10077e-03 (0.00210077)
a[51] = 2.14125e-03 (0.00214125)
a[50] = 2.18332e-03 (0.00218332)
a[49] = 2.22708e-03 (0.00222708)
a[48] = 2.27263e-03 (0.00227263)
a[47] = 2.32009e-03 (0.00232009)
a[46] = 2.36956e-03 (0.00236956)
a[45] = 2.42120e-03 (0.0024212)
a[44] = 2.47513e-03 (0.00247513)
a[43] = 2.53152e-03 (0.00253152)
a[42] = 2.59054e-03 (0.00259054)
a[41] = 2.65237e-03 (0.00265237)
a[40] = 2.71723e-03 (0.00271723)
a[39] = 2.78535e-03 (0.00278535)
a[38] = 2.85696e-03 (0.00285696)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)

a[ 71] = 1.00000e+00 (1)
a[70] = -1.23239e-01 (-0.123239)
a[69] = 1.71906e-02 (0.0171906)
a[68] = -3.37236e-04 (-0.000337236)
a[67] = 1.88039e-03 (0.00188039)
a[66] = 1.63062e-03 (0.00163062)
a[65] = 1.69011e-03 (0.00169011)
a[64] = 1.71181e-03 (0.00171181)
a[63] = 1.73915e-03 (0.00173915)
a[62] = 1.76673e-03 (0.00176673)
a[61] = 1.79529e-03 (0.00179529)
a[60] = 1.82477e-03 (0.00182477)
a[59] = 1.85524e-03 (0.00185524)
a[58] = 1.88674e-03 (0.00188674)
a[57] = 1.91933e-03 (0.00191933)
a[56] = 1.95307e-03 (0.00195307)
a[55] = 1.98801e-03 (0.00198801)
a[54] = 2.02423e-03 (0.00202423)
a[53] = 2.06179e-03 (0.00206179)
a[52] = 2.10077e-03 (0.00210077)
a[51] = 2.14125e-03 (0.00214125)
a[50] = 2.18332e-03 (0.00218332)
a[49] = 2.22708e-03 (0.00222708)
a[48] = 2.27263e-03 (0.00227263)
a[47] = 2.32009e-03 (0.00232009)
a[46] = 2.36956e-03 (0.00236956)
a[45] = 2.42120e-03 (0.0024212)
a[44] = 2.47513e-03 (0.00247513)
a[43] = 2.53152e-03 (0.00253152)
a[42] = 2.59054e-03 (0.00259054)
a[41] = 2.65237e-03 (0.00265237)
a[40] = 2.71723e-03 (0.00271723)
a[39] = 2.78535e-03 (0.00278535)
a[38] = 2.85696e-03 (0.00285696)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)

bobbym · 2010-11-27 12:23:12

Each iteration picks up almost 1 correct digit.

Remember you start off like this. In this case e is large.
n + e

e/8, e/64, e/512, e/4096...

Pretty obvious that the error is chopped up. Small error = accurate answer!

Remember in the forward direction you needed 32 digits for ln(9/8) and for each iteration, to get 5 or 6 that are accurate. You lost 25 digits. Only natural to say backwards you pick up 25 digits of accuracy. You even can start with a negative number! As long as you iterate enough, any amount of precision is possible in this direction for this problem.

I think this idea is first attributed to Wilkinson, a really brilliant guy. It does seem like magic the first time you see it! It is part of the beauty of numerical analysis!

qmech · 2010-11-28 05:29:47

Thank you for showing that to me. Yes, it is beautiful.

In essence, the iterations generate the Taylor series expansion for ln(1+(1/8)). Is this correct? If so, this explains why you have the 1/n term. For exp(x), that term would have to be 1/n!, and you shouldn't need the minus sign....?

bobbym · 2010-11-28 08:10:28

Hi qmech;

In essence, the iterations generate the Taylor series expansion for ln(1+(1/8)). Is this correct?

That I cannot say. As for the 1 / n term, when I first saw this problem I had no inkling it could be part of a larger problem. Recently I found the excellent paper "An example of Round Off Error," by Bruno Guerrieri. Like him I first saw the problem in Cleve Moler's book. Only Bruno went much further with it. Apparently it comes from the attempt to integrate:

The recurrence arises naturally from the attempts to get an accurate answer! That is where the 1 / n comes from. I had to spend a lot of time in remembering and finding that particular recurrence. You see modern compilers and math packages have mechanisms to prevent errors due to smearing. To illustrate the point I had to find one that I was sure would defeat Borland's excellent compiler!

I assume you are using the Taylor expansion of ln(1 + x)?

qmech · 2010-11-29 06:19:42

If a[n-1] = (1/8) * ( (1/n) - a[n] ) and I start from a[3] and iterate to a[0], here's what I get:

a[0] = (1/8) * ( 1 - (1/8) * ( (1/2) - (1/8) * ( (1/3) - a[3] ))))))...)))))

The coefficient of the "1" is 1/8.
The coefficient of the "1/2" is -(1/8)^2
The coefficient of the "1/3" is (1/8)^3

The Taylor expansion of ln(1+x) = x - (1/2)*x^2 + (1/3)*x^3...

So when x=1/8, for a[0] these 2 series are the same.

bobbym · 2010-11-29 09:58:21

Hi qmech;

Have you tried plugging in 1 / 8 into what you are getting and into the Taylor series of ln(1 + x )?

qmech · 2010-11-30 07:47:21

I compared the values (1/2), (1/3)..., (1/20) and got (no surprises) the values of
log(1+(1/2)), log(1+(1/3))...log(1+(1/20)) for a[0].

It will be interesting to see if this works for other Taylor expansions (e.g., sin(x), exp(x), etc.).

bobbym · 2010-11-30 10:03:49

Hi qmech;

Good luck my friend. You are in the undiscovered territory at least as far as I know. If I ever can understand why that is working I will pass it on.

Math Is Fun Forum

#26 2010-11-24 08:02:35

Re: Rolling dice: getting "33" vs "34"

#27 2010-11-24 10:24:44

Re: Rolling dice: getting "33" vs "34"

#28 2010-11-24 10:49:48

Re: Rolling dice: getting "33" vs "34"

#29 2010-11-24 11:35:39

Re: Rolling dice: getting "33" vs "34"

#30 2010-11-24 14:29:52

Re: Rolling dice: getting "33" vs "34"

#31 2010-11-25 09:20:19

Re: Rolling dice: getting "33" vs "34"

#32 2010-11-25 09:58:29

Re: Rolling dice: getting "33" vs "34"

#33 2010-11-25 12:04:43

Re: Rolling dice: getting "33" vs "34"

#34 2010-11-25 14:44:06

Re: Rolling dice: getting "33" vs "34"

#35 2010-11-26 05:00:06

Re: Rolling dice: getting "33" vs "34"

#36 2010-11-26 08:56:39

Re: Rolling dice: getting "33" vs "34"

#37 2010-11-27 05:29:44

Re: Rolling dice: getting "33" vs "34"

#38 2010-11-27 07:42:56

Re: Rolling dice: getting "33" vs "34"

#39 2010-11-27 10:27:40

Re: Rolling dice: getting "33" vs "34"

#40 2010-11-27 11:11:53

Re: Rolling dice: getting "33" vs "34"

#41 2010-11-27 12:03:26

Re: Rolling dice: getting "33" vs "34"

#42 2010-11-27 12:23:12

Re: Rolling dice: getting "33" vs "34"

#43 2010-11-28 05:29:47

Re: Rolling dice: getting "33" vs "34"

#44 2010-11-28 08:10:28

Re: Rolling dice: getting "33" vs "34"

#45 2010-11-29 06:19:42

Re: Rolling dice: getting "33" vs "34"

#46 2010-11-29 09:58:21

Re: Rolling dice: getting "33" vs "34"

#47 2010-11-30 07:47:21

Re: Rolling dice: getting "33" vs "34"

#48 2010-11-30 10:03:49

Re: Rolling dice: getting "33" vs "34"

Board footer