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

@alex_01 : Here's some more explanation on point (ii).

According to (1) (with s=1), the probability that there is exactly one success in k trials with probability p of success in each one is k*p*(1-p)^(k-1).  I call this f(k). 

> f:= k -> k*p*(1-p)^(k-1);
   assume(p>0, p<1);

The maximum of f(k) for positive integers k will be obtained for the least k for which f(k+1) <= f(k).

> simplify(f(k+1)/f(k));

> solve(% <= 1, k);

,

So we want the least positive integer k >= (1-p)/p, i.e. ceil((1-p)/p).

@alex_01 : Here's some more explanation on point (ii).

According to (1) (with s=1), the probability that there is exactly one success in k trials with probability p of success in each one is k*p*(1-p)^(k-1).  I call this f(k). 

> f:= k -> k*p*(1-p)^(k-1);
   assume(p>0, p<1);

The maximum of f(k) for positive integers k will be obtained for the least k for which f(k+1) <= f(k).

> simplify(f(k+1)/f(k));

> solve(% <= 1, k);

,

So we want the least positive integer k >= (1-p)/p, i.e. ceil((1-p)/p).

This might be a bug.  But it might help if you uploaded a file with the function you're trying to convert to Fortran.

gkokovidis's solution is not quite right: if the implicitplot has more than one component, there will be lines with four numbers (corresponding to where one component ends and the other starts).  For example, with

p1:= implicitplot(x^2 - y^2 - 1, x = -2 .. 2, y = -2 .. 2):

one line in the middle of the file is

-2.000000 1.730455 2.000000 1.730455

To avoid that, you might try something like this:

> Q:=[op([1,1..-2], p1)]:
   for i from 1 to nops(Q) do
     fprintf("c:/temp/implotdata.txt","%f",Q[i]);
     fprintf("c:/temp/implotdata.txt","\n");
   end do:
   fclose("c:/temp/implotdata.txt");

gkokovidis's solution is not quite right: if the implicitplot has more than one component, there will be lines with four numbers (corresponding to where one component ends and the other starts).  For example, with

p1:= implicitplot(x^2 - y^2 - 1, x = -2 .. 2, y = -2 .. 2):

one line in the middle of the file is

-2.000000 1.730455 2.000000 1.730455

To avoid that, you might try something like this:

> Q:=[op([1,1..-2], p1)]:
   for i from 1 to nops(Q) do
     fprintf("c:/temp/implotdata.txt","%f",Q[i]);
     fprintf("c:/temp/implotdata.txt","\n");
   end do:
   fclose("c:/temp/implotdata.txt");

Unfortunately isolve doesn't always work even for quadratics.  For example:

> Q := x^2 + 5*x*y + 2*y^2 + 2*x + 3*y;
    isolve(Q = 27);

# returns nothing

Note that x=1, y=2 is a solution.

Unfortunately isolve doesn't always work even for quadratics.  For example:

> Q := x^2 + 5*x*y + 2*y^2 + 2*x + 3*y;
    isolve(Q = 27);

# returns nothing

Note that x=1, y=2 is a solution.

@Chris : You want to get from abs(z-3) <= 2*abs(z+3) to (x + 5)^2 + y^2 >= 16

> ineq:= abs(z-3) <= 2*abs(z+3);
   map(`^`,ineq,2);  # squaring both sides
   evalc(eval(%,z=x+I*y)); # go to xy coordinates, and get rid of absolute values;
  

   
> % - lhs(%);  # subtract left side from both sides
    student[completesquare](%,x); # complete the square using variable x
    %/3 + 16;

@Chris : You want to get from abs(z-3) <= 2*abs(z+3) to (x + 5)^2 + y^2 >= 16

> ineq:= abs(z-3) <= 2*abs(z+3);
   map(`^`,ineq,2);  # squaring both sides
   evalc(eval(%,z=x+I*y)); # go to xy coordinates, and get rid of absolute values;
  

   
> % - lhs(%);  # subtract left side from both sides
    student[completesquare](%,x); # complete the square using variable x
    %/3 + 16;

A simple work-around in this case (though it makes the 2D math pretty) is to use the prefix form `.`(a,b) rather than the infix a . b.

Puzzle: why doesn't the following work?

subs(Typesetting:-delayDotProduct = :-`.`, eval(f));

That's strange: it also works for me today, but I'm sure it didn't when I tried it the other day.
Oh, oh: I hope I didn't just copy the ":" from the original posting. 

That's strange: it also works for me today, but I'm sure it didn't when I tried it the other day.
Oh, oh: I hope I didn't just copy the ":" from the original posting. 

As I said, I suspect the coefficients will be extremely complicated: perhaps so complicated that in practice it is impossible to solve the equations.  I tried this:

eqs := subs(diff(w(r, t), r$2) = Wrr, diff(w(r, t), r) = Wr, diff(u(r, t), r) = Ur, w(r, t) = W, u(r, t) = U,
 map(numer, [T1, T2, T3, T4]));
with(Groebner);
Basis(eqs, plex(Mr, Nr, `Mθ`, `Nθ`));

but after some time (with Memory showing as 761.11M and Time as 776.01s) I got

Error, (in Groebner:-Basis) object too large

IIRC a polynomial is allowed to have up to 2^25 terms, and violating this limit is what's likely causing this error.

I suppose it's possible that this is just an intermediate step and the final results would be much simpler, but I doubt it.  So this raises the question, what would you do with such a gigantic object if you could obtain it?

As I said, I suspect the coefficients will be extremely complicated: perhaps so complicated that in practice it is impossible to solve the equations.  I tried this:

eqs := subs(diff(w(r, t), r$2) = Wrr, diff(w(r, t), r) = Wr, diff(u(r, t), r) = Ur, w(r, t) = W, u(r, t) = U,
 map(numer, [T1, T2, T3, T4]));
with(Groebner);
Basis(eqs, plex(Mr, Nr, `Mθ`, `Nθ`));

but after some time (with Memory showing as 761.11M and Time as 776.01s) I got

Error, (in Groebner:-Basis) object too large

IIRC a polynomial is allowed to have up to 2^25 terms, and violating this limit is what's likely causing this error.

I suppose it's possible that this is just an intermediate step and the final results would be much simpler, but I doubt it.  So this raises the question, what would you do with such a gigantic object if you could obtain it?