Well, I think there is no whole or circular argument anymore, but
it will depend on what you are using as defintion and properties
for exp(x) and/or a^x.

1. Version: like I did above, l'Hospital filling up arguments after
Alex's comments. The only thing missing there is exp(x) ' = exp(x).
That follows if you know the (absolute convergent) series for exp
and you know/accept differentiation of such series, since for term
n you have ( x^(n+1)/(n+1)! )' = (n+1)*x^n/(n+1)! = x^n/n!.

2. Version: if you are willing to see the task as difference version
for f'(x) in x=0, f:=a^x, you can even omit the dubious 'Hospital
and use the argumentation given there after that.

3. Version: I looked that up in the nice little book Otto Forster "Analysis1".

He introduces exp as series (and by quotient criterion it is seen absolute
convergent), then gives remainder estimate by geometric series, (10 lines)
error cutting after the N-th term <= 2*|x|^(N+1)/(N+1)!, for |x| <= 1+N/2.

Then |exp(x) - (1+x)| <= x^2 for |x |<= 1 + N/2, N=1 using that estimate.

Now dividing by 0 < |x| < 1 gives what you need

|  e^x - 1     |   | e^x - (1+x) |
|  ------- - 1 | = | ----------- |  <= |x|
|     x        |   |     x       |
 

  '[numer(f1),denom(f1)]': 'eval(%,h=0)': '%' =%;

                [numer(f1), denom(f1)]|      = [0, 0]
                                      |h = 0

  '[diff(numer(f1),h),diff(denom(f1),h)]': 'eval(%,h=0)': '%' =%;

           d             d            |
          [-- numer(f1), -- denom(f1)]|      = [ln(a), 1]
           dh            dh           |h = 0

so it is ln(a) by l'Hospital's rule.

Over time I simplify got used to the fact, that Maple does not always
allow to input things, as one would it expect from typing Mathematics.

Even if it mainly allows much intuitive handling - much better then
what I have seen from Mathematica (where I find it quite odd to write
Exp[x] for exp(x) (exp is lower case in Math and functions have round
brackets - besides the uncomfortable handling on European keyboards).
Or Pari or others (but my experience is rather limited on other CAS). 

And that handling for me was on major point to decide for Maple (after
using Mupad, which handled similar).

Part of 'exceptions' are indexed variables, so I simply try to avoid
them as far as possible.

In many cases x(n) instead of x[n] would do it as well (and for Math
it is the same, indexing more or less just means to use functions on
integer intervals)


Now for your problem (where the 'trick' with x(0) would not work either).

If you do not use the arrow notation but 'unapply' then it works and
you even can see, how Maple beats it into its needed shape:

  42;                   # write down the function expression
  y:=unapply(%, x[0]);  # make it a function by telling it the variable(s)

                            y := x_0 -> 42

One sees, that Maple constructs a function, but replaces x[0] by a symbol
(and thus silently solved its problem with x[0] internally).

Which mathematically is correct: by using the Math notation x[0] has to
be understood as bound variable (in some domain), but of course that does
not depend on the name of the variable.

Now you can use it, it takes indexed variables as input:

  'y( x[0] )': '%'=%;
                             y(x[0]) = 42




As others already mentioned neither x[0] nor x(0) are considered as symbols
(and that's what the error message for the intially entered command wants
to tell us). 

  x[0]; whattype(%); type(%%,symbol);

                               indexed
                                false


  x(0); whattype(%); type(%%,symbol);

                               function
                                false

However that's not what I would expect and my guessing is, that here the
type of x is given, not that of x[0] or x(0).

At least the help shows type( f(x,y,z), function ) as 'true', while in
Math I would say type( f, function ) = true.

But for all of the following ways one gets 'false':

 type( sin, function );

 X:= x -> x; type( X, function );

 x; X:= unapply(%,x); type( X, function );


However using 'type( %, procedure )' gives 'true'.


Which is not quite consistent with intuitive Math notation ... but
usually I can live with those handicaps ...

If computing diagonalizations for symmetric matrices: does Maple go through Eigenvalues/Eigenvectors only (needs roots of polynomials of may be higher order) or does it use other ways (have not tested it beyond dim=3 or 4 in concrete examples)?

try to either use 'formatted' in posting or just put it in a file (prefered: *.mws) and upload using the almsot invisible green arrow at the end of the menu bar of the board editor (visible, if logged in)

if you would use common notations, then your f = 1/(cos(x)^2),
so at Pi/2 ... hence you should not write down the integral

something like

for k  from 1 to 100 by  5 do
  cat(convert(k,string),".jpg");
en do;

Digits:=14;

numQ := proc(y, t1, t2, t3, t4)
local lnConst, e1, e2, p, lnc0, c0, c1, c2, aleph, q, 
  result, t12, t13, t15, t19, t22, t6, t7, t9, lnK5, ln_aleph;
  lnConst := -4681.92777205978276;
  e1 := 9.10323944425767683;
  e2 := 6.56420567208544207;
  p := 945.427297455447429;
  lnc0 := -5.34824719351398036;
  c1 := 16.5500650823099687;
  c2 := .763123640889490148e-1;
  t6 := .217003416872837899e-1*t3;
  t7 := .77256843031344e-1*t4;
  t9 := exp(-.22850024468343e-1*t1+t6-t7+e1);
  t12 := .908028439918035485*t3;
  t13 := .184630769011192971e-2*t4;
  t15 := exp(-.74140364919166*t2+t12+t13+e2);
  t19 := exp(.228500244683431746e-1*t1+t6-t7+e1);
  t22 := exp(.741403649191661292*t2+t12+t13+e2);
  aleph := y*t9+y*t15+t19+t22;
  ln_aleph := ln(aleph);
  lnK5 := lnc0+ln(y)*c1+ln_aleph*p-t4*t4-t2*t2-t3*t3-t1*t1-
    47.7572820906515031*ln_aleph*ln_aleph-ln(-y+1.0)*c2;
  q := lnConst+lnK5;
  result := exp(q);
  return result
end proc

# plot an example
plot('numQ'(0.1,tau,tau,tau,tau), tau= -2 .. 2, title="y=0.1");

The original task should be done by integrating Q, t in IR^4,
and y in ]0 ... 1[, for which t in ]-6 .. 6[^4 might be enough,
may be up to 8 is better. Have not checked it. And not even
tried to think about error estimates seriously. But claim it
is OK ... well, I am not forced to publish ...

  gc();
  infolevel[`evalf/int`]:=5:
  st:=time():
  
  yRange:= 0 .. 1;
  tRange:= -8 .. 8;
  relError:= 1.0 * 10^(-3);
  
  [t1 = -infinity..infinity, t2 = -infinity..infinity, 
    t3 = -infinity..infinity, t4 = -infinity..infinity]:
  subs(-infinity= op(1,tRange), infinity= op(-1,tRange),%):
  Z1:=op(%):
  
  Int('numQ6', [yRange, tRange$4], method=_cuhre, epsilon = relError);
  evalf(%);
  
  `seconds needed` = time()-st;
  infolevel[`evalf/int`]:=0:

It needs ~ 13 sec on my PC and ~ 10 min for relError:= 1.0 * 10^(-6),
where in the latter case 57704829 evaluation points have been used,
6*10^7 function evaluations ... which is a lot.

No better timings if 'compile' is used, the integrator does not like it.


Just a summary on a fairly non-trivial task:

The above comes through processing the original task symbolically
for z and the last 4 variables, hence only 5 variables of the 10
survive. After that a 4-dim transform is found, which makes it an
integral ~ f(y,t)*pdfN(4,t), the later a 4-dim normal pdf with all
correlations = 0, hence centered at 0 (otherwise ~ (-10,100,0,0)
should be choosen, depending on choosing integration order) and it
dominates f for large t's. And finally all can be brought into the
posted procedure using codegen/optimze after resolving overflow issues.


That is my last post on that topic & thread.

Cheers

To me it looks as a system of quadrics and generically they intersect
discrete due to dim = #equations. But the only thing I can 'see' or
guess is that the first 5 are equivalent to "Mittelpunktsflächen"
(= "centered hypersurfaces" ?) and the others seem parabolic. Not sure,
but if one would have the normal form for them, then it may be clearer.

The reason for the question might certainly help ... sqrt(3) smells
a bit like something based on Algebra or Number Theory.

And if the solution is not a discrete set, then it might matter over
which ring this is considered (not that I want to solve it ... never
used that corner of Maple).

So what are your findings so far?

Do you get 'reasonable' results, if you use ~ 0.0001398 as result? I think, it could/should be of that size.

showstat(`evalf/exp`) or showstat(`exp`)shows the code for floats or symbolics

May be, that 'tryhard' has problems with uninitialized variables?
I do not know, but simplify/size is quite good - I included that.

You do not always need 'tryhard': if it is OK without than it is
OK. And do not forget, that finally your (Fortran) compiler will
also do an optimization, so it will be altered in any case.

Also the dimensions of arrays FORTR and stampa do not fit for
assigning, just correct it.

The uploaded sheet contains the modifications, searching for my
initials "AVt" shows some comments.

For the question "codegen or CodeGeneration?": if Maple kills all
the old stuff (at least from the help pages), then it really should
provide equivalents - and they dont (makeproc). What sloppiness.

Hope it helps ...

Download 102_codegen_Fortran_Kram.mws
View file details

More or less - if one wants it numerical (since there seems to be
no easy chance for an analytic solution ... but it is only a model,
to be justified, so perhaps a modification would do as well ...),
then more has to invested (due to numerical cancellation problems
through extreme magnitudes).

Here is a sketch: one can do the integration wrt z (Maple seems to
hang up in simplifications, due to concrete figures), contributing
a multiplicative constant (in terms of Gamma).

The result is like y^r*(-y+1)^(-s)*F*E, where E=exp(quadric in x)
and F is the rest (powers + log).

Pulling out the constant of the quadric it can be seen to be affine
equivalent to sum( x_i ^2 ) + some constant (i.e. a sphere).

Thus (up to constants) we have y^r*(-y+1)^(-s)*F*exp(sum( x_i ^2 ))
by a multi-dimensional change of coordinates in x (note that one
need not to care for domains, since affines map IR^n onto itself
and the needed det(Jacobian) is a constant [y is not affected]).

The last term is a product of normal distributions, hence centered
at zero - which gives a guess, where to integrate: the normal pdf's
should dominate the other terms.

Now carefully collecting all constants *may* give a chance, even if
it seems the integrand is quite small ( << 1).

But I do not want to carry out that plan. And may be some further
simplifications are possible (looking at 'F' and y and log(sum exp)).

However my objection is, that if it is already that difficult to
find the normalizing constant (and that might be somewhat incorrect),
then a subsequent use has a chance to be numerical difficult as well.


PS: first though the poster is actually 'Mike' ... that would be
typical for him ...

1st: avoid floats like C=1.000 if you have no real reason, try exact values.

2nd: try do use functions, it gives better structure

3rd: if having doubts and using test values, then try simple ones or just one
(which can be typed in without running a script).

4th: avoid 'displayprecision' if there is no need for it, try evalf[3]().

Hm ... and using Array (upper case) is more concurrent ...


Now again:

  a := 1: C := 1: L := 4: CS := C*L: # avoid floats ...
  
  # your formula:
  Product(
    GAMMA(1+CS)*GAMMA(1+CS+(L-n+1)*r)*Hypergeom([1+CS, 1], [1+C], r)/
      (GAMMA(r+1+CS)*GAMMA(1+CS+(L-n)*r)), n = 1 .. L-1)*
    Hypergeom([1+CS, 1], [1+C], r);  
  #simplify(%); # not really needed, can be done later ...
  value(%);
  convert(%,StandardFunctions); 
  factor(%);
  f:=unapply(%,r);

                             4   2           4
                   54 (r - 2)  (r  - 2 r + 2)  GAMMA(5 + 4 r)
         f := r -> ------------------------------------------
                                    16             4
                            (-1 + r)   GAMMA(r + 5)


Your values:

  theR:=Array(1 .. 10):
  
  ind := 1:
  for rangedb from -5 by 5 to 40 do 
    y := 10^((1/10)*rangedb); 
    R := C/y; 
    r := R/(a+L*R); #r:='r':
    theR[ind] := r; 
    ind := ind+1 
  end do:
  r:='r':

Do not forget to reset your variables ... I would do it another way.

  theR:=evala(theR);

This gives the exact values for your various r. For all of them you
claim that different results are to be given. So take one, the 2nd
is quite simple, it is r = 1/5.

  f(1/5);        # strange, what Maple does here? No, due to Gamma.
  simplify(%);
  evalf[200](%): # use high precision in case of doubts
  evalf(%);      # cut down do working precision

                             11.07390057

Thx for the paper, however I do not see your formula and the notations
in the article are quite technical, and I do not want to step into the
details there. Well, I do not even see, which coefficients and inputs
you have to use F_A in eq (1) of the paper.

However I am sure your code is not correct, i.e. Product()* 2F1 is not
what you want.


For the convergence problem (approaching abs = 1): this already occures
if working with hypergeometric series, thus Maples does not use it in
a direct, simple way.