1#####################################################################
2# this file contains some useful procedures for the computation of atan in
3# DDE
4
5restart:
6Digits := 200:
7
8with (numapprox):with(orthopoly):
9read "double-extended.mpl";
10read "common-procedures.mpl";
11
12e := 2^(-6.3):
13marge := 2^(-30);
14P19 := convert(series(arctan(x),x=0,20),polynom):
15P18 := (P19-x)/x;
16Q18 := polyExact2Ext(P18, 9);
17Q9 := -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;
18Qprime := poly_exact( -x /3 + x^2/5 - x^3/7 + x^4 /9);
19#log2 (infnorm( (arctan(x)-x*(1+Q18))/x, x=0..e));
20
21#function taken from the "old" coef_atan to compute interval parameters.
22maxx := bi -> simplify( solve( (x-bi) / (1+x*bi) = e ,x) ):
23minx := bi -> simplify(solve( (x-bi) /(1+x*bi)=-e ,x)):
24nextbi := proc (x) evalf( max(solve( minx(bi) = x ,bi) )*(1-marge)); end:
25allbi := proc (n)
26local xi,nbi,x,i,j;
27global b,a, nbofai, nbofbi, valuefordicho;
28    x := e;
29    nbi := 0;
30    i := 0;
31    while(i<n and nbi < 1/e and nbi >= 0) do
32        nbi := nearest ( nextbi(x) );
33        b[i] := evalf( nbi );
34        a[i] := x;
35        x   := evalf(maxx(nbi));
36        i:=i+1;
37    od;
38    j:=0;
39    while ( 2^j < i ) do j:=j+1 od:
40    nbofai := i;
41    nbofbi := i;
42    b[i-1] := nearest(1/e+4):
43    valuefordicho := 2^j;
44    return i,b[0],b[i-1];
45end:
46allbi(100);
47#-------------
48
49#-----------------------
50#quick error calc :
51#-----Reduction:-------
52XredEpsilon := 2^(-64) + 2^(-64) + 2^(-64) + 2^(-64):
53Xred2Epsilon := 2*XredEpsilon + 2^(-64);
54log[2.](Xred2Epsilon);
55errlist:=errlist_quickphase_horner(degree(Qprime),0,0,Xred2Epsilon, 2^(-64));
56errorr := compute_horner_rounding_error(Qprime,x,e,errlist,true);
57qEpsilon := errorr[1]:
58log[2.](qEpsilon);
59deltaApprox := e^11/11:
60polyDelta := e^3*qEpsilon + 2^(-64)*e+ 3*2^(-64)*e^3:
61EpsilonFinal := polyDelta/arctan(e)+2*2^(-64)+deltaApprox/e:
62log[2.](EpsilonFinal);
63#------no reduction-------
64Errlist:=errlist_quickphase_horner(degree(Qprime),0,0,0, 2**(-64)):
65qEpsilon:= compute_horner_rounding_error(Qprime,x,e,errlist,true)[1]:
66deltaApprox := e^11/5:
67EpsilonFinalNoRed := infnorm( (qEpsilon.x^3+2^(-64).x+2^(-64).x^3)/arctan(x),x=2^(-27)..2^(-6) ) + deltaApprox + 2^(-64):
68log[2.](EpsilonFinalNoRed);
69#-----------------------
70
71# Output :
72
73doth:=proc()
74 local filename, fd, i, hi, lo;
75   filename := "TEMPATAN/atan_ext.h":
76    fd := fopen (filename,WRITE,TEXT):
77    fprintf(fd,"/* file generated by atan_ext.mpl*/\n\n"):
78    fprintf(fd,"#include \"double_ext.h\"\n"):
79
80    fprintf(fd, "static const double  HALFPI = %1.50e;\n", nearest(Pi/2)):
81
82    fprintf(fd, "#define MIN_REDUCTION_NEEDED %s\n",  printDoubleAsHexInt(e)):
83    fprintf(fd,"#define A 0\n"):
84    fprintf(fd,"#define B 1\n"):
85    fprintf(fd,"#define ATAN_BHI 0\n"):
86    fprintf(fd,"#define ATAN_BLO 1\n"):
87
88    fprintf(fd,"#define epsilon %1.50e\n", EpsilonFinal):
89    fprintf(fd,"#define epsilon_no_red %1.50e\n",EpsilonFinalNoRed):
90    fprintf(fd,"#define TWO_M_64 %1.50e\n",2^(-64)):
91    fprintf(fd,"#define TWO_10 %1.50e\n",2^10):
92
93    if 1+1=3 then
94        fprintf(fd, "__declspec(align(16)) static const unsigned long long int a_table[%d] = {\n", nbofai );
95        for i from 0 to nbofai - 1 do
96            fprintf(fd, " /*a[%d] ~= %1.20e   */ ",i,a[i]):
97            fprintf(fd, "  %s,\n", printDoubleAsULL(a[i]));
98        od:
99        fprintf(fd,"\n};\n"):
100
101        fprintf(fd, "static const double b_table[%d] = {\n", nbofai );
102        for i from 0 to nbofai - 1 do
103            fprintf(fd, "   /*b[%d] = */   %1.50e,\n", i, b[i]):
104        od:
105        fprintf(fd,"\n};\n"):
106
107    else
108        fprintf(fd, "static const struct{long long int a; double b;} ab_table[%d] = {\n", nbofai);
109        for i from 0 to nbofai - 1 do
110            fprintf(fd, " { /*a[%d] ~= %1.20e   */ ",i,a[i]):
111            fprintf(fd, "  %s,\n", printDoubleAsULL(a[i]));
112            fprintf(fd, "   /*b[%d] = */   %1.50e},\n", i, b[i]):
113        od:
114        fprintf(fd,"\n};\n"):
115
116    fi:
117
118    fprintf(fd, "#define atanb_table ((const XC_FLOAT_TYPE *)_atanb_table)\n");
119    fprintf(fd, "__declspec(align(16)) static const unsigned short _atanb_table[] = {\n", nbofai );
120    for i from 0 to nbofai - 1 do
121        fprintf(fd, " /*atan_b[%d] ~= %1.10e*/\n ",i, arctan(b[i])):
122        hi:=nearestExt(arctan(b[i]));
123        lo:=arctan(b[i])-nearestExt(arctan(b[i]));
124        fprintf(fd, " %s, \n", printDoubleExtAsShort(hi));
125        fprintf(fd, "   %s, \n", printDoubleExtAsShort(lo));
126    od:
127    fprintf(fd,"};\n"):
128
129    printPolyExt(fd,Q9,4,"coef_poly"):
130    fclose(fd):
131end proc;
132
133doth();
134
135