Posts Tagged ‘Data.Ratio’

Generalizing Heron’s method

November 25, 2012

Suppose we are given a positive rational number y. We would like to compute an approximation of the n-th root of y, which we denote by \sqrt[n]{y}. If n = 2, we can use Heron’s method to approximate \sqrt{y}, as we saw a few days ago. But, what if we have that n > 2? Is it possible to generalize Heron’s method?

Indeed, it is. Let us introduce a function f_{n,y} : \mathbb{R} \to \mathbb{R}, defined by f_{n,y} (x) = x^n - y, where n \geq 2 and y > 0. By construction, the positive zero of function f_{n,y} is \sqrt[n]{y}. Do note that if n is odd, then -\sqrt[n]{y} is not a zero of f_{n,y}. Recall that Newton’s method uses the recurrence relation x_{k+1} = g_{n,y} ( x_k ), where function g_{n,y} : \mathbb{R} \to \mathbb{R} is defined by

g_{n,y} (x) = x - \displaystyle\frac{f_{n,y} (x)}{f_{n,y}' (x)}

where f_{n,y}' is the first derivative of f_{n,y}. Hence, we obtain

g_{n,y} (x) = x - \displaystyle\frac{f_{n,y} (x)}{f_{n,y}' (x)} = x - \displaystyle\left(\frac{x^n - y}{n x^{n-1}}\right) = \displaystyle\frac{1}{n} \left( (n-1) x + \frac{y}{x^{n-1}}\right)

Thus, we have the generalized Heron’s method

x_{k+1} = \displaystyle\frac{1}{n} \left( (n-1) x_k + \frac{y}{x_k^{n-1}}\right)

If we make n = 2, we obtain the standard Heron’s method, as we expected. The following Haskell script implements this generalized Heron’s method using arbitrary-precision rational numbers (included in the Data.Ratio library) to represent y and the sequence of approximations (x_k)_{k \in \mathbb{N}_0}:

import Data.Ratio

-- define generalized Heron map
g :: Integer -> Rational -> Rational -> Rational
g n y x = (1 % n) * ((fromIntegral (n-1) * x) + (y / (x^(n-1))))

-- initial approximation
x0 :: Rational
x0 = 1 % 1

-- list of approximations of the n-th root of y
roots :: Integer -> Rational -> [Rational]
roots n y | n  > 1 = iterate (g n y) x0
          | n <= 1 = error "Invalid value of n!!"

We load this script into GHCi. Let us now perform some numerical experiments.

__________

Experiment #1

We would like to compute a rational approximation of \sqrt{2}. We have n = 2 and y = 2. Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 2 2)
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
 886731088897 % 627013566048]

We obtain the following rational approximation after five iterations

x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}

which is the same approximation I obtained twice before.

__________

Experiment #2

We would like to compute a rational approximation of \sqrt{5}. We have n = 2 and y = 5. Here is a very brief GHCi session:

*Main Data.Ratio> take 7 (roots 2 5)
[1 % 1,3 % 1,7 % 3,47 % 21,2207 % 987,4870847 % 2178309,
 23725150497407 % 10610209857723]

We obtain the following rational approximation after six iterations

x_6 = \displaystyle\frac{23725150497407}{10610209857723} \approx \sqrt{5}

How good is this approximation? Let us check:

*Main Data.Ratio> (23725150497407 % 10610209857723)^2 - 5
4 % 112576553224922323902744729
*Main Data.Ratio> 4 / 112576553224922323902744729
3.5531377408652764e-26

Not bad! We could now obtain a rational approximation of the famous golden ratio \varphi := \frac{1 + \sqrt{5}}{2}, as follows:

*Main Data.Ratio> let xs = roots 2 5
*Main Data.Ratio> let phis = map (/2) (map (+1) xs)
*Main Data.Ratio> take 7 phis
[1 % 1,2 % 1,5 % 3,34 % 21,1597 % 987,3524578 % 2178309,
 17167680177565 % 10610209857723]

We obtain the following rational approximation after six iterations

\tilde{\varphi}_6 = \displaystyle\frac{17167680177565}{10610209857723} \approx \varphi

How good is this approximation? Let us recall that the golden ratio is one of the solutions of the equation x^2 - x - 1 = 0. Hence, we can check how close to zero \tilde{\varphi}_6^2 - \tilde{\varphi}_6 - 1 is:

*Main Data.Ratio> let phi6 = 17167680177565 % 10610209857723
*Main Data.Ratio> phi6^2 - phi6 - 1
1 % 112576553224922323902744729
*Main Data.Ratio> 1 / 112576553224922323902744729
8.882844352163191e-27

I would say that is fairly close.

__________

Experiment #3

We now would like to compute a rational approximation of \sqrt[3]{10}. We have n = 3 and y = 10. Here is a very brief GHCi session:

*Main Data.Ratio> take 6 (roots 3 10)
[1 % 1,4 % 1,23 % 8,4909 % 2116,55223315303 % 25495981298,
83759169926117983945469262167029 % 38876457805393768546966848104041]

We obtain the following rational approximation after five iterations

x_5 = \displaystyle\frac{83759169926117983945469262167029}{38876457805393768546966848104041} \approx \sqrt[3]{10}

Is this approximation any good? Let us check:

*Main Data.Ratio> let x5 = (roots 3 10) !! 5 
*Main Data.Ratio> x5^3
5876207108075025254603649658428809744645894398375582704
87986732898174481543266076676033790365389 % 58757060813
2677736272189042625165443146909261971794220978628909441
43919908404893015241356540921
*Main Data.Ratio> 5876207108075025254603649658428809744
6458943983755827048798673289817448154326607667603379036
5389 / 587570608132677736272189042625165443146909261971
79422097862890944143919908404893015241356540921
10.000852709004354

This is somewhat disappointing. The rational approximation x_5 is a ratio of two enormously long integers, but we do not have such a small error. Note that x_5^3 - 10 \approx 10^{-3}. We can do some quick error analysis. Let f_{3,10} (x) = x^3 - 10. The 1st order Taylor approximation of f_{3,10} around x_5 evaluated at x = \sqrt[3]{10} is

f_{3,10} (\sqrt[3]{10}) \approx f_{3,10} (x_5) + f_{3,10}' (x_5) (\sqrt[3]{10} - x_5)

and, since f_{3,10} (\sqrt[3]{10}) = 0, we obtain the estimate of the error

x_5 - \sqrt[3]{10} \approx \displaystyle\frac{f_{3,10} (x_5)}{f_{3,10}' (x_5)}

In other words, the output deviation f_{3,10} (x_5) - f_{3,10} (\sqrt[3]{10}) should be divided by the slope f_{3,10}' (x_5) to obtain an estimate of the magnitude of the approximation error x_5 - \sqrt[3]{10}. Let us use GHCi yet once again:

*Main Data.Ratio> let error = (x5^3 - 10) / (3 * x5^2)
*Main Data.Ratio> error
1670089160826306272530773923851043922672595525468316978
5941152245094153072382174540074985393 % 272741620880842
8590518540423731512512814773109630109912907276686930689
12248623218095571481624481
*Main Data.Ratio> 1670089160826306272530773923851043922
6725955254683169785941152245094153072382174540074985393
 / 2727416208808428590518540423731512512814773109630109
91290727668693068912248623218095571481624481
6.123338108179484e-5

This is still quite a huge error, I would say.

__________

Bonus Experiment

Let us go crazy and attempt to compute a rational approximation of \sqrt[100]{100}. We have n = 100 and y = 100. We expect the rational numbers to quickly become ratios of astronomically long integers. Therefore, we will not be showing any rational approximations on GHCi. Here is something scary:

*Main Data.Ratio> let x3 = (roots 100 100) !! 3
*Main Data.Ratio> let error = (x3^100 - 100) / (100 * x3^99)
*Main Data.Ratio> let error_float = fromIntegral (numerator error) 
                                  / fromIntegral (denominator error)
*Main Data.Ratio> error_float

After some 20 minutes without an output, I decided to abort this experiment. If you want to show x3, you can, but you will see thousands of digits. Thousands! You have been warned.

Tags:, , , , , , , ,
Posted in Haskell, Numerical Methods | 4 Comments »

Heron’s method using integer arithmetic

November 22, 2012

Suppose we are given a positive rational number y, whose square root \sqrt{y} we would like to compute. To be more precise, we would like to compute an approximation of \sqrt{y}.

Let x_k \neq 0 be an approximation of \sqrt{y}. Hence, we have that \sqrt{y} is the geometric mean of x_k and y / x_k. If we replace the geometric mean of these two numbers by their arithmetic mean, we hopefully obtain a “better” approximation x_{k+1}, as follows

x_{k+1} = \displaystyle\frac{1}{2} \left( x_k + \frac{y}{x_k}\right)

If all goes well, then x_N will be “sufficiently close” to \sqrt{y} for a “sufficiently large” integer N. This algorithm is the famous Heron’s method, devised by Heron of Alexandria some 2000 years ago. Some call it the Babylonian method, which suggests that the algorithm may be some 4000 years old. It’s not merely vintage, it’s antique.

Nonetheless, this is not a post on archaeology. We don’t want to sit around and talk about Heron’s method, we want to actually implement it. What programming language should we use? Let us use Haskell. How are we going to represent y and the sequence of approximations (x_k)_{k \in \mathbb{N}_0}?

__________

Heron’s method using floating-point arithmetic

Though millions of people have already implemented Heron’s method using floating-point arithmetic (I implemented it in C in late 2000), one more implementation would do no harm:

y :: Floating a => a
y = 2

-- define Heron map
g :: Floating a => a -> a
g x = 0.5 * (x + y/x)

-- initial approximation
x0:: Floating a => a
x0 = 1

-- list of approximations
xs :: Floating a => [a]
xs = iterate g x0

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[1.0,1.5,1.4166666666666665,1.4142156862745097,
 1.4142135623746899,1.414213562373095]
*Main> let x5 = xs !! 5
*Main> x5**2 - 2
-4.440892098500626e-16

We thus have the following approximation after five iterations

x_5 = 1.414213562373095 \approx \sqrt{2}

__________

Heron’s method using rational arithmetic

Let us now use rational numbers of arbitrary precision. The following script uses the Data.Ratio library:

import Data.Ratio

y :: Rational
y = 2

-- define Heron map
g :: Rational -> Rational
g x = 0.5 * (x + y/x)

-- initial approximation
x0:: Rational
x0 = 1 % 1

-- list of approximations
xs :: [Rational]
xs = iterate g x0

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
 886731088897 % 627013566048]
*Main> let x5 = xs !! 5
*Main> x5*x5 - 2
1 % 393146012008229658338304
*Main> 1 / 393146012008229658338304
2.5435842395854372e-24

Note that the error is eight orders of magnitude smaller than in the previous case (where we used finite-precision floating-point numbers). We pay for more accuracy with more computation. Thus, after five iterations, we have the following rational approximation

x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}

This is probably not news to you, my dear reader. Last April, I did, in fact, implement Heron’s method under the disguise of Newton’s method (see the appendix).

__________

Heron’s method using integer arithmetic

A rational number y \in \mathbb{Q} can be expressed as a fraction a / b, where a and b are integers and b \neq 0. Note that a / b can (obviously?) be represented as a pair of integers (a,b). Since Haskell has arbitrary-precision integers, we can implement Heron’s method using pairs of arbitrary-precision integers to represent y and the sequence of approximations (x_k)_{k \in \mathbb{N}_0}.

Let y := a / b and x_k := p_k / q_k. The recurrence relation

x_{k+1} = \displaystyle\frac{1}{2} \left( x_k + \frac{y}{x_k}\right)

can thus be rewritten in the following form

\displaystyle\frac{p_{k+1}}{q_{k+1}} = \frac{1}{2} \left( \frac{p_k}{q_k} + \frac{a}{b} \frac{q_k}{p_k} \right) = \frac{b \, p_k^2 + a \, q_k^2}{2 b \, p_k \, q_k}

We then have two coupled recurrence relations

\left[\begin{array}{c} p_{k+1}\\ q_{k+1}\end{array}\right] = \left[\begin{array}{c} b \, p_k^2 + a \, q_k^2\\ 2 b \, p_k \, q_k\end{array}\right]

where a, b, p_k and q_k are integers. Here is a Haskell script that implements this vector recurrence relation:

a, b :: Integer
a = 2
b = 1

-- define Heron map
g :: (Integer,Integer) -> (Integer,Integer)
g (p,q) = (b * p^2 + a * q^2, 2 * b * p * q)

-- initial approximation
p0, q0 :: Integer
p0 = 1
q0 = 1

-- list of approximations
xs :: [(Integer,Integer)]
xs = iterate g (p0,q0)

Let us load this script into GHCi. Here’s a GHCi session:

*Main> take 6 xs
[(1,1),(3,2),(17,12),(577,408),(665857,470832),
 (886731088897,627013566048)]
*Main> let x5 = xs !! 5
*Main> x5
(886731088897,627013566048)

Hence, after five iterations, we obtain the very same rational approximation we obtained in the previous case (where we used rational numbers of arbitrary-precision)

x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}

It should be noted, however, that this implementation using arbitrary-precision integers is different in one very important aspect from the one using arbitrary-precision rational numbers. What is the difference? Let GHCi answer this question:

*Main> import Data.Ratio
*Main Data.Ratio> 2 % 2
1 % 1
*Main Data.Ratio> 4 % 4
1 % 1
*Main Data.Ratio> 3 % 6
1 % 2

The main difference between the two implementations is that when we use arbitrary-precision rational numbers, there is automatic reduction of every fraction to its irreducible form, which requires computing the greatest common divisor (GCD) of the numerator and denominator of each rational number at each iteration. This has a computational cost, of course, but it also has the benefit of avoiding reducible fractions (which are inefficient representations of rational numbers). When we use pairs of arbitrary-precision integers, we do not compute the GCD (although we could, if we wanted to), but we may have to work with pairs of unnecessarily long integers (that correspond to reducible fractions), which is potentially dangerous.

For example, suppose that we replace the first three lines in the script above with the following three lines:

a, b :: Integer
a = 2000
b = 1000

We run the script again. Here’s yet another GHCi session:

 
*Main> take 6 xs
[(1,1),(3000,2000),(17000000000,12000000000),
(577000000000000000000000,408000000000000000000000),
(665857000000000000000000000000000000000000000000000,
470832000000000000000000000000000000000000000000000),
(8867310888970000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000,
62701356604800000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000)]
*Main> let x5 = xs !! 5
*Main> x5
(8867310888970000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000,
62701356604800000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000)

Frankly, I am shocked! Multiplying a and b by 1000 results in p_5 and q_5 being multiplied by 10^{90}. I may have made a mistake counting all those zeros, but it’s a lot more zeros than I expected.

__________

Appendix: from Heron to Newton

Interestingly, Heron’s method can be derived from the no less famous Newton’s method. Let us introduce a function f : \mathbb{R} \to \mathbb{R}, defined by f (x) = x^2 - y. By construction, the zeros of function f are \sqrt{y} and -\sqrt{y}. Recall that Newton’s method uses the recurrence relation x_{k+1} = g ( x_k ), where function g : \mathbb{R} \to \mathbb{R} is defined by

g (x) = x - \displaystyle\frac{f (x)}{f' (x)}

where f' is the first derivative of f. Hence, we obtain the following

g (x) = x - \displaystyle\frac{f (x)}{f' (x)} = x - \displaystyle\left(\frac{x^2 - y}{2 x}\right) = \displaystyle\frac{1}{2} \left( x + \frac{y}{x}\right)

which is the recurrence relation used in Heron’s method.

Tags:, , , , , , , ,
Posted in Haskell, Numerical Methods | 3 Comments »

Towards the Cantor set

August 3, 2012

Today we will “construct” the standard (ternary) Cantor set [1]. We start with the closed unit interval \mathcal{C}_0 := [0,1], and then remove its open middle third ]1/3, 2/3[ to obtain \mathcal{C}_1 := [0, 1/3] \cup [2/3, 1]. We then remove the open middle thirds of each of these two intervals to obtain

\mathcal{C}_2 := [0, 1/9] \cup [2/9, 3/9] \cup [6/9, 7/9] \cup [8/9, 1]

and, continuing this process ad infinitum, we obtain a nested sequence of compact sets \mathcal{C}_0 \supset \mathcal{C}_1 \supset \mathcal{C}_2 \supset \dots. Note that \mathcal{C}_n is the union of 2^n intervals, each of length 3^{-n}. We refer to the following set

\mathcal{C} := \displaystyle\bigcap_{n=0}^{\infty} \mathcal{C}_n

as the standard (ternary) Cantor set. This set is most interesting, indeed, since it is uncountable and has Lebesgue measure zero [1].

__________

In Haskell

In Haskell, a closed interval [a, b] can be represented by an ordered pair (a,b). Each set \mathcal{C}_n can be represented by a list of ordered pairs, where each pair represents a closed interval. We create a function f that takes \mathcal{C}_n and returns \mathcal{C}_{n+1} = f ( \mathcal{C}_n ). We will work with arbitrary-precision rational numbers, not floating-point numbers.

The following Haskell script lazily generates the sequence of sets:

import Data.Ratio

type Interval = (Rational, Rational)

-- remove the middle third of an interval
removeMiddleThird :: Interval -> [Interval]
removeMiddleThird (a,b) = [(a,b'),(a',b)] 
                     	  where b' = (2%3)*a + (1%3)*b
                                a' = (1%3)*a + (2%3)*b

-- define function f
f :: [Interval] -> [Interval]
f intervals = concat $ map removeMiddleThird intervals

-- create list of sets
sets :: [[Interval]]
sets = iterate f [(0,1)]

-- define Lebesgue measure
measure :: Interval -> Rational
measure (a,b) = b - a

Note that we used function iterate to generate the sequence of sets. Here is a GHCi session:

*Main> -- take first 4 sets
*Main> take 4 sets
[[(0 % 1,1 % 1)],[(0 % 1,1 % 3),(2 % 3,1 % 1)],
[(0 % 1,1 % 9),(2 % 9,1 % 3),(2 % 3,7 % 9),(8 % 9,1 % 1)],
[(0 % 1,1 % 27),(2 % 27,1 % 9),(2 % 9,7 % 27),(8 % 27,1 % 3),
(2 % 3,19 % 27),(20 % 27,7 % 9),(8 % 9,25 % 27),(26 % 27,1 % 1)]]
*Main> -- compute measure of C_0
*Main> sum $ map measure (sets !! 0) 
1 % 1
*Main> -- compute measure of C_1
*Main> sum $ map measure (sets !! 1) 
2 % 3
*Main> -- compute measure of C_2
*Main> sum $ map measure (sets !! 2) 
4 % 9
*Main> -- compute measure of C_3
*Main> sum $ map measure (sets !! 3) 
8 % 27

This is arguably the most useless Haskell script ever written.

__________

References

[1] Charles Chapman Pugh, Real Mathematical Analysis, Springer-Verlag, New York, 2002.

Tags:, ,
Posted in Haskell, Mathematics | Leave a Comment »

Newton’s method using rational arithmetic

April 2, 2012

The floating-point number system puts further restraints on what can be expected. We can at best hope to keep the relative error under control, and we cannot expect to find zeros far from the origin with great absolute accuracy.

– Richard W. Hamming [1]

I am somewhat disenchanted with floating-point arithmetic. We all have implemented numerical methods using floating-point numbers. But, how many of us have implemented numerical methods using rational numbers of arbitrary precision? I certainly have not. Hence, let us give it a try. In this post I will implement Newton’s method [1] (also known as the Newton-Raphson method) in Haskell using arbitrary-precision rational numbers.

Given a continuous (real) function f : \mathbb{R} \to \mathbb{R} and x_0, an initial approximation of the zero of function f, we then obtain the next approximation using the following recurrence relation

x_{k+1} = x_k - \displaystyle\frac{f (x_k)}{f' (x_k)}

where f' is the first derivative of f. We will use arbitrary-precision rational numbers (represented as a ratio of two Integer values)

type Rational = Ratio Integer

which are defined in the Data.Ratio library.

__________

Example

Suppose we would like to find a rational approximation of \sqrt{2}. We create a real function defined by f (x) = x^2 - 2, whose zeros are  \sqrt{2} and -\sqrt{2}. Note that the first derivative of f is f' (x) = 2 x. The initial approximation is, say, x_0 = 1 / 1, which is closer to \sqrt{2} than to -\sqrt{2}. Hopefully, the sequence \{x_k\} will converge to \sqrt{2} instead of converging to -\sqrt{2}. We are performing a numerical experiment, and success is not guaranteed!

The following Haskell script implements Newton’s method:

import Data.Ratio

-- define function f
f :: Rational -> Rational
f x = x*x - 2

-- define 1st derivative of function f
f' :: Rational -> Rational
f' x = 2*x

-- define Newton map
g :: Rational -> Rational
g x = x - (f x / f' x)

-- initial approximation
x0 :: Rational
x0 = 1 % 1

-- list of rational approximations
xs :: [Rational]
xs = iterate g x0

where we used (higher-order) function iterate to create the list of rational approximations. We run the script above in GHCi and use take to compute the first few rational approximations in the list:

*Main> take 6 xs
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
886731088897 % 627013566048]

Hence, we have, for instance, the following rational approximation

x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}.

How good is this rational approximation? Let us see:

*Main> let x5 = xs !! 5
*Main> x5*x5 - 2
1 % 393146012008229658338304
*Main> 886731088897 / 627013566048
1.4142135623730951
*Main> sqrt 2
1.4142135623730951

Fairly good indeed! And all it took was five iterations!

__________

References

[1] Richard W. Hamming, Numerical Methods for Scientists and Engineers, Dover Publications, 1973.

Tags:, , , , , , ,
Posted in Haskell, Numerical Methods | 12 Comments »

Fair play

March 31, 2012

Last February’s Ponder This challenge was the following:

To decide who would play first in a human vs. computer game, the human suggested to flip a coin it produced from its pocket, stating “heads you win, tails I win.” The computer replied, “Humans are very devious. I know that your coin is not fair, and that the probability of heads on that coin is slightly less than one half. To make this perfectly fair, flip the coin repeatedly, and I win if the coin first shows heads on toss number 1, 14, 15, 19, 20, or 23. You win if it first shows heads on any other toss.”

If the computer is correct and p is the probability that the coin shows heads, give the best rational approximation to p with a denominator having fewer than 10 digits.

Here is the solution.

__________

My solution #1:

Let \mathbb{P} (H) = p be the probability that the coin will land with the heads side up. If the coin is biased in favor of the human player, then p < 1/2. The probability that the coin toss will first produce “heads” on the n-th trial (where n \geq 1) is (1-p)^{n-1} p. Let us define \mathcal{I} := \{1, 14, 15, 19, 20, 23\}. The probability that the computer will win is thus

\displaystyle \sum_{n \in \mathcal{I}} (1-p)^{n-1} p.

We now assume that the computer knows the exact value of p and, thus, is able to compute its probability of winning. Since the computer is interested in making this game “perfectly fair”, then his probability of winning should equal the human’s probability of winning, i.e., it should equal 1/2. Hence, we obtain the equation

p + (1-p)^{13} p + (1-p)^{14} p + (1-p)^{18} p + (1-p)^{19} p +(1-p)^{22} p = \frac{1}{2}

which we can solve using the following Python + SymPy script:

from sympy import *

# create symbolic variable
p = Symbol('p')

# solve polyomial equation
f = p + p*(1-p)**13+p*(1-p)**14+p*(1-p)**18+p*(1-p)**19+p*(1-p)**22
Solutions = solve(f-0.5,p)

# find (purely) real solutions
RealSols = filter(lambda z: im(z)==0, Solutions)

# print solutions
print "\nAll solutions:"
for sol in Solutions:
    print sol

# print real solutions
print "\nReal solutions:"
for sol in RealSols:
	print sol

This script outputs the following list of solutions:

All solutions:
0.334938061741437 - 1.03291350457235*I
0.334938061741437 + 1.03291350457235*I
1.97217279050061 + 0.143849563303172*I
1.50662666735084 - 0.797180194029411*I
1.92419089868505 + 0.454390584787364*I
0.0780115260218319 + 0.464383156442821*I
0.00197767511752256 - 0.101887310197153*I
1.27721993255011 + 0.951650156864179*I
1.27721993255184 - 0.951650156863701*I
0.247036637335883 - 0.624592516289942*I
0.625456509765409 - 0.845921918004056*I
0.0780115260218319 - 0.464383156442821*I
0.00197767511752256 + 0.101887310197153*I
0.991590970537258 - 0.923668626671437*I
1.79082570892033 - 0.599036471538201*I
1.92419089870934 - 0.454390584808249*I
1.79082570891289 + 0.599036471546724*I
0.499905242928390
0.6254565097654 + 0.845921918004058*I
0.991590970537348 + 0.923668626671824*I
0.247036637335883 + 0.624592516289942*I
1.50662666734992 + 0.797180194030292*I
1.97217279048996 - 0.14384956326771*I

Real solutions:
0.499905242928390

Hence, we have that the probability of heads is

p = 0.499905242928390.

We can convert this number to a ratio of integers in Haskell using function approxRational in the library Data.Ratio:

import Data.Ratio

-- probability of heads
p = 0.499905242928390

-- create (infinite) list of epsilons
epsilons = map (10**) [-1,-2..]

-- create (infinite) list of ratios
ratios = map (approxRational p) epsilons

-- list of ratios whose denominators have less than 10 digits
ratios' = takeWhile (\x -> denominator x < 10^10) ratios

Running the script above in GHCi, we obtain a list of ratios:

*Main Data.Ratio> ratios'
[1 % 2,1 % 2,1 % 2,1 % 2,2386 % 4773,2611 % 5223,2636 % 5273,
2638 % 5277,13189 % 26383,44843 % 89703,166183 % 332429,
393036 % 786221,1405961 % 2812455,6243733 % 12489833,
14906352 % 29818355,69693988 % 139414397]

and we finally obtain the following rational approximation

p \approx \displaystyle\frac{69693988}{139414397}

whose denominator has nine digits. According to the webmaster, the rational approximation I presented above is not the best one. Hence, let us now try to obtain a better one.

__________

My solution #2:

Most likely, the reason why I did not obtain the best rational approximation is that the real solution I obtained was not precise enough. Hence, let us use the following MATLAB command to solve the polynomial equation:

solve('p + p*(1-p)^13+p*(1-p)^14+p*(1-p)^18+p*(1-p)^19+p*(1-p)^22 = 0.5')

This command produces the following list of solutions:

ans =

 .19776751175225606651411342151624e-2-.10188731019715301414082803578298*i
 .19776751175225606651411342151624e-2+.10188731019715301414082803578298*i
 .78011526021831885333432975572559e-1-.46438315644282071005250209738747*i
 .78011526021831885333432975572559e-1+.46438315644282071005250209738747*i
    .24703663733588276238861533742100-.62459251628994218281562354280746*i
    .24703663733588276238861533742100+.62459251628994218281562354280746*i
    .33493806174143698737516495548579-1.0329135045723528674007841220342*i
    .33493806174143698737516495548579+1.0329135045723528674007841220342*i
                                        .49990524292838950452962194522516
    .62545650976540488339104798827398-.84592191800405765193777444233924*i
    .62545650976540488339104798827398+.84592191800405765193777444233924*i
    .99159097053730697830960408610308-.92366862667162608535284662817351*i
    .99159097053730697830960408610308+.92366862667162608535284662817351*i
    1.2772199325509563050241902871670-.95165015686389448554291677056666*i
    1.2772199325509563050241902871670+.95165015686389448554291677056666*i
    1.5066266673486985232104204448828-.79718019402868789409470957216463*i
    1.5066266673486985232104204448828+.79718019402868789409470957216463*i
    1.7908257089204728277307542114304-.59903647154154201043833983995659*i
    1.7908257089204728277307542114304+.59903647154154201043833983995659*i
    1.9241908987061350990877558129099-.45439058479360516453156947604463*i
    1.9241908987061350990877558129099+.45439058479360516453156947604463*i
    1.9721727904901564352190617939256-.14384956327517564484908725760533*i
    1.9721727904901564352190617939256+.14384956327517564484908725760533*i

We now obtain a much more precise real solution

p =0.49990524292838950452962194522516.

We modify the Haskell script, using the new value of p and a much finer list of epsilons (the precisions):

import Data.Ratio

-- probability of heads
p = 0.49990524292838950452962194522516

-- create (infinite) list of epsilons
epsilons = map ((1.005)**) [-1,-2..]

-- create (infinite) list of ratios
ratios = map (approxRational p) epsilons

-- list of ratios whose denominators have less than 10 digits
ratios' = takeWhile (\x -> denominator x < 10^10) ratios

Which yields the following rational approximation

p \approx \displaystyle\frac{61031369}{122085875}

whose denominator has nine digits. Although this rational approximation does not match the one in the official solution, I still got my name on the “people who answered correctly” list. Surprisingly, the value of p which I obtained using MATLAB is quite different from the value of p used by the webmaster

p = 0.499905242928506776611191007553.

Which value of p is “more correct”? I do not know. Do you?

Tags:, , , , , ,
Posted in Haskell, Puzzles, Python | Leave a Comment »

Follow

Get every new post delivered to your Inbox.

Join 77 other followers