So what do *actually* need ? May be it would help to say more
about the context instead of restricting to technical details only

@GPekov : it should fit, 8 decimals is ~ single precision

  Digits:=15;

  # w =XVal, xHalf = X11, c = C12 
  F:=proc(X,                       # variable (complex) 
    arrReturn,                     # array storing results (complex)
    CI0P,CI1P,CK0P,CK1P)           # coefficients for polynomials (Reals)
  local xHalf, w, C11, c, NST, FI0M, FI1M, FK0M, FK0, FK1M, FK1, C13, FK2,
     x, byX, absxSquare, i; 
  #  global CI0P,CI1P,CK0P,CK1P,CK0G,CK1G; 
  
    x:=evalf(X); # in case of symbolic input, like 2*Pi*I
    byX:= 1.0/x; 
    w:= 0.04*x*x;
    absxSquare:= Re(x)*Re(x) + Im(x)*Im(x);
    #if 0.0 = absxSquare then error end if;
    if absxSquare < 25.0 then
      xHalf:= 0.5*x;
      c:=ln(xHalf);
      if   absxSquare < 0.790678161e-9 then NST:= 2; # sqrt(absxSquare) = abs(x) < 0.28119e-4;
      elif 0.7905877225 < absxSquare   then NST:=15; # 0.88915 < sqrt(absxSquare) = abs(x)
      else 
        C11:= 1.6094379124341003746 - ln(absxSquare)*0.5; # = (ln(5/XAbs));
        NST:=ceil(max(2,min(15,24.177/C11)));
      end if; 
  
  ### FI0M:=add(CI0P[i]*w^(i-1),i=2..NST) etc, Horner form for that polynomials
  FI0M:= CI0P[NST]; FI1M:=CI1P[NST]; FK0M:=CK0P[NST]; FK1M:=CK1P[NST];
  for i from NST-1 to 2 by -1 do
    FI0M := w*FI0M + CI0P[i];
    FI1M := w*FI1M + CI1P[i];
    FK0M := w*FK0M + CK0P[i];
    FK1M := w*FK1M + CK1P[i];
  end do;
  FI0M := w*FI0M; FI1M := w*FI1M; FK0M := w*FK0M; FK1M := w*FK1M;
  ###
  
      FI1M:= x*FI1M; 
      FK0M:=-FI0M*c+FK0M; 
      FK1M:=(FI1M*x*c+FK1M)*byX+(CK0P[1]+0.5)*xHalf;
  
      FK0:= FK0M-c+CK0P[1];
      FK1:=FK1M+byX+(c-CK0P[1]-0.5)*xHalf; 
  
    elif 25.0      C13:=(exp(-x)/sqrt(x)); 
  #  FK0:= (CK0G[1]+CK0G[2]*w^(2-1))*C13;
  #  FK1:= (CK1G[1]+CK1G[2]*w^(2-1))*C13;
      FK0:= (1.25331410000000 -.313328530000000e-1*w^(2-1))*C13;
      FK1:= (1.25331410000000 +.939985600000000e-1*w^(2-1))*C13;
    else
      error;
    end if;
    
    FK2:=(FK0+FK1*2*byX);
  
  arrReturn[1]:=FK0;
  arrReturn[2]:=FK1;
  arrReturn[3]:=FK2; 
  end proc; 
  maplemint(F);

Then you define an Array for the 'output' like

  arrResult:=Array([0.0, 0.0, 0.0], datatype=complex[8]);

and after calling the procedure the results are stored there (it avoids, that
in every call an Array has to be constructed), a convenient usage would be

  G:=proc(X) 
  global CI0P,CI1P,CK0P,CK1P, arrResult;
  F(X, 
    arrResult,
    CI0P,CI1P,CK0P,CK1P); 
  return arrResult 
  end proc;

  zTst:=-2.9+I*0; #zTst:=-Pi+exp(1)*I;
  GBIKF012(zTst): evalf(%);
  G(zTst): simplify(%, zero);

Testing (at least in parts) can be done like this:

  checkP:=proc(t) GBIKF012(t)[1] end proc;
  checkV:=proc(t) evalhf(G(t))[1] end proc;

  plot( checkP - checkV, 1e-6 .. 5 - 1e-6);   # absolut differences
  plot( 1 - checkP/checkV, 1e-6 .. 5 - 1e-6); # relative differences

  plot( (Re@checkP) - (Re@checkV), 5 + 1e-6 .. 8);   # absolut differences
  plot( 1 - (Re@checkP)/(Re@checkV), 5 + 1e-6 .. 8); # relative differences


However I have troubles to compile F (even if providing complete informations
about the input Arrays), but I think that I will stop here.

I hope you find some variant which is useful for you.

Pekov_optimize.mws

Edited: I recognize, that this board software still is not able to correctly
treat code, which contains 'less than' etc. Just take the sheet. Sigh &  Amen.

@GPekov : it should fit, 8 decimals is ~ single precision

  Digits:=15;

  # w =XVal, xHalf = X11, c = C12 
  F:=proc(X,                       # variable (complex) 
    arrReturn,                     # array storing results (complex)
    CI0P,CI1P,CK0P,CK1P)           # coefficients for polynomials (Reals)
  local xHalf, w, C11, c, NST, FI0M, FI1M, FK0M, FK0, FK1M, FK1, C13, FK2,
     x, byX, absxSquare, i; 
  #  global CI0P,CI1P,CK0P,CK1P,CK0G,CK1G; 
  
    x:=evalf(X); # in case of symbolic input, like 2*Pi*I
    byX:= 1.0/x; 
    w:= 0.04*x*x;
    absxSquare:= Re(x)*Re(x) + Im(x)*Im(x);
    #if 0.0 = absxSquare then error end if;
    if absxSquare < 25.0 then
      xHalf:= 0.5*x;
      c:=ln(xHalf);
      if   absxSquare < 0.790678161e-9 then NST:= 2; # sqrt(absxSquare) = abs(x) < 0.28119e-4;
      elif 0.7905877225 < absxSquare   then NST:=15; # 0.88915 < sqrt(absxSquare) = abs(x)
      else 
        C11:= 1.6094379124341003746 - ln(absxSquare)*0.5; # = (ln(5/XAbs));
        NST:=ceil(max(2,min(15,24.177/C11)));
      end if; 
  
  ### FI0M:=add(CI0P[i]*w^(i-1),i=2..NST) etc, Horner form for that polynomials
  FI0M:= CI0P[NST]; FI1M:=CI1P[NST]; FK0M:=CK0P[NST]; FK1M:=CK1P[NST];
  for i from NST-1 to 2 by -1 do
    FI0M := w*FI0M + CI0P[i];
    FI1M := w*FI1M + CI1P[i];
    FK0M := w*FK0M + CK0P[i];
    FK1M := w*FK1M + CK1P[i];
  end do;
  FI0M := w*FI0M; FI1M := w*FI1M; FK0M := w*FK0M; FK1M := w*FK1M;
  ###
  
      FI1M:= x*FI1M; 
      FK0M:=-FI0M*c+FK0M; 
      FK1M:=(FI1M*x*c+FK1M)*byX+(CK0P[1]+0.5)*xHalf;
  
      FK0:= FK0M-c+CK0P[1];
      FK1:=FK1M+byX+(c-CK0P[1]-0.5)*xHalf; 
  
    elif 25.0      C13:=(exp(-x)/sqrt(x)); 
  #  FK0:= (CK0G[1]+CK0G[2]*w^(2-1))*C13;
  #  FK1:= (CK1G[1]+CK1G[2]*w^(2-1))*C13;
      FK0:= (1.25331410000000 -.313328530000000e-1*w^(2-1))*C13;
      FK1:= (1.25331410000000 +.939985600000000e-1*w^(2-1))*C13;
    else
      error;
    end if;
    
    FK2:=(FK0+FK1*2*byX);
  
  arrReturn[1]:=FK0;
  arrReturn[2]:=FK1;
  arrReturn[3]:=FK2; 
  end proc; 
  maplemint(F);

Then you define an Array for the 'output' like

  arrResult:=Array([0.0, 0.0, 0.0], datatype=complex[8]);

and after calling the procedure the results are stored there (it avoids, that
in every call an Array has to be constructed), a convenient usage would be

  G:=proc(X) 
  global CI0P,CI1P,CK0P,CK1P, arrResult;
  F(X, 
    arrResult,
    CI0P,CI1P,CK0P,CK1P); 
  return arrResult 
  end proc;

  zTst:=-2.9+I*0; #zTst:=-Pi+exp(1)*I;
  GBIKF012(zTst): evalf(%);
  G(zTst): simplify(%, zero);

Testing (at least in parts) can be done like this:

  checkP:=proc(t) GBIKF012(t)[1] end proc;
  checkV:=proc(t) evalhf(G(t))[1] end proc;

  plot( checkP - checkV, 1e-6 .. 5 - 1e-6);   # absolut differences
  plot( 1 - checkP/checkV, 1e-6 .. 5 - 1e-6); # relative differences

  plot( (Re@checkP) - (Re@checkV), 5 + 1e-6 .. 8);   # absolut differences
  plot( 1 - (Re@checkP)/(Re@checkV), 5 + 1e-6 .. 8); # relative differences


However I have troubles to compile F (even if providing complete informations
about the input Arrays), but I think that I will stop here.

I hope you find some variant which is useful for you.

Pekov_optimize.mws

Edited: I recognize, that this board software still is not able to correctly
treat code, which contains 'less than' etc. Just take the sheet. Sigh &  Amen.

I go with Reiman: roughly evalf can not commute with functions (and would like
to here MORE from acer, may be as very own contribution - please!).

The system works with relative precision (as it was already said), but that
function 'frac' does not work that way - only up to scaling (to ~ 1) and that
is just what you do in an 'inverse' way, you scale off.

By the way (and that is possibly covered by Reiman's post) you would observe
similar, if writing a difference a - b (say a=Pi and b=rational approximate)
of same magnitude, same leading decimal places:

A direct evalf would give zero, while it may be false, one pain in Numerics.


Concerning MMA (and the mentioned discussion):  nobody knows how many 'Digits'
that system uses for a certain task.

so: what is your input (type, range)?

and: it is somewhat puzzling, that your Arrays of
constants have only 8 decimal places ... why?

hello paula,

I use M12 for *.mw sheets (for some other reasons) and get an error
about the alphabeth A - Z

I should have said that in the worksheet the jump in 0 is defined through

_EnvUseHeavisideAsUnitStep := true

I should have said that in the worksheet the jump in 0 is defined through

_EnvUseHeavisideAsUnitStep := true

Here a somewhat more simple version for the problem

data:= [c=1];
U:=Int(Int(Heaviside(t-x-y/c),y = 1+x .. c*x),x = 0 .. 1);
V:= eval(U, data);

u:=value(U):
uData:=eval(u, data); # = symbolic solution of the integral, now evaluated in parameter c=1

v:=value(V); # = first substitute parameter c=1, then evaluate to a symbolic solution 

subs(t=1, V); evalf[8](%); # a test value in t=1, if the c=1 was substituted first
  - 0.25

eval(uData, t=1);
  1/2

eval(v, t=1);
  -1/4

Here a somewhat more simple version for the problem

data:= [c=1];
U:=Int(Int(Heaviside(t-x-y/c),y = 1+x .. c*x),x = 0 .. 1);
V:= eval(U, data);

u:=value(U):
uData:=eval(u, data); # = symbolic solution of the integral, now evaluated in parameter c=1

v:=value(V); # = first substitute parameter c=1, then evaluate to a symbolic solution 

subs(t=1, V); evalf[8](%); # a test value in t=1, if the c=1 was substituted first
  - 0.25

eval(uData, t=1);
  1/2

eval(v, t=1);
  -1/4

evaluating parameters and integration do not commute here ?

Changing variables is ok, I think it is not limited to continous
integrands, but the transformation itself has to be continous.

evaluating parameters and integration do not commute here ?

Changing variables is ok, I think it is not limited to continous
integrands, but the transformation itself has to be continous.

Markiyan Hirnyk pointed out to me, that I assigned tau twice in the sheet,
at the very beginning.

It was tau_1 and tau_t in the original post and I hate 'decorated variables',
but the latter is not used in the sheet, so I did not recognize the double
assignment (the 2nd assignment would be ignored, there is no tau after
the 1st one).

So it still works (despite my sloppiness)

Markiyan Hirnyk pointed out to me, that I assigned tau twice in the sheet,
at the very beginning.

It was tau_1 and tau_t in the original post and I hate 'decorated variables',
but the latter is not used in the sheet, so I did not recognize the double
assignment (the 2nd assignment would be ignored, there is no tau after
the 1st one).

So it still works (despite my sloppiness)

I guess, that Markiyan Hirnyk says: 

It is not a continous function in parameters in general.