1
2#---------------------------------------------------------------------
3# ieeedouble converts a number to IEEE double extended format.
4# returns sign (-1 or 1), exponent between -16383 and 16383, mantissa as a fraction between 0.5 and 1.
5
6# TODO : use JMM procedure; check subnormals etc
7ieeedoubleExt:=proc(xx)
8local x, sign, logabsx, exponent, mantissa, infmantissa;
9    x:=evalf(xx):
10    if (x=0) then
11        sign,exponent,mantissa := 0,0,0;
12    else
13        if (x<0) then sign:=-1:
14        else sign:=1:
15        fi:
16        exponent := floor(log2(sign*x));
17        if (exponent>16383) then mantissa:=infinity: exponent:=16383:
18        elif (exponent< -16382) then
19            # denorm
20            exponent := -16383
21        fi:
22        infmantissa := sign*x*2^(63-exponent);
23        if frac(infmantissa) <> 0.5 then mantissa := round(infmantissa)
24        else
25            mantissa := floor(infmantissa);
26            if type(mantissa,odd) then mantissa := mantissa+1 fi;
27        fi;
28        mantissa := mantissa*2^(-63);
29    fi;
30    sign,exponent,mantissa;
31end:
32
33nearestExt := proc(x)
34    local sign, exponent, mantissa;
35
36    sign, exponent, mantissa := ieeedoubleExt(x);
37    sign*mantissa*2^(exponent);
38end:
39
40
41
42ieeehexaExt:= proc(x)
43local resultat, sign, exponent, mantissa, t;
44
45    if(x=0) then resultat:=["0000","0000","0000","0000","0000"];
46    elif(x=-0) then resultat:=["8000","0000","0000","0000","0000"];
47    else
48        sign,exponent,mantissa := ieeedoubleExt(x);
49        t := 2**80 +  (exponent+16383)*2^64 + mantissa*2^63;
50        if (sign=-1) then
51            t := t + 2**79;
52        fi:
53   t := convert(t, hex);
54   t:=convert(t, string):
55
56   resultat:=[substring(t, 2..5),
57              substring(t, 6..9 ),  substring(t, 10..13 ),
58              substring(t, 14..17 ),  substring(t, 18..21 )];
59
60    end if:
61    resultat;
62end proc:
63
64
65
66printDoubleAsShort:=proc(x)
67    local ss;
68    ss:=ieeehexa(x);
69    cat( "DOUBLE_HEX(",
70         substring(ss[1], 1..4),  ", " ,
71         substring(ss[1], 5..8),  ", " ,
72         substring(ss[2], 1..4),  ", " ,
73         substring(ss[2], 5..8)) ;
74end proc:
75
76printDoubleExtAsShort:=proc(x)
77    local ss;
78    ss:=ieeehexaExt(x);
79    cat( "LDOUBLE_HEX(", ss[1],  ", ", ss[2],  ", ",ss[3],  ", ", ss[4],  ", ", ss[5],  ")");
80end proc:
81
82printDoubleAsULL:=proc(x)
83    local ss;
84    ss:=ieeehexa(x);
85    cat( "ULL(", ss[1], ss[2],  ")");
86end proc:
87
88printDoubleAsHexInt:=proc(x)
89    local ss;
90    ss:=ieeehexa(x);
91    cat(ss[1], ss[2]);
92end proc:
93
94
95
96
97#---------------------------------------------------------------------
98# hi_lo takes an arbitrary precision number x and returns two doubles such that:
99# x ~ x_hi + x_lo
100hiloExt:= proc(x)
101local x_hi, x_lo, res:
102    x_hi:= nearestExt(evalf(x)):
103    res:=x-x_hi:
104    if (res = 0) then
105        x_lo:=0:
106    else
107        x_lo:=nearestExt(evalf(res)):
108    end if;
109    x_hi,x_lo;
110end:
111
112
113#---------------------------------------------------------------------
114# Like poly_exact, but the n first coefficients are exactly representable as the sum of two doubles.
115# (to actually get the two doubles, use procedure hi_lo)
116
117polyExact2Ext:=proc(P,n)
118local deg,i, coef, coef_hi, coef_lo, Q:
119    Q:= 0:
120    convert(Q, polynom):
121    deg:=degree(P,x):
122    for i from 0 to deg do
123        coef :=coeff(P,x,i):
124        coef_hi, coef_lo:=hiloExt(coef):
125        Q:= Q + coef_hi*x^i:
126        if(i<n) then
127            Q := Q + coef_lo*x^i:
128        fi:
129    od:
130    return(Q);
131end:
132
133printPolyExt := proc(fd,P,n, name_of_poly)
134local deg,i, coef, coef_hi, coef_lo;
135    convert(Q, polynom):
136    deg:=degree(P,x):
137    fprintf(fd, " static const long double %s[%d][2] = {\n", name_of_poly, deg+1);
138    for i from 0 to deg do
139        coef :=coeff(P,x,i):
140        coef_hi, coef_lo:=hiloExt(coef):
141
142        fprintf(fd,"{ %1.50eL, ",coef_hi);
143
144        if(i<n) then
145            fprintf(fd," %1.50eL},\n",coef_lo);
146        else
147            fprintf(fd,"0},\n");
148        fi:
149    od:
150    fprintf(fd,"}; \n");
151end:
152