#####################################################################
# this file contains some useful procedures for the computation of atan in
# DDE

restart:
Digits := 200:

with (numapprox):with(orthopoly):
read "double-extended.mpl";
read "common-procedures.mpl";

e := 2^(-6.3):
marge := 2^(-30);
P19 := convert(series(arctan(x),x=0,20),polynom):
P18 := (P19-x)/x;
Q18 := polyExact2Ext(P18, 9);
Q9 := -1/3 +1/5*x - 1/7*x^2 + 1/9*x^3 - 1/11*x^4 + 1/13*x^5 - 1/15*x^6 + 1/17*x^7- 1/19*x^8;
Qprime := poly_exact( -x /3 + x^2/5 - x^3/7 + x^4 /9);
#log2 (infnorm( (arctan(x)-x*(1+Q18))/x, x=0..e));

#function taken from the "old" coef_atan to compute interval parameters.
maxx := bi -> simplify( solve( (x-bi) / (1+x*bi) = e ,x) ):
minx := bi -> simplify(solve( (x-bi) /(1+x*bi)=-e ,x)):
nextbi := proc (x) evalf( max(solve( minx(bi) = x ,bi) )*(1-marge)); end:
allbi := proc (n)
local xi,nbi,x,i,j;
global b,a, nbofai, nbofbi, valuefordicho;
    x := e;
    nbi := 0;
    i := 0;
    while(i<n and nbi < 1/e and nbi >= 0) do
        nbi := nearest ( nextbi(x) );
        b[i] := evalf( nbi );
        a[i] := x;
        x   := evalf(maxx(nbi));
        i:=i+1;
    od;
    j:=0;
    while ( 2^j < i ) do j:=j+1 od:
    nbofai := i;
    nbofbi := i;
    b[i-1] := nearest(1/e+4):
    valuefordicho := 2^j;
    return i,b[0],b[i-1];
end:
allbi(100);
#-------------

#-----------------------
#quick error calc :
#-----Reduction:-------
XredEpsilon := 2^(-64) + 2^(-64) + 2^(-64) + 2^(-64):
Xred2Epsilon := 2*XredEpsilon + 2^(-64);
log[2.](Xred2Epsilon);
errlist:=errlist_quickphase_horner(degree(Qprime),0,0,Xred2Epsilon, 2^(-64));
errorr := compute_horner_rounding_error(Qprime,x,e,errlist,true);
qEpsilon := errorr[1]:
log[2.](qEpsilon);
deltaApprox := e^11/11:
polyDelta := e^3*qEpsilon + 2^(-64)*e+ 3*2^(-64)*e^3:
EpsilonFinal := polyDelta/arctan(e)+2*2^(-64)+deltaApprox/e:
log[2.](EpsilonFinal);
#------no reduction-------
Errlist:=errlist_quickphase_horner(degree(Qprime),0,0,0, 2**(-64)):
qEpsilon:= compute_horner_rounding_error(Qprime,x,e,errlist,true)[1]:
deltaApprox := e^11/5:
EpsilonFinalNoRed := infnorm( (qEpsilon.x^3+2^(-64).x+2^(-64).x^3)/arctan(x),x=2^(-27)..2^(-6) ) + deltaApprox + 2^(-64):
log[2.](EpsilonFinalNoRed);
#-----------------------

# Output :

doth:=proc()
 local filename, fd, i, hi, lo;
   filename := "TEMPATAN/atan_ext.h":
    fd := fopen (filename,WRITE,TEXT):
    fprintf(fd,"/* file generated by atan_ext.mpl*/\n\n"):
    fprintf(fd,"#include \"double_ext.h\"\n"):

    fprintf(fd, "static const double  HALFPI = %1.50e;\n", nearest(Pi/2)):

    fprintf(fd, "#define MIN_REDUCTION_NEEDED %s\n",  printDoubleAsHexInt(e)):
    fprintf(fd,"#define A 0\n"):
    fprintf(fd,"#define B 1\n"):
    fprintf(fd,"#define ATAN_BHI 0\n"):
    fprintf(fd,"#define ATAN_BLO 1\n"):

    fprintf(fd,"#define epsilon %1.50e\n", EpsilonFinal):
    fprintf(fd,"#define epsilon_no_red %1.50e\n",EpsilonFinalNoRed):
    fprintf(fd,"#define TWO_M_64 %1.50e\n",2^(-64)):
    fprintf(fd,"#define TWO_10 %1.50e\n",2^10):

    if 1+1=3 then
        fprintf(fd, "__declspec(align(16)) static const unsigned long long int a_table[%d] = {\n", nbofai );
        for i from 0 to nbofai - 1 do
            fprintf(fd, " /*a[%d] ~= %1.20e   */ ",i,a[i]):
            fprintf(fd, "  %s,\n", printDoubleAsULL(a[i]));
        od:
        fprintf(fd,"\n};\n"):

        fprintf(fd, "static const double b_table[%d] = {\n", nbofai );
        for i from 0 to nbofai - 1 do
            fprintf(fd, "   /*b[%d] = */   %1.50e,\n", i, b[i]):
        od:
        fprintf(fd,"\n};\n"):

    else
        fprintf(fd, "static const struct{long long int a; double b;} ab_table[%d] = {\n", nbofai);
        for i from 0 to nbofai - 1 do
            fprintf(fd, " { /*a[%d] ~= %1.20e   */ ",i,a[i]):
            fprintf(fd, "  %s,\n", printDoubleAsULL(a[i]));
            fprintf(fd, "   /*b[%d] = */   %1.50e},\n", i, b[i]):
        od:
        fprintf(fd,"\n};\n"):

    fi:

    fprintf(fd, "#define atanb_table ((const XC_FLOAT_TYPE *)_atanb_table)\n");
    fprintf(fd, "__declspec(align(16)) static const unsigned short _atanb_table[] = {\n", nbofai );
    for i from 0 to nbofai - 1 do
        fprintf(fd, " /*atan_b[%d] ~= %1.10e*/\n ",i, arctan(b[i])):
        hi:=nearestExt(arctan(b[i]));
        lo:=arctan(b[i])-nearestExt(arctan(b[i]));
        fprintf(fd, " %s, \n", printDoubleExtAsShort(hi));
        fprintf(fd, "   %s, \n", printDoubleExtAsShort(lo));
    od:
    fprintf(fd,"};\n"):

    printPolyExt(fd,Q9,4,"coef_poly"):
    fclose(fd):
end proc;

doth();