@Markiyan Hirnyk 

The coefficients of F(n,a,1) and F(n,a,b) seem to be the same, so it should be easy to obtain them.
 

@Markiyan Hirnyk 

It was Maple + "incomplete induction":

with(LinearAlgebra):
F:=(n,a,b) -> Matrix(n, (i, j) ->binomial(a*i+b*j, j) );
factor(expand(Determinant(F(6,a,1))));
R:=(n,a)-> 1/n!*Product(1+k*a,k=1..n)*a^(n*(n-1)/2); # guess for b=1

1.  For the convenience here is a compact statement of the theorem:

2. It is remarkable that is is fast but the simplification of eval(Circle, P) to 0  seems to need some work e.g.

simplify(convert(%,exp));

and is not fast at all. Probably there are better transforms to obtain this faster.

Download str.mw

No, it confirms it!

@Markiyan Hirnyk 

@Markiyan Hirnyk 
Consider

f:=z -> sqrt(1/z)*erf(sqrt(z));  # our function
g:=z -> 1/sqrt(z)*erf(sqrt(z)); 

g is an entire function (0 is a removable singularity) and f equals g except on the negative real axis.

f is discontinuous on the negative axis but this is difficult to see in a standard plot (practically it is g which is plotted).

f(0):=1:
A:=Matrix(40,40, (i,j)-> evalf(Re(f((i-20)/10 + I*(j-20)/10)))):
plots:-surfdata(A);

    Probably such "guessing" problems should be declared "mathematically valid" only if they ask for the shortest Maple expression (as in Kolmogorov complexity) and two such minimal expressions do not exist.
    Of course Maple could be replaced with another language, but Maple is the best :-)

Sorry, I don't see any bug.

@Markiyan Hirnyk 

I don't understand what workaround you are talking about.
The only issue is that in the help _C is not mentioned; probably it is considered obvious.

@Markiyan Hirnyk 

@Preben Alsholm 

But the comparison is irrelevant because RK4 has a fixed step while in RK45 it is adaptive (as implemented in Maple).

@one pound 

Formally the code works for q := Int(1/ln(t), t = 0 .. x)   too
but it is mathematically incorrect (the new integrals do not exist!).

@one pound 

You may use integration by parts but for
q := Int(1/ln(t), t = exp(1) ..  x);

and then get rid of the constant (which does not matter).

@rlopez 

The problem with this approach is that all the integrals (except q itself) are divergent (even as Cauchy Principal Value) for x>1.

@Kitonum 

No, it is in the source code of GAMMA.

@Kitonum 

This is by design because GAMMA increases very fast.
If you really want such huge numbers then set (not recommended):

`GAMMA/magic`:=100,100,infinity,50;

and use s(100);