1####################################################################### 2# This file is part of the crlibm library, and is distributed under 3# the LGPL. 4# To use: 5# restart; read "log2-td.mpl"; 6Digits := 120: 7 8interface(quiet=true): 9 10read "common-procedures.mpl": 11read "triple-double.mpl": 12mkdir("TEMPLOG"): 13 14 15# We want log2h + log2m + log2l + delta = log(2) such that 16# log2h and log2m have at least 11 trailing zeros 17# in order to have an exact multiplication with E, which is lower than 1024 in 18# magnitude 19# The resting accuracy is enough for both quick and accurate phases. 20 21log2acc := log(2): 22log2h := round(log2acc * 2**(floor(-log[2](abs(log2acc))) + (53 - 11))) / 2**(floor(-log[2](abs(log2acc))) + (53 - 11)): 23log2m := round((log2acc - log2h) * 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11))) / 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11)): 24log2l := log2acc - (log2h + log2m): 25 26 27L := 7: # number of bits used to address the table 28 29MAXINDEX := round(2^L * (sqrt(2)-1)): 30 31for i from 0 to MAXINDEX-1 do 32 center[i] := 1 + i*2^(-L): # center[i] in [1, 2[ 33 t := evalf(1/center[i]): 34 r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23): 35 (logih[i], logim[i], logil[i]) := hi_mi_lo(evalf(-log(r[i]))): 36od: 37for i from MAXINDEX to 2^L do 38 # y has been divided by two, center[i] in [0.5, 1[ 39 center[i]:=(1 + i*2^(-L)) / 2: 40 t := evalf(1/center[i]): 41 r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23): 42 (logih[i], logim[i], logil[i]) := hi_mi_lo(evalf(-log(r[i]))): 43od: 44 45 46 47 48#Computation of ZMax. 49for i from 0 to MAXINDEX-1 do 50 __x := center[i] + 2^(-L-1) : 51 zmax[i] := (__x*r[i]-1) : 52 __x := center[i] - 2^(-L-1) : 53 zmin[i] := (__x*r[i]-1) : 54od: 55for i from MAXINDEX to 2^L do 56 __x := center[i] + 2^(-L-2) : 57 zmax[i] := (__x*r[i]-1) : 58 __x := center[i] - 2^(-L-2) : 59 zmin[i] := (__x*r[i]-1) : 60od: 61 62zmaxmax:=0: 63zminmin:=0: 64for i from 0 to 2^L do 65 tabulated_value := logih[i] + logim[i] + logil[i]: 66 poly_approx_min := evalf(log(1+zmin[i])): 67 poly_approx_max := evalf(log(1+zmax[i])): 68 69 # Test if we have a case where we cancellate a lot 70 # i.e. the polynomial approximation value is greater than half the tabulated value 71 # the tabulated value is not exactly zero and we are of opposite sign 72 73 if ((abs(poly_approx_min) > 0.75*abs(tabulated_value)) and (tabulated_value <> 0.0) and (poly_approx_min * tabulated_value < 0)) then 74 printf("Polynomial approximation is greater in magnitude in zmin[%d] than half the tabluated value\n",i): 75 76 printf("The tabulated value is %1.50e\n",tabulated_value): 77 if (tabulated_value <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(tabulated_value)))) fi: 78 printf("The value of polynomial in zmin[%d] is %1.50e\n",i,poly_approx_min): 79 if (poly_approx_min <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(poly_approx_min)))) fi: 80 81 summe := poly_approx_min + tabulated_value: 82 printf("The exponent of the sum of both is %d\n",floor(log2(abs(summe)))): 83 84 85 fi: 86 87 if ((abs(poly_approx_max) > 0.75*abs(tabulated_value)) and (tabulated_value <> 0.0) and (poly_approx_max * tabulated_value <0)) then 88 printf("Polynomial approximation is greater in magnitude in zmax[%d] than half the tabluated value\n",i): 89 90 printf("The tabulated value is %1.50e\n",tabulated_value): 91 if (tabulated_value <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(tabulated_value)))) fi: 92 printf("The value of polynomial in zmax[%d] is %1.50e\n",i,poly_approx_max): 93 if (poly_approx_max <> 0.0) then printf("i.e. the value has the exponent %d\n",floor(log2(abs(poly_approx_max)))) fi: 94 95 summe := poly_approx_max + tabulated_value: 96 printf("The exponent of the sum of both is %d\n",floor(log2(abs(summe)))): 97 98 99 fi: 100 101 102 if zmax[i] > zmaxmax then zmaxmax := zmax[i]: fi: 103 if zmin[i] < zminmin then zminmin := zmin[i]: fi: 104od: 105printf("zminmin = -2^(%2f) zmaxmax = 2^(%2f)\n", log2(-zminmin), log2(zmaxmax) ) : 106 107PolyDegreeQuick:=7: 108 109printf(" degree of the polynomial used in the quick phase is %d\n",PolyDegreeQuick); 110 111DDNumberQuick:=3: 112 113printf(" number of double doubles used for the coefficients is %d\n",DDNumberQuick); 114 115#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof 116#and constrain the first two coefficients to 1 and -1/2 in order to save up a full multiplication and a rounding error 117polyQuick:= poly_exact2(x*(1+x*(-0.5+x*(numapprox[minimax]((((log(1+x)/x)-1)/x+0.5)/x, 118 x=-zmaxmax..zmaxmax, [PolyDegreeQuick-3,0], 1 , 'deltaApprox')))), DDNumberQuick): 119 120#Try to verify the bound for using double double arithmetic. 121#Idea: compare the maximum absolute value of the polynomial in zmaxmax (the polynomial has its maxima at the limits) 122#with the maximum value of the first term which is calculated in double precision only 123 124p := unapply(polyQuick,x): 125printf(" using only %d double doubles should be fine since the hi z ulp should affect the result starting from bit %f\n", 126 DDNumberQuick,evalf(53 + log2(p(zmaxmax)) - log2(zmaxmax^(DDNumberQuick)),5)): 127 128 129epsilonApproxQuick := numapprox[infnorm]( 1-polyQuick/log(1+x), x=zminmin..zmaxmax): 130printf(" approximation rel error for the quick phase is 2^(%2f)\n", log2(epsilonApproxQuick) ) : 131deltaApproxQuick := numapprox[infnorm]( polyQuick-log(1+x), x=zminmin..zmaxmax): 132printf(" approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxQuick) ) : 133 134 135PolyDegreeAccurate:=14: 136 137printf(" degree of the polynomial used in the accurate phase is %d\n",PolyDegreeAccurate); 138 139DDNumberAccu:=7: 140TDNumberAccu:=3: 141 142printf(" number of triple doubles used for the coefficients is %d\n",TDNumberAccu); 143printf(" number of double doubles used for the coefficients is %d\n",DDNumberAccu); 144 145 146#Keep -zmaxmax..zmaxmax to keep c1=1, which is useful in the proof 147polyAccurate:= poly_exact32(x*(1+x*(-0.5+x*(numapprox[minimax]((((log(1+x)/x)-1)/x+0.5)/x, 148 x=-zmaxmax..zmaxmax, [PolyDegreeAccurate-3,0], 1 , 'deltaApprox')))), 149 TDNumberAccu, DDNumberAccu): 150 151#Try to verify the bound for using double double arithmetic. 152#Idea: compare the maximum absolute value of the polynomial in zmaxmax (the polynomial has its maxima at the limits) 153#with the maximum value of the first term which is calculated in double precision only 154 155pp := unapply(polyAccurate,x): 156printf(" using only %d triple doubles should be fine since the mi z ulp should affect the result starting from bit %f\n", 157 TDNumberAccu,evalf(106 + log2(p(zmaxmax)) - log2(zmaxmax^(TDNumberAccu)),5)): 158printf(" using only %d double doubles should be fine since the hi z ulp should affect the result starting from bit %f\n", 159 DDNumberAccu,evalf(53 + log2(p(zmaxmax)) - log2(zmaxmax^(TDNumberAccu + DDNumberAccu)),5)): 160 161 162epsilonApproxAccurate := numapprox[infnorm]( 1-polyAccurate/log(1+x), x=zminmin..zmaxmax): 163printf(" approximation rel error for the accurate phase is 2^(%2f)\n", log2(epsilonApproxAccurate) ) : 164deltaApproxAccurate := numapprox[infnorm]( polyAccurate-log(1+x), x=zminmin..zmaxmax): 165printf(" approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxAccurate) ) : 166 167 168 169#Compute now the inverse of ln(10) for the final multiplication of ln(x) with this constant for obtaining log10(x) 170#Compute also the relative error of the constant stored as a triple double. 171Log10inv := evalf(1 / log(10)): 172 173(log10invh, log10invm, log10invl) := hi_mi_lo(Log10inv): 174 175Log10invhml := log10invh + log10invm + log10invl: 176 177epsilonLog10invhml := evalf(abs((Log10invhml - Log10inv) / Log10inv)): 178 179printf(" Log10inv = 1 / ln(10) stored as a triple-double is exact with a relative error of 2^(%2f)\n", 180 evalf(log[2](epsilonLog10invhml))): 181 182 183#Compute now the constant needed for the unfiltered directed final rounding of the triple-double result 184#This constant is supposed to be the reciprocal of the critical accuracy of the function. 185#We suppose this critical accuracy to be 2^(-120) because the worst case (for RD) we know of is x = 403ce41d 8fa665fa 186#We can easily spend some guard bits since we want simply filter out cases with a theoretical worst case accuracy 187#of -infty (exact floating point images) and as we know that we are exact to at least 122 bits. 188 189wca := 2^(122): 190 191 192#------------------------------------------------------------------- 193# Output 194 195 196filename:="TEMPLOG/log10-td.h": 197fd:=fopen(filename, WRITE, TEXT): 198 199fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"): 200 201fprintf(fd, "\n/*File generated by maple/log10-td.mpl*/\n"): 202 203fprintf(fd, "\n\#define L %d\n\n",L): 204fprintf(fd, "\#define MAXINDEX %d\n\n", MAXINDEX): 205fprintf(fd, "\#define INDEXMASK %d\n", 2^L-1): 206fprintf(fd, "\#define two52 %1.50e\n", 2^(52)): 207fprintf(fd, "\#define log2h %1.50e\n", log2h): 208fprintf(fd, "\#define log2m %1.50e\n", log2m): 209fprintf(fd, "\#define log2l %1.50e\n", log2l): 210fprintf(fd, "\#define log10invh %1.50e\n",log10invh): 211fprintf(fd, "\#define log10invm %1.50e\n",log10invm): 212fprintf(fd, "\#define log10invl %1.50e\n",log10invl): 213fprintf(fd, "\#define WORSTCASEACCURACY %1.50e\n",wca): 214 215epsilon_quick_1 := 2^(-61): # The Gappa proof will show this bound 216epsilon_quick_2 := 2^(-61): # The Gappa proof will show this bound 217fprintf(fd, "\#define ROUNDCST1 %1.50e\n", compute_rn_constant(epsilon_quick_1)): 218fprintf(fd, "\#define ROUNDCST2 %1.50e\n", compute_rn_constant(epsilon_quick_2)): 219fprintf(fd, "\#define RDROUNDCST1 %1.50e\n", epsilon_quick_1): 220fprintf(fd, "\#define RDROUNDCST2 %1.50e\n", epsilon_quick_2): 221 222 223fprintf(fd, "\n\n"): 224 225 226# Print the defines for the define statements 227 228for i from 3 to PolyDegreeQuick do 229 fprintf(fd, "\#define c%d %1.50e\n",i,coeff(polyQuick,x,i)): 230od: 231 232fprintf(fd, "\n\n"): 233 234for i from 3 to (DDNumberAccu + TDNumberAccu -1) do 235 (hi,lo) := hi_lo(coeff(polyAccurate,x,i)): 236 fprintf(fd, "\#define accPolyC%dh %1.50e\n",i,hi): 237 fprintf(fd, "\#define accPolyC%dl %1.50e\n",i,lo): 238od: 239 240for i from (DDNumberAccu + TDNumberAccu) to PolyDegreeAccurate do 241 fprintf(fd, "\#define accPolyC%d %1.50e\n",i,coeff(polyAccurate,x,i)): 242od: 243 244fprintf(fd, "\n\n"): 245 246# Print the table 247fprintf(fd, "typedef struct rri_tag {float ri; double logih; double logim; double logil;} rri; \n"): 248fprintf(fd, "static const rri argredtable[%d] = {\n", 2^L): 249for i from 0 to 2^L-1 do 250 fprintf(fd, " { \n"): 251 fprintf(fd, " %1.50e, /* r[%d] */ \n", r[i], i): 252 fprintf(fd, " %1.50e, /* logih[%d] */ \n", logih[i], i): 253 fprintf(fd, " %1.50e, /* logim[%d] */ \n", logim[i], i): 254 fprintf(fd, " %1.50e, /* logil[%d] */ \n", logil[i], i): 255 fprintf(fd, " } "): 256 if(i<2^L-1) then fprintf(fd, ", \n"): fi 257od: 258fprintf(fd, "}; \n \n"): 259 260fclose(fd): 261 262for j from 0 to 2^L-1 do 263 filename:=cat("TEMPLOG/log10-td_",j,".sed"): 264 fd:=fopen(filename, WRITE, TEXT): 265 fprintf(fd, " s/_log2h/%1.50e/g\n", log2h): 266 fprintf(fd, " s/_log2m/%1.50e/g\n", log2m): 267 fprintf(fd, " s/_log2l/%1.50e/g\n", log2l): 268 fprintf(fd, " s/_logih/%1.50e/g\n", logih[j]): 269 fprintf(fd, " s/_logim/%1.50e/g\n", logim[j]): 270 fprintf(fd, " s/_logil/%1.50e/g\n", logil[j]): 271 fprintf(fd, " s/_zmin/%1.50e/g\n", zmin[j]): 272 fprintf(fd, " s/_zmax/%1.50e/g\n", zmax[j]): 273 for i from 3 to PolyDegreeQuick do 274 fprintf(fd, " s/_c%d/%1.50e/g\n", i, coeff(polyQuick,x,i)): 275 od: 276 fprintf(fd, " s/_epsilonApproxQuick/%1.50e/g\n", epsilonApproxQuick): 277 fprintf(fd, " s/_epsilonLog10invhml/%1.50e/g\n", epsilonLog10invhml): 278 fclose(fd): 279 od: 280 281for j from 0 to 2^L-1 do 282 filename:=cat("TEMPLOG/log10-td-accurate_",j,".sed"): 283 fd:=fopen(filename, WRITE, TEXT): 284 fprintf(fd, " s/_log2h/%1.50e/g\n", log2h): 285 fprintf(fd, " s/_log2m/%1.50e/g\n", log2m): 286 fprintf(fd, " s/_log2l/%1.50e/g\n", log2l): 287 fprintf(fd, " s/_logih/%1.50e/g\n", logih[j]): 288 fprintf(fd, " s/_logim/%1.50e/g\n", logim[j]): 289 fprintf(fd, " s/_logil/%1.50e/g\n", logil[j]): 290 fprintf(fd, " s/_zmin/%1.50e/g\n", zmin[j]): 291 fprintf(fd, " s/_zmax/%1.50e/g\n", zmax[j]): 292 for i from 3 to (DDNumberAccu + TDNumberAccu -1) do 293 (hi,lo) := hi_lo(coeff(polyAccurate,x,i)): 294 fprintf(fd, " s/_accPolyC%dh/%1.50e/g\n", i, hi): 295 fprintf(fd, " s/_accPolyC%dl/%1.50e/g\n", i, lo): 296 od: 297 for i from (DDNumberAccu + TDNumberAccu) to PolyDegreeAccurate do 298 fprintf(fd, " s/_accPolyC%d/%1.50e/g\n", i, coeff(polyAccurate,x,i)): 299 od: 300 fprintf(fd, " s/_epsilonApproxAccurate/%1.50e/g\n", epsilonApproxAccurate): 301 fprintf(fd, " s/_epsilonLog10invhml/%1.50e/g\n", epsilonLog10invhml): 302 fclose(fd): 303 od: 304 305# A shell script to use them 306filename:="run-log10-td-proof.sh": 307fd:=fopen(filename, WRITE, TEXT): 308fprintf(fd, "#!/bin/sh\n"): 309fprintf(fd, "# You probably need to edit the path to the gappa executable\n"): 310fprintf(fd, "GAPPA=~/ble/gappa-0.4.5/src/gappa\n"): 311fprintf(fd, "# Test all the possible table value for E=1\n"): 312fprintf(fd, "for num in `seq 0 %d`; do\n", 2^L-1): 313fprintf(fd, " echo $num, E=1:\n"): 314fprintf(fd, " sed -f ./TEMPLOG/log10-td_$num.sed log10-td.gappa | $GAPPA > /dev/null\n"): 315fprintf(fd, " echo\n"): 316fprintf(fd, "done\n"): 317fprintf(fd, "# For the case E=0 we first handle the cases 0 and %d using log10-td-E0-logir0.gappa\n", 2^L): 318fprintf(fd, "echo 0 and %d, E=0:\n", 2^L): 319fprintf(fd, "sed -f log10-td_0.sed log10-td-E0-logir0.gappa | $GAPPA > /dev/null\n"): 320fprintf(fd, "# then the other cases where logirh <> 0\n"): 321fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1): 322fprintf(fd, " echo $num, E=0:\n"): 323fprintf(fd, " sed -f ./TEMPLOG/log10-td_$num.sed log10-td-E0.gappa | $GAPPA > /dev/null\n"): 324fprintf(fd, " echo\n"): 325fprintf(fd, "done\n"): 326fprintf(fd, "# Accurate phase: Test all the possible table value for E=1\n"): 327fprintf(fd, "for num in `seq 0 %d`; do\n", 2^L-1): 328fprintf(fd, " echo Accurate phase: $num, E=1:\n"): 329fprintf(fd, " sed -f ./TEMPLOG/log10-td-accurate_$num.sed log10-td-accurate.gappa | $GAPPA > /dev/null\n"): 330fprintf(fd, " echo\n"): 331fprintf(fd, "done\n"): 332fprintf(fd, "# Accurate phase: For the case E=0 we first handle the cases 0 and %d using log10-td-accurate-E0-logir0.gappa\n", 2^L): 333fprintf(fd, "echo 0 and %d, E=0:\n", 2^L): 334fprintf(fd, "sed -f ./TEMPLOG/log10-td-accurate_0.sed log10-td-accurate-E0-logir0.gappa | $GAPPA > /dev/null\n"): 335fprintf(fd, "# Accurate phase: then the other cases where logirh <> 0\n"): 336fprintf(fd, "for num in `seq 1 %d`; do\n", 2^L-1): 337fprintf(fd, " echo $num, E=0:\n"): 338fprintf(fd, " sed -f ./TEMPLOG/log10-td-accurate_$num.sed log10-td-accurate-E0.gappa | $GAPPA > /dev/null\n"): 339fprintf(fd, " echo\n"): 340fprintf(fd, "done\n"): 341fclose(fd): 342 343printf("----DONE---\n") : 344 345 346