1####################################################################### 2# This file is part of crlibm, the correctly rounded mathematical library, 3# which has been developed by the Arénaire project at École normale supérieure 4# de Lyon. 5# 6# Copyright (C) 2004-2011 David Defour, Catherine Daramy-Loirat, 7# Florent de Dinechin, Matthieu Gallet, Nicolas Gast, Christoph Quirin Lauter, 8# and Jean-Michel Muller 9# 10# This library is free software; you can redistribute it and/or 11# modify it under the terms of the GNU Lesser General Public 12# License as published by the Free Software Foundation; either 13# version 2.1 of the License, or (at your option) any later version. 14# 15# This library is distributed in the hope that it will be useful, 16# but WITHOUT ANY WARRANTY; without even the implied warranty of 17# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18# Lesser General Public License for more details. 19# 20# You should have received a copy of the GNU Lesser General Public 21# License along with this library; if not, write to the Free Software 22# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 23 24# To use: 25# restart; read "exp-td.mpl"; 26Digits := 145: 27 28interface(quiet=true): 29 30read "common-procedures.mpl": 31read "triple-double.mpl": 32read "gal.mpl": 33mkdir("TEMPACOS"): 34 35with(orthopoly,T): 36 37 38truncPoly := proc(p,k) 39local i, q: 40q := 0: 41convert(q,polynom): 42for i from 0 to min(degree(p,x),k) do 43 q := q + coeff(p,x,i) * x^i: 44od: 45return (q): 46end: 47 48intervals := 10: 49 50for i from 1 to intervals do 51 epsAccurate[i] := 0.5: 52 epsQuick[i] := 0.5: 53 polyAccurate[i] := 1: 54 polyQuick[i] := 1: 55od: 56epsAccurateSpecial := 0.5: 57epsQuickExtra := 0.5: 58 59lowestIntervalMax := 0.185: 60highestIntervalMin := 0.78: 61 62polyAccurateTDCoeffsLowest := 13: 63polyAccurateDDCoeffsLowest := 12: 64polyAccurateDCoeffsLowest := 16: 65 66polyAccurateDegreeLowest := polyAccurateTDCoeffsLowest + polyAccurateDDCoeffsLowest + polyAccurateDCoeffsLowest: 67 68polyQuickDDCoeffsLowest := 12: 69polyQuickDCoeffsLowest := 10: 70 71polyQuickDegreeLowest := polyQuickDDCoeffsLowest + polyQuickDCoeffsLowest: 72 73polyAccurateTDCoeffsMiddle := 7: 74polyAccurateDDCoeffsMiddle := 9: 75polyAccurateDCoeffsMiddle := 19: 76 77polyAccurateDegreeMiddle := polyAccurateTDCoeffsMiddle + polyAccurateDDCoeffsMiddle + polyAccurateDCoeffsMiddle: 78 79polyQuickDDCoeffsMiddle := 7: 80polyQuickDCoeffsMiddle := 7: 81 82polyQuickDegreeMiddle := polyQuickDDCoeffsMiddle + polyQuickDCoeffsMiddle: 83 84 85polyAccurateTDCoeffsHighest := 9: 86polyAccurateDDCoeffsHighest := 9: 87polyAccurateDCoeffsHighest := 11: 88 89polyAccurateDegreeHighest := polyAccurateTDCoeffsHighest + polyAccurateDDCoeffsHighest + polyAccurateDCoeffsHighest: 90 91polyQuickDDCoeffsHighest := 9: 92polyQuickDCoeffsHighest := 9: 93 94polyQuickDegreeHighest := polyQuickDDCoeffsHighest + polyQuickDCoeffsHighest: 95 96extrabound := hexa2ieee(["3F500000","00000000"]): 97polyQuickDegreeExtra := 5: 98 99bound[0] := 0: 100b := nearest(lowestIntervalMax): 101he := ieeehexa(b): 102bound[1] := hexa2ieee([he[1],"00000000"]): 103b := nearest(highestIntervalMin): 104he := ieeehexa(b): 105bound[intervals - 1] := hexa2ieee([he[1],"00000000"]): 106bound[intervals] := 1: 107linearWidth := (highestIntervalMin - lowestIntervalMax) / (intervals - 2): 108scaleSum := 0: 109scale[2] := 1.5: 110scale[3] := 1.35: 111scale[4] := 1.18: 112scale[5] := 1.00: 113scale[6] := 0.915: 114scale[7] := 0.74: 115scale[8] := 0.595: 116scale[9] := 0.50: 117for i from 2 to (intervals-1) do 118 scaleSum := scaleSum + scale[i]: 119od: 120for i from 2 to (intervals-1) do 121 scale[i] := scale[i] / scaleSum; 122od: 123current := lowestIntervalMax: 124for i from 2 to (intervals-2) do 125 b := nearest(current + (highestIntervalMin - lowestIntervalMax) * scale[i]): 126 current := b: 127 he := ieeehexa(b): 128 bound[i] := hexa2ieee([he[1],"00000000"]): 129od: 130printf("Using %d intervals with bounds:\n",intervals): 131for i from 0 to intervals do 132 printf("bound[%d] = %1.30e\n",i,bound[i]): 133od: 134printf("Using an extra bound for truncating the quick poly to deg. %d in interval #0 for small args up to %f (2^(%f))\n", 135 polyQuickDegreeExtra,extrabound,log[2](abs(extrabound))); 136 137printf("Computing Gal's accurate table values for interval midpoints\n"): 138for i from 2 to (intervals-1) do 139 m := nearest((bound[i] + bound[i-1]) / 2): 140 mhe := ieeehexa(m): 141 printf("Interval %d: accurate table research start value %f (%s%s)\n",i,m,mhe[1],mhe[2]): 142 midpointFloat[i] := galDoubleToDoubleDouble(nearest(m),arcsin,2^(-121),2^(15)): 143od: 144 145 146printf("Using the following floating point midpoints for intervals 2-%d:\n",intervals-2): 147for i from 2 to (intervals-1) do 148 mhe := ieeehexa(midpointFloat[i]): 149 printf("midpointFloat[%d] = %f (%s%s)\n",i,midpointFloat[i],mhe[1],mhe[2]): 150od: 151 152printf("The reduced argument z is therefore bounded by:\n"): 153printf("Interval 1: |z| < 2^(%f)\n", 154 log[2](abs(bound[1]))): 155for i from 2 to (intervals-1) do 156 printf("Interval %d: |z| < 2^(%f)\n",i, 157 log[2](max(abs(midpointFloat[i] - bound[i-1]),abs(midpointFloat[i] - bound[i])))): 158od: 159printf("Interval %d: |z| < 2^(%f)\n",intervals, 160 log[2](abs(1 - bound[intervals-1]))): 161 162 163 164printf("Using a %d degree polynomial for lowest interval (1) (accurate phase)\n", 165 polyAccurateDegreeLowest): 166printf("with %d triple-double, %d double-double and %d double coefficients\n", 167 polyAccurateTDCoeffsLowest, polyAccurateDDCoeffsLowest, polyAccurateDCoeffsLowest): 168printf("Using a %d degree polynomial for lowest interval (1) (quick phase)\n", 169 polyQuickDegreeLowest): 170printf("with %d double-double and %d double coefficients\n", 171 polyQuickDDCoeffsLowest, polyQuickDCoeffsLowest): 172printf("Using a %d degree polynomial for middle intervals (2-%d) (accurate phase)\n", 173 polyAccurateDegreeMiddle,intervals-2): 174printf("with %d triple-double, %d double-double and %d double coefficients\n", 175 polyAccurateTDCoeffsMiddle, polyAccurateDDCoeffsMiddle, polyAccurateDCoeffsMiddle): 176printf("Using a %d degree polynomial for middle intervals (2-%d) (quick phase)\n", 177 polyQuickDegreeMiddle,intervals-2): 178printf("with %d double-double and %d double coefficients\n", 179 polyQuickDDCoeffsMiddle, polyQuickDCoeffsMiddle): 180printf("Using a %d degree polynomial for highest interval (%d) (accurate phase)\n", 181 polyAccurateDegreeHighest,intervals): 182printf("with %d triple-double, %d double-double and %d double coefficients\n", 183 polyAccurateTDCoeffsHighest, polyAccurateDDCoeffsHighest, polyAccurateDCoeffsHighest): 184printf("Using a %d degree polynomial for highest interval (%d) (quick phase)\n", 185 polyQuickDegreeHighest,intervals): 186printf("with %d double-double and %d double coefficients\n", 187 polyQuickDDCoeffsHighest, polyQuickDCoeffsHighest): 188 189if true then 190 191printf("Computing polynomials for interval 1 ([%f;%f])\n",bound[0],bound[1]): 192 193f := unapply(convert(series((arcsin(sqrt(x))/sqrt(x) - 1)/x,x=0,300),polynom),x): 194 195p := unapply(eval(numapprox[chebyshev](f(x),x=(bound[0])^2..(bound[1])^2,2^(-127))),x): 196 197polyAccuExact[1] := truncPoly(subs(X=x^2,p(X))*x^3+x,polyAccurateDegreeLowest): 198 199polyAccurate[1] := poly_exact32(polyAccuExact[1],polyAccurateTDCoeffsLowest,polyAccurateDDCoeffsLowest): 200 201epsAccurate[1] := numapprox[infnorm]((polyAccurate[1]/arcsin(x))-1, x=bound[0]..bound[1]): 202epsAccurateSpecial := numapprox[infnorm]((polyAccurate[1]/arcsin(x))-1, x=bound[0]..evalf(sin(2^(-18)))): 203 204epsAccurateUnknown := numapprox[infnorm]((polyAccurate[1]/arcsin(x))-1, x=bound[0]..evalf(cos(12867/8192))): 205 206polyQuickExact[1] := truncPoly(polyAccuExact[1],polyQuickDegreeLowest): 207polyQuick[1] := poly_exact2(polyQuickExact[1],polyQuickDDCoeffsLowest): 208 209epsQuick[1] := numapprox[infnorm]((polyQuick[1]/arcsin(x))-1, x=bound[0]..bound[1]): 210 211polyQuickExtraExact := truncPoly(polyAccuExact[1],polyQuickDegreeExtra): 212polyQuickExtra := poly_exact2(polyQuickExtraExact,polyQuickDegreeExtra): 213epsQuickExtra := numapprox[infnorm]((polyQuickExtra/arcsin(x))-1, x=bound[0]..extrabound): 214 215end if: 216 217for i from 2 to (intervals-1) do 218 219printf("Computing polynomials for interval %d ([%f;%f])\n",i,bound[i-1],bound[i]): 220printf("Reduced argument z will be in interval [%1.8e;%1.8e] ([-2^%f,2^%f]),\n", 221 bound[i-1]-midpointFloat[i],bound[i]-midpointFloat[i], 222 log[2](abs(bound[i-1]-midpointFloat[i])),log[2](abs(bound[i]-midpointFloat[i]))): 223 224f := unapply(arcsin(x+midpointFloat[i]),x): 225 226fhelp := unapply(convert(series((f(x)-arcsin(midpointFloat[i]))/x,x=0,polyAccurateDegreeMiddle*3),polynom),x): 227 228polyAccuExact[i] := numapprox[minimax](fhelp(x), 229 x=bound[i-1]-midpointFloat[i]..bound[i]-midpointFloat[i], 230 [polyAccurateDegreeMiddle,0],1,'err')*x + nearestDD(evalf(arcsin(midpointFloat[i]))): 231 232polyAccurate[i] := poly_exact32(polyAccuExact[i],polyAccurateTDCoeffsMiddle,polyAccurateDDCoeffsMiddle): 233 234epsAccurate[i] := numapprox[infnorm]((polyAccurate[i]/f(x))-1, 235 x=bound[i-1]-midpointFloat[i]..bound[i]-midpointFloat[i]): 236 237polyQuickExact[i] := truncPoly(polyAccuExact[i],polyQuickDegreeMiddle): 238polyQuick[i] := poly_exact2(polyQuickExact[i],polyQuickDDCoeffsMiddle): 239 240epsQuick[i] := numapprox[infnorm]((polyQuick[i]/f(x))-1, 241 x=bound[i-1]-midpointFloat[i]..bound[i]-midpointFloat[i]): 242 243 244od: 245 246 247if true then 248 249printf("Computing polynomials for interval %d ([%f;%f])\n",intervals,bound[intervals-1],bound[intervals]): 250 251 252g := unapply(((arcsin(1 - x) - Pi/2)/sqrt(2*x)),x): 253f := unapply(convert(series(((g(x)+1)/x),x=0,polyAccurateDegreeHighest*4),polynom),x): 254 255 256polyAccuExact[intervals] := numapprox[minimax](f(x),x=(1-bound[intervals]+2^(-53))..(1-bound[intervals-1]), 257 [polyAccurateDegreeHighest-1,0],1,'err')*x-1: 258 259polyAccurate[intervals] := poly_exact32(polyAccuExact[intervals],polyAccurateTDCoeffsHighest,polyAccurateDDCoeffsHighest): 260 261 262epsAccurate[intervals] := numapprox[infnorm](((unapply(polyAccurate[intervals],x)(x)/g(x))-1), 263 x=(1-bound[intervals]+2^(-53))..(1-bound[intervals-1])): 264 265 266polyQuickExact[intervals] := truncPoly(polyAccuExact[intervals],polyQuickDegreeHighest): 267polyQuick[intervals] := poly_exact2(polyQuickExact[intervals],polyQuickDDCoeffsHighest): 268 269epsQuick[intervals] := numapprox[infnorm](((unapply(polyQuick[intervals],x)(x)/g(x))-1), 270 x=(1-bound[intervals]+2^(-53))..(1-bound[intervals-1])): 271 272printf("Checking if the polynomial for interval %d is exactly -1 in z = %f...\n",intervals,1-bound[intervals]); 273 274if (unapply(polyAccurate[intervals],x)(1-bound[intervals]) = -1) then 275 printf(" Check passed!\n"): 276else 277 printf(" Check failed!\n"): 278end if: 279 280end if: 281 282 283 284for i from 1 to intervals do 285 printf("Relative error for accurate phase polynomial in interval %d ([%f;%f]) is 2^(%f)\n", 286 i,bound[i-1],bound[i],log[2](abs(epsAccurate[i]))): 287 printf("Relative error for quick phase polynomial in interval %d ([%f;%f]) is 2^(%f)\n", 288 i,bound[i-1],bound[i],log[2](abs(epsQuick[i]))): 289od: 290printf("Relative error for accurate phase polynomial #1 in special interval [0;sin(2^(-18))]) is 2^(%f)\n", 291 log[2](abs(epsAccurateSpecial))): 292printf("Relative error for accurate phase polynomial #1 in \"unknown\" interval [0;cos(12867/8192)]) is 2^(%f)\n", 293 log[2](abs(epsAccurateUnknown))): 294 295 296printf("Relative error for quick phase extra case truncated polynomial #1 in special interval [0;%f]) is 2^(%f)\n", 297 extrabound,log[2](abs(epsQuickExtra))): 298 299 300(PiHalfH, PiHalfM, PiHalfL) := hi_mi_lo(evalf(Pi/2)): 301(PiH, PiM, PiL) := hi_mi_lo(evalf(Pi)): 302 303epsPiDD := evalf(((PiHalfH + PiHalfM) - Pi/2)/(Pi/2)): 304epsPiTD := evalf(((PiHalfH + PiHalfM + PiHalfL) - Pi/2) / (Pi/2)): 305 306printf("Relative error for storing Pi/2 as a double-double is 2^(%f)\n",evalf(log[2](abs(epsPiDD)))): 307printf("Relative error for storing Pi/2 as a triple-double is 2^(%f)\n",evalf(log[2](abs(epsPiTD)))): 308 309# Ce qui suit est pifometrique et doit etre prouve en Gappa ensuite 310 311arithmeticalErrorQuick[1] := 2^(-61): 312arithmeticalErrorQuick[2] := 2^(-61): 313arithmeticalErrorQuick[3] := 2^(-61): 314arithmeticalErrorQuick[4] := 2^(-61): 315arithmeticalErrorQuick[5] := 2^(-61): 316arithmeticalErrorQuick[6] := 2^(-61): 317arithmeticalErrorQuick[7] := 2^(-61): 318arithmeticalErrorQuick[8] := 2^(-61): 319arithmeticalErrorQuick[9] := 2^(-61): 320arithmeticalErrorQuick[10] := 2^(-68): 321 322arithmeticalErrorQuickExtra := 2^(-80): 323 324for i from 1 to intervals do 325 estimatedOverallEpsQuick[i] := abs(epsQuick[i]) + 326 abs(arithmeticalErrorQuick[i]) + 327 abs(epsQuick[i]) * abs(arithmeticalErrorQuick[i]): 328 printf("Relative quick phase overall error bound to show in Gappa for interval %d ([%f;%f]) is 2^(%f)\n", 329 i,bound[i-1],bound[i],log[2](abs(estimatedOverallEpsQuick[i]))): 330od: 331 332estimatedOverallEpsQuickExtra := abs(epsQuickExtra) + 333 abs(arithmeticalErrorQuickExtra) + 334 abs(epsQuickExtra) * abs(arithmeticalErrorQuickExtra): 335printf("Relative quick phase overall error bound to show for extra truncted poly in interval [0;%f]) is 2^(%f)\n", 336 extrabound,log[2](abs(estimatedOverallEpsQuickExtra))): 337 338 339printf("Write tables...\n"): 340 341filename:="TEMPACOS/acos-td.h": 342fd:=fopen(filename, WRITE, TEXT): 343 344fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"): 345 346fprintf(fd, "\n/* File generated by maple/acos-td.mpl */\n\n"): 347 348printf("Write table with bounds\n"): 349 350fprintf(fd, "/* High order words of interval bounds (low order word is 0) */\n"): 351for i from 1 to (intervals-1) do 352 heb := ieeehexa(bound[i]): 353 fprintf(fd, "\#define BOUND%d 0x%s\n",i,heb[1]): 354od: 355 356heb := ieeehexa(extrabound): 357fprintf(fd, "\#define EXTRABOUND 0x%s\n",heb[1]): 358 359printf("Write additional constants\n"): 360 361fprintf(fd, "\n\n/* Pi/2 as a triple-double*/\n"): 362fprintf(fd, "\#define PIHALFH %1.50e\n",PiHalfH): 363fprintf(fd, "\#define PIHALFM %1.50e\n",PiHalfM): 364fprintf(fd, "\#define PIHALFL %1.50e\n",PiHalfL): 365fprintf(fd, "\n\n/* Pi as a triple-double*/\n"): 366fprintf(fd, "\#define PIH %1.50e\n",PiH): 367fprintf(fd, "\#define PIM %1.50e\n",PiM): 368fprintf(fd, "\#define PIL %1.50e\n",PiL): 369fprintf(fd, "\#define PIHALFDOUBLERN %1.50e\n",nearest(evalf(Pi/2))): 370fprintf(fd, "\#define PIHALFDOUBLERU %1.50e\n",roundUp(evalf(Pi/2))): 371fprintf(fd, "\#define PIHALFDOUBLERD %1.50e\n",roundDown(evalf(Pi/2))): 372 373printf("Write table with midpoints and polynomial coefficients\n"): 374k := 0: 375for i from 0 to polyAccurateDegreeLowest do 376 (hi,mi,lo) := hi_mi_lo(coeff(polyAccurate[1],x,i)): 377 if ((abs(hi) = 1.0) and (mi = 0) and (lo = 0)) then 378 printf( 379 "Coefficient %d of interval 1 polynomial is exactly %f and will not be stored in the table\n",i,hi): 380 else 381 g := 0; 382 if (hi <> 0) then 383 if (i <= polyQuickDegreeLowest) then 384 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dh, quickPolyLowC%dh */",hi,k,i,i): 385 else 386 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dh */",hi,k,i): 387 end if: 388 k := k + 1: 389 end if: 390 if (mi <> 0) then 391 if (i <= polyQuickDDCoeffsLowest) then 392 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dm, quickPolyLowC%dl */",mi,k,i,i): 393 else 394 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dm */",mi,k,i): 395 end if: 396 k := k + 1: 397 end if: 398 if (lo <> 0) then 399 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dl */",lo,k,i): 400 k := k + 1: 401 end if: 402 end if: 403od: 404tbl[k] := sprintf("%1.50e, \t/* %d, RN rounding constant quick poly low*/", 405 compute_rn_constant(estimatedOverallEpsQuick[1]),k): 406k := k + 1: 407tbl[k] := sprintf("%1.50e, \t/* %d, RD rounding constant quick poly low*/", 408 estimatedOverallEpsQuick[1],k): 409k := k + 1: 410printf("Table for interval 1 written\n"): 411for l from 2 to (intervals-1) do 412 tbl[k] := sprintf("%1.50e, \t/* %d, midpoint in interval %d*/",midpointFloat[l],k,l): 413 tblidx[l] := k; 414 k := k + 1; 415 for i from 0 to polyAccurateDegreeMiddle do 416 (hi,mi,lo) := hi_mi_lo(coeff(polyAccurate[l],x,i)): 417 if ((abs(hi) = 1.0) and (mi = 0) and (lo = 0)) then 418 printf( 419 "Coefficient %d of interval %d polynomial is exactly %f and will not be stored in the table\n",i,l,hi): 420 else 421 g := 0; 422 if (hi <> 0) then 423 if (i <= polyQuickDegreeMiddle) then 424 tbl[k] := sprintf( 425 "%1.50e, \t/* %d, accPolyMid%dC%dh, quickPolyMid%dC%dh */",hi,k,l,i,l,i): 426 else 427 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyMid%dC%dh */",hi,k,l,i): 428 end if: 429 k := k + 1: 430 end if: 431 if (mi <> 0) then 432 if (i <= polyQuickDDCoeffsMiddle) then 433 tbl[k] := sprintf( 434 "%1.50e, \t/* %d, accPolyMid%dC%dm, quickPolyMid%dC%dl */",mi,k,l,i,l,i): 435 else 436 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyMid%dC%dm */",mi,k,l,i): 437 end if: 438 k := k + 1: 439 end if: 440 if (lo <> 0) then 441 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyMid%dC%dl */",lo,k,l,i): 442 k := k + 1: 443 end if: 444 end if: 445 od: 446 tbl[k] := sprintf("%1.50e, \t/* %d, RN rounding constant quick poly middle %d*/", 447 compute_rn_constant(estimatedOverallEpsQuick[l]),k,l): 448 k := k + 1: 449 tbl[k] := sprintf("%1.50e, \t/* %d, RD rounding constant quick poly middle %d*/", 450 estimatedOverallEpsQuick[l],k,l): 451 k := k + 1: 452 printf("Table for interval %d written\n",l): 453od: 454tblidx[intervals] := k: 455for i from 0 to polyAccurateDegreeHighest do 456 (hi,mi,lo) := hi_mi_lo(coeff(polyAccurate[intervals],x,i)): 457 if ((abs(hi) = 1.0) and (mi = 0) and (lo = 0)) then 458 printf( 459 "Coefficient %d of interval %d polynomial is exactly %f and will not be stored in the table\n",i,intervals,hi): 460 else 461 g := 0; 462 if (hi <> 0) then 463 if (i <= polyQuickDegreeHighest) then 464 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dh, quickPolyHighC%dh */",hi,k,i,i): 465 else 466 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dh */",hi,k,i): 467 end if: 468 k := k + 1: 469 end if: 470 if (mi <> 0) then 471 if (i <= polyQuickDDCoeffsHighest) then 472 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dm, quickPolyHighC%dl */",mi,k,i,i): 473 else 474 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dm */",mi,k,i): 475 end if: 476 k := k + 1: 477 end if: 478 if (lo <> 0) then 479 tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dl */",lo,k,i): 480 k := k + 1: 481 end if: 482 end if: 483od: 484tbl[k] := sprintf("%1.50e, \t/* %d, RN rounding constant quick poly high*/", 485 compute_rn_constant(estimatedOverallEpsQuick[intervals]),k): 486k := k + 1: 487tbl[k] := sprintf("%1.50e, \t/* %d, RD rounding constant quick poly high*/", 488 estimatedOverallEpsQuick[intervals],k): 489k := k + 1: 490printf("Table for interval %d written\n",intervals): 491tbllen := k: 492printf("The whole table has %d entries, so uses %d bytes of memory\n",tbllen,tbllen*8): 493fprintf(fd,"\n\n/* Indices to the following table */\n"): 494for i from 2 to intervals do 495 fprintf(fd,"\#define TBLIDX%d %d\n",i,tblidx[i]): 496od: 497fprintf(fd, "\n\n/* Table with midpoints and polynomial coefficients */\n"): 498fprintf(fd, "static const double tbl[%d] = {\n",tbllen): 499for i from 0 to (tbllen - 1) do 500 fprintf(fd, "%s\n",tbl[i]): 501od: 502fprintf(fd, "};\n\n"): 503 504fclose(fd): 505 506printf("----DONE---\n"): 507 508 509