81 lines · plain
1// Implement WorstCase functions to compute the worst case for x mod C, with2// the exponent of x ranges from emin to emax, and precision of x is p.3// Adapted to Sollya from the Maple function in4// J-M. Muller, "Elementary Functions", 3rd ed, Section 11.3.2.5//6// Some examples:7//8// 1) Worst case for trig range reduction fast passes:9//10// Single precision11// > WorstCase(24, -6, 32, pi/32, 128);12// numbermin : 1074188713// expmin : 714// Worst case: 0x1.47d0fep3015// numberofdigits : -32.88816//17// Double precision18// > WorstCase(53, -8, 32, pi/128, 256);19// numbermin : 641102796277577420// expmin : -5321// Worst case: 0x1.6c6cbc45dc8dep-122// numberofdigits : -66.486723//24// 2) Worst case for exponential function range reduction:25//26// Single precision27// > WorstCase(24, 1, 8, log(2), 128);28// numbermin : 290727029// expmin : -2230// Worst case: 0x1.62e43p-131// numberofdigits : -28.967832//33// Double precision34// > WorstCase(53, 0, 11, log(2), 128);35// numbermin : 780414346020669936// expmin : -5137// Worst case: 0x1.bb9d3beb8c86bp138// numberofdigits : -57.493139//40verbosity=0;41procedure WorstCase(p, emin, emax, C, ndigits) {42 epsilonmin = 12345.0;43 Digits = ndigits;44 45 powerofBoverC = 2^(emin - p) / C;46 for e from emin - p + 1 to emax - p + 1 do {47 powerofBoverC = 2 * powerofBoverC;48 a = floor(powerofBoverC);49 Plast = a;50 r = round(1/(powerofBoverC - a), ndigits, RN);51 a = floor(r);52 Qlast = 1;53 Q = a;54 P = Plast * a + 1;55 while (Q < 2^p - 1) do {56 r = round(1/(r - a), ndigits, RN);57 a = floor(r);58 NewQ = Q * a + Qlast;59 NewP = P * a + Plast;60 Qlast = Q;61 Plast = P;62 Q = NewQ;63 P = NewP;64 };65 epsilon = C * abs(Plast - Qlast * powerofBoverC);66 if (epsilon < epsilonmin) then {67 epsilonmin = epsilon;68 numbermin = Qlast;69 expmin = e;70 };71 };72 display=decimal!;73 print("numbermin : ", numbermin);74 print("expmin : ", expmin);75 display=hexadecimal!;76 print("Worst case: ", numbermin * 2^expmin);77 display=decimal!;78 ell = round(log2(epsilonmin), ndigits, RN);79 print("numberofdigits : ", ell);80};81