I would worry about how much apparent randomness is introduced by the encoding process.  For example, suppose you spoiled the randomness of the lottery by making sure that every draw included the number 25.  Would your tests be able to detect that?

@gotamo : If there's a slight gradient (always with the same sign), you're talking about a single maximum at one of the ends of the interval.  So in the above context you might take the average of the two x[i] values that you know the maximum is between.

@gotamo : If there's a slight gradient (always with the same sign), you're talking about a single maximum at one of the ends of the interval.  So in the above context you might take the average of the two x[i] values that you know the maximum is between.

@hirnyk :  OK, you want Maple to take a picture as input and count the regions.  I like a challenge!  You could start by putting the picture into a Matrix.  Let's say this is a black-and-white image with the lines in black.  Then make a Graph with vertices the white pixels, connected to their nearest neighbours that are white.  Finally, ConnectedComponents in the GraphTheory package would produce the connected components of this. Unfortunately it might not work very well, because pixelation effects may produce some small connected components, especially when two intersecting lines are almost parallel.

For example:

# make a picture
V:= [[1,1],[1,3],[2,5], [5,4], [5,2], [4,6],[3,1]];
P:=plots[display](seq(seq(plot([V[i],V[j]],colour=black),j=i+1..nops(Vertices)),
  i=1..nops(Vertices)-1), axes=none, view=[0.5..5.5,0.5 .. 6.5]):
# write it to a jpeg 
plotsetup(jpeg,plotoutput="c:/test/testpic.jpg",plotoptions="height=240,width=240"); 
P; plotsetup(default);


# now read it in to Maple
R:= ImageTools:-Read("c:/test/testpic.jpg");

                   [ 1..240 x 1..240 x 1..3 3-D Array ]
              R := [ Data Type: float[8]              ]
                   [ Storage: rectangular             ]
                   [ Order: C_order                   ]

# only one colour component needed, say red
R1:= R[1..,1..,1];
# N is the number of rows and columns (in this case the same)
N := 240:
# load GraphTheory, find the vertices (white pixels) and edges.  I'll encode 
# the pixel (i,j) as vertex N*i+j.
with(GraphTheory):
vertices:= map(t -> N*t[1]+t[2],select(t -> (R1[op(t)] > 0.5), 
     [seq(seq([i,j],i=1..N),j=1..N)])):
vedges:= map(t -> {N*t[1]+t[2],N*t[1]+t[2]+1}, 
     select(t -> (R1[t[1],t[2]] > 0.5 and R1[t[1],t[2]+1] > 0.5), 
       {seq(seq([i,j],i=1..N),j=1..N-1)})):
hedges:= map(t -> {N*t[1]+t[2],N*(t[1]+1)+t[2]},
     select(t -> (R1[t[1],t[2]] > 0.5 and R1[t[1]+1,t[2]]>0.5),
       {seq(seq([i,j],i=1..N-1),j=1..N)})):
# construct the graph
G:= Graph(vertices, vedges union hedges);

  G := Graph 1: an undirected unweighted graph with 55362 vertices\
         and 107230 edge(s)

# now find the connected components
C:= ConnectedComponents(G):
nops(C);

                                  61

# Oops: should be 50 I think (49 inside the heptagon plus one outside; note 
# that in this case we have three diagonals intersecting at one interior point).  
# However, there are some very small connected components, 
# of one or two pixels, which I think accounts for the discrepancy.

@hirnyk :  OK, you want Maple to take a picture as input and count the regions.  I like a challenge!  You could start by putting the picture into a Matrix.  Let's say this is a black-and-white image with the lines in black.  Then make a Graph with vertices the white pixels, connected to their nearest neighbours that are white.  Finally, ConnectedComponents in the GraphTheory package would produce the connected components of this. Unfortunately it might not work very well, because pixelation effects may produce some small connected components, especially when two intersecting lines are almost parallel.

For example:

# make a picture
V:= [[1,1],[1,3],[2,5], [5,4], [5,2], [4,6],[3,1]];
P:=plots[display](seq(seq(plot([V[i],V[j]],colour=black),j=i+1..nops(Vertices)),
  i=1..nops(Vertices)-1), axes=none, view=[0.5..5.5,0.5 .. 6.5]):
# write it to a jpeg 
plotsetup(jpeg,plotoutput="c:/test/testpic.jpg",plotoptions="height=240,width=240"); 
P; plotsetup(default);


# now read it in to Maple
R:= ImageTools:-Read("c:/test/testpic.jpg");

                   [ 1..240 x 1..240 x 1..3 3-D Array ]
              R := [ Data Type: float[8]              ]
                   [ Storage: rectangular             ]
                   [ Order: C_order                   ]

# only one colour component needed, say red
R1:= R[1..,1..,1];
# N is the number of rows and columns (in this case the same)
N := 240:
# load GraphTheory, find the vertices (white pixels) and edges.  I'll encode 
# the pixel (i,j) as vertex N*i+j.
with(GraphTheory):
vertices:= map(t -> N*t[1]+t[2],select(t -> (R1[op(t)] > 0.5), 
     [seq(seq([i,j],i=1..N),j=1..N)])):
vedges:= map(t -> {N*t[1]+t[2],N*t[1]+t[2]+1}, 
     select(t -> (R1[t[1],t[2]] > 0.5 and R1[t[1],t[2]+1] > 0.5), 
       {seq(seq([i,j],i=1..N),j=1..N-1)})):
hedges:= map(t -> {N*t[1]+t[2],N*(t[1]+1)+t[2]},
     select(t -> (R1[t[1],t[2]] > 0.5 and R1[t[1]+1,t[2]]>0.5),
       {seq(seq([i,j],i=1..N-1),j=1..N)})):
# construct the graph
G:= Graph(vertices, vedges union hedges);

  G := Graph 1: an undirected unweighted graph with 55362 vertices\
         and 107230 edge(s)

# now find the connected components
C:= ConnectedComponents(G):
nops(C);

                                  61

# Oops: should be 50 I think (49 inside the heptagon plus one outside; note 
# that in this case we have three diagonals intersecting at one interior point).  
# However, there are some very small connected components, 
# of one or two pixels, which I think accounts for the discrepancy.

I have only one thing to add to Pagan's excellent answer: if your i is supposed to be a square root of -1, in Maple that is I rather than i.

I have only one thing to add to Pagan's excellent answer: if your i is supposed to be a square root of -1, in Maple that is I rather than i.

Also interesting is the "tip of the horn", a point in the (r,k) plane where the number of real zeros of the discriminant d of P[3](x) changes (regarding d as a function of r).  That happens at k = sqrt(27), r = sqrt(27)/8.

Axel seems to have missed the line

k :=  1/(x*h11 - 2*y*h1 + z*h);

Thus his solution is for arbitrary k.

There is a problem with your solution, though, because it makes x*h11 - 2*y*h1+z*h, which is supposed to be 1/k, 0:

> normal(eval(x*h11 - 2*y*h1 + z*h, sol));

0

So these are not really solutions, but rather values that make the equations undefined.

Groebner basis methods confirm that (at least for "generic" values of the parameters) there are no solutions that don't make this 0.

> with(Groebner):
   G:= map(numer @ lhs, eqs);
   Solve(G, [x,y,z]);

So to get the numerators zero you either have x*h11 - 2*y*h1+z*h = 0 with z <> 0, or x,y,z all 0 (which again makes x*h11 - 2*y*h1+z*h = 0).

Axel seems to have missed the line

k :=  1/(x*h11 - 2*y*h1 + z*h);

Thus his solution is for arbitrary k.

There is a problem with your solution, though, because it makes x*h11 - 2*y*h1+z*h, which is supposed to be 1/k, 0:

> normal(eval(x*h11 - 2*y*h1 + z*h, sol));

0

So these are not really solutions, but rather values that make the equations undefined.

Groebner basis methods confirm that (at least for "generic" values of the parameters) there are no solutions that don't make this 0.

> with(Groebner):
   G:= map(numer @ lhs, eqs);
   Solve(G, [x,y,z]);

So to get the numerators zero you either have x*h11 - 2*y*h1+z*h = 0 with z <> 0, or x,y,z all 0 (which again makes x*h11 - 2*y*h1+z*h = 0).

Actually, if _C1 = 8*ln(2), we don't know what the limit will be.  It might be anything.  For example:

> Q:= A + (_C1 - 8*ln(2))*exp(t);
   limit(Q,t=infinity);

> eval(Q, _C1 = 8*ln(2));

> limit(%, t=infinity);

Actually, if _C1 = 8*ln(2), we don't know what the limit will be.  It might be anything.  For example:

> Q:= A + (_C1 - 8*ln(2))*exp(t);
   limit(Q,t=infinity);

> eval(Q, _C1 = 8*ln(2));

> limit(%, t=infinity);

@alex_01 : For example, f(1) = p (this is just the probability that the last trial is a success), while f(2) = 2*p*(1-p)  (the probability that there is exactly one success in the last two trials, i.e. either (failure, success) or (success, failure)).  For k = 1 to be optimal you need f(1) >= f(2), which is equivalent to p >= 1/2.  Your coin-flip is p = 1/2.    And indeed for p = 1/2, (1-p)/p = 1.

@alex_01 : For example, f(1) = p (this is just the probability that the last trial is a success), while f(2) = 2*p*(1-p)  (the probability that there is exactly one success in the last two trials, i.e. either (failure, success) or (success, failure)).  For k = 1 to be optimal you need f(1) >= f(2), which is equivalent to p >= 1/2.  Your coin-flip is p = 1/2.    And indeed for p = 1/2, (1-p)/p = 1.

Further to this: that condition k >= (1-p)/p can be written as k*p/(1-p) >= 1.  If the probability of success is p, the odds in favour of success are p to 1-p, or p/(1-p) to 1. 

More generally, suppose the j'th trial (starting from the end) has probability p[j] of success.  Let f(k) be the probability that there is exactly one success in the last k trials.  Then it turns out that the condition to have f(k+1) <= f(k) is
sum(p[j]/(1-p[j]), j=1..k) >= 1.  So, as the article says, you add up the "odds" p[j]/(1-p[j]), starting from the last trial, until the sum exceeds 1. 

This can be seen as follows.  Let Q(k) = product(1-p[j],j=1..k) be the probability that all k trials are failures.  Then to have exactly one success in the last k+1 trials, either trial k+1 is a success and the others are failures, or trial k+1 is a failure and there is one success in the last k.  Thus

f(k+1) = p[k+1]*Q(k) + (1-p[k+1])*f(k)

f(k+1)/f(k) = p[k+1]*Q(k)/f(k) + 1 - p[k+1]
                 <= 1 if and only if Q(k) <= f(k)

But it is easy to see that f(k) = Q(k)*sum(p[j]/(1-p[j]), j=1..k), since the probability that trial number i is a success and the others are failures is obtained from the expression for Q(k) by replacing the factor 1-p[i] by p[i].