I have no reasonable idea, but for dim=6 at least I can get integral = 2.744072e-003 with relative error 1e-6 (#evals = 684,228,711), not using Maple. Is that 'correct'?

https://en.wikipedia.org/wiki/7-Zip, for any MS Win

Es ist eine Reihe, keine Summe einer Reihe (oder eben die Summe über eine Folge).

Es gibt keine "symbolische" Lösung für dieses Problem. Daher "evalf", um ein numerisches Ergebis zu erhalten.

trials and patients ... essentially it is there ... and lost patients for the curve options (thx vv Kitonum!)

Nice, since it also gives insight if being printed. I modified it a bit for 
my personal taste. What I miss: how to modify options for the path (like the
color and the thickness, the default is a bit ugly)?

F:= t -> <cos(2*Pi*t), sin(2*Pi*t)>; # vector, unit interval

p:=Student[VectorCalculus][TangentVector](F(t), 
  range = 0 .. 1, 
  output = plot, #output = animation, #normalized = true,
  vectors=10, 
  vectoroptions = [shape=arrow, thickness = 2, color=blue, 
                   length= 1e-4, head_length = 1/16, head_width =1/16]):

p0:=plots[pointplot](F(0 + 1e-3*0), 
  style=point, symbol=solidcircle, symbolsize=20, color=blue):
p1:=plots[pointplot](F(1 - 1e-3), 
  style=point, symbol=solidcircle, symbolsize=20, color=magenta):
 
plots[display](p, p0,p1);

Yes, any factor is smaller ... stupid me ...

Hm ... I see that 181 is prime, I do not see why d must be prime and would have modified " irem(n,d)<>0 or not(isprime(d)) "

@Christopher2222 

Yes, that what I meant. If your task is not "fitting" then you may consider to filter them away, a kind of pre-processing. But that depends on what you want.

@Christopher2222

Not sure whether it makes sense for you, but you may look for a somewhat
dated book, Carol Alexander, Market Models and IIRC there have been papers
on her site covering similar for asian rice markets (though it is long ago
I cared for those things).

For your data and desired 'model' I would stare at the data and would guess
the rough period from the peaks x=77 and x~90 and the 'level' from average
between 45 and 29, the shift from 0 also can be done. Then you have some
initail guess for fitting and may be able to set restrictions as well.

That way you can relate data and model at least a bit.

Having a plot like yours I would try a second pass, considering some data
as 'outliers'. At least I would remove them at x=77 and 85, possibly also
at x=79. After that it goes more to 'explain' why these are exceptions.

Of course it is not as sophisticated as Mac Dude's view.

@vv 

Could it be that the recursion may fail to stop in finitely many steps?

My feeling is: if x < 1/3 in base 3 has a series without any coefficient=1
(i.e. only 0 or 2) it might be the case, like x = 1/12 (originally x=1/4).

Numerically - for 'precision' - that is not a point.

A brute way to close those many windows might be the task manager (in MS Win) by killing mserver. Though the practical question also might be: having a dozend only - which one to take ... a better way may be to plot to files, I think it is possible to do it for *.eps.

@vv Thank you. 

Correcting eps to 1/(10^40000*Pi) gives k0 = 83837, keeping my notations.

Due to your hint I can reproduce your result:

Using F(3*('eps')*3^k) = 2^k*F(3*('eps')) for k = k0 - 1 and using your
function for the l.h.s. it gives 3/4 as solution, already for n=3 (and
it needs only 4 recursive calls).

Thus one gets 1 / 2 - sol/eval(2^k, k=k0-1)' / 2 = 1/2 - 3/2^(k0+2),
as you said.


PS: there is a small typo in the last Cf(...), 1/3 - .. instead of 1/2 ?

The 'minimalistic' function Cf gives an approximate answer and depends
on computational precision, as vv said. 

For the example one can estimate an error intervall for the true function:

The function is monotone increasing, F(0)=0, F(1-x)=1-F(x), F(x)=F(3*x)/2
(the latter for x <= 1/3) and thus is determined on 0 <= x <= 1/3.

Estimate F(x) for  x:= 1/3 - eps, eps:= 3/(10^40000*Pi) ~ 1E-4000.

F(x) = F(3*x)/2 = F(1 - 3*eps)/2 = (1 - F(3*eps))/2 = 1/2 - F(3*eps)/2.

Solving, expanding and taking floor for 3*eps = 1/3^k one gets that for
k0 := 83835 one has 1/3^(k0+1) < 3*eps < 1/3^(k0).

Since F(1/3^k) = 1/2^k one gets the estimate

1/2 - delta <= F(x) <= 1/2 - delta/2, delta:= 1/2^(k0+1) ~ 0.7*1e-25237.

Speed is not my point (I know that vv is faster). I say that the Chebyshev formulation allows to reduce work by then known formulae. For the re-stated task of course one can transform back to coordinates used by vv.

Is the mathematical solution known for your task?

Pondering on it the only suggestion I have: using the representation of G(n,k) in terms of Chebyshev polynomials allows to reduce it to the diagonal n=k, since the "rest" can be written in terms of explicite representation (of the Cheb poly) without recursion. And for the diagonal a solution is known.

I do not say that a CAS will do it in any case and for high values or far of the diagonal - it is a suggestion. For me it is just unclear why a CAS can provide a solution here.