1 /*
2 * This function computes log, correctly rounded,
3 * using triple double arithmetics
4 *
5 * Copyright (C) 2004-2011 David Defour, Catherine Daramy-Loirat,
6 * Florent de Dinechin, Matthieu Gallet, Nicolas Gast, Christoph Quirin Lauter,
7 * and Jean-Michel Muller
8 *
9 * Author: Christoph Lauter
10 * christoph.lauter at ens-lyon.fr
11 *
12 * This file is part of crlibm, the correctly rounded mathematical library,
13 * which has been developed by the Arénaire project at École normale supérieure
14 * de Lyon.
15 *
16 * This library is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU Lesser General Public
18 * License as published by the Free Software Foundation; either
19 * version 2.1 of the License, or (at your option) any later version.
20 *
21 * This library is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24 * Lesser General Public License for more details.
25 *
26 * You should have received a copy of the GNU Lesser General Public
27 * License along with this library; if not, write to the Free Software
28 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
29 */
30
31
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include "crlibm.h"
35 #include "crlibm_private.h"
36 #include "triple-double.h"
37 #include "log-td.h"
38 #ifdef BUILD_INTERVAL_FUNCTIONS
39 #include "interval.h"
40 #endif
41
42
43 #define AVOID_FMA 0
44
45
46
log_td_accurate(double * logh,double * logm,double * logl,int E,double ed,int index,double zh,double zl,double logih,double logim)47 void log_td_accurate(double *logh, double *logm, double *logl, int E, double ed, int index, double zh, double zl, double logih, double logim) {
48 double highPoly, t1h, t1l, t2h, t2l, t3h, t3l, t4h, t4l, t5h, t5l, t6h, t6l, t7h, t7l, t8h, t8l, t9h, t9l, t10h, t10l, t11h, t11l;
49 double t12h, t12l, t13h, t13l, t14h, t14l, zSquareh, zSquarem, zSquarel, zCubeh, zCubem, zCubel, higherPolyMultZh, higherPolyMultZm;
50 double higherPolyMultZl, zSquareHalfh, zSquareHalfm, zSquareHalfl, polyWithSquareh, polyWithSquarem, polyWithSquarel;
51 double polyh, polym, polyl, logil, logyh, logym, logyl, loghover, logmover, loglover, log2edhover, log2edmover, log2edlover;
52 double log2edh, log2edm, log2edl;
53
54
55 #if EVAL_PERF
56 crlibm_second_step_taken++;
57 #endif
58
59
60 /* Accurate phase:
61
62 Argument reduction is already done.
63 We must return logh, logm and logl representing the intermediate result in 118 bits precision.
64
65 We use a 14 degree polynomial, computing the first 3 (the first is 0) coefficients in triple double,
66 calculating the next 7 coefficients in double double arithmetics and the last in double.
67
68 We must account for zl starting with the monome of degree 4 (7^3 + 53 - 7 >> 118); so
69 double double calculations won't account for it.
70
71 */
72
73 /* Start of the horner scheme */
74
75 #if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
76 highPoly = FMA(FMA(FMA(FMA(accPolyC14,zh,accPolyC13),zh,accPolyC12),zh,accPolyC11),zh,accPolyC10);
77 #else
78 highPoly = accPolyC10 + zh * (accPolyC11 + zh * (accPolyC12 + zh * (accPolyC13 + zh * accPolyC14)));
79 #endif
80
81 /* We want to write
82
83 accPolyC3 + zh * (accPoly4 + zh * (accPoly5 + zh * (accPoly6 + zh * (accPoly7 + zh * (accPoly8 + zh * (accPoly9 + zh * highPoly))))));
84 ( t14 t13 t12 t11 t10 t9 t8 t7 t6 t5 t4 t3 t2 t1 )
85
86 with all additions and multiplications in double double arithmetics
87 but we will produce intermediate results labelled t1h/t1l thru t14h/t14l
88 */
89
90 Mul12(&t1h, &t1l, zh, highPoly);
91 Add22(&t2h, &t2l, accPolyC9h, accPolyC9l, t1h, t1l);
92 Mul22(&t3h, &t3l, zh, zl, t2h, t2l);
93 Add22(&t4h, &t4l, accPolyC8h, accPolyC8l, t3h, t3l);
94 Mul22(&t5h, &t5l, zh, zl, t4h, t4l);
95 Add22(&t6h, &t6l, accPolyC7h, accPolyC7l, t5h, t5l);
96 Mul22(&t7h, &t7l, zh, zl, t6h, t6l);
97 Add22(&t8h, &t8l, accPolyC6h, accPolyC6l, t7h, t7l);
98 Mul22(&t9h, &t9l, zh, zl, t8h, t8l);
99 Add22(&t10h, &t10l, accPolyC5h, accPolyC5l, t9h, t9l);
100 Mul22(&t11h, &t11l, zh, zl, t10h, t10l);
101 Add22(&t12h, &t12l, accPolyC4h, accPolyC4l, t11h, t11l);
102 Mul22(&t13h, &t13l, zh, zl, t12h, t12l);
103 Add22(&t14h, &t14l, accPolyC3h, accPolyC3l, t13h, t13l);
104
105 /* We must now prepare (zh + zl)^2 and (zh + zl)^3 as triple doubles */
106
107 Mul23(&zSquareh, &zSquarem, &zSquarel, zh, zl, zh, zl);
108 Mul233(&zCubeh, &zCubem, &zCubel, zh, zl, zSquareh, zSquarem, zSquarel);
109
110 /* We can now multiplicate the middle and higher polynomial by z^3 */
111
112 Mul233(&higherPolyMultZh, &higherPolyMultZm, &higherPolyMultZl, t14h, t14l, zCubeh, zCubem, zCubel);
113
114 /* Multiply now z^2 by -1/2 (exact op) and add to middle and higher polynomial */
115
116 zSquareHalfh = zSquareh * -0.5;
117 zSquareHalfm = zSquarem * -0.5;
118 zSquareHalfl = zSquarel * -0.5;
119
120 Add33(&polyWithSquareh, &polyWithSquarem, &polyWithSquarel,
121 zSquareHalfh, zSquareHalfm, zSquareHalfl,
122 higherPolyMultZh, higherPolyMultZm, higherPolyMultZl);
123
124 /* Add now zh and zl to obtain the polynomial evaluation result */
125
126 Add233(&polyh, &polym, &polyl, zh, zl, polyWithSquareh, polyWithSquarem, polyWithSquarel);
127
128 /* Reconstruct now log(y) = log(1 + z) - log(ri) by adding logih, logim, logil
129 logil has not been read to the time, do this first
130 */
131
132 logil = argredtable[index].logil;
133
134 Add33(&logyh, &logym, &logyl, logih, logim, logil, polyh, polym, polyl);
135
136 /* Multiply log2 with E, i.e. log2h, log2m, log2l by ed
137 ed is always less than 2^(12) and log2h and log2m are stored with at least 12 trailing zeros
138 So multiplying naively is correct (up to 134 bits at least)
139
140 The final result is thus obtained by adding log2 * E to log(y)
141 */
142
143 log2edhover = log2h * ed;
144 log2edmover = log2m * ed;
145 log2edlover = log2l * ed;
146
147 /* It may be necessary to renormalize the tabulated value (multiplied by ed) before adding
148 the to the log(y)-result
149
150 If needed, uncomment the following Renormalize3-Statement and comment out the copies
151 following it.
152 */
153
154 /* Renormalize3(&log2edh, &log2edm, &log2edl, log2edhover, log2edmover, log2edlover); */
155
156 log2edh = log2edhover;
157 log2edm = log2edmover;
158 log2edl = log2edlover;
159
160 Add33(&loghover, &logmover, &loglover, log2edh, log2edm, log2edl, logyh, logym, logyl);
161
162 /* Since we can not guarantee in each addition and multiplication procedure that
163 the results are not overlapping, we must renormalize the result before handing
164 it over to the final rounding
165 */
166
167 Renormalize3(logh,logm,logl,loghover,logmover,loglover);
168
169 }
170
171
172
173 /*************************************************************
174 *************************************************************
175 * ROUNDED TO NEAREST *
176 *************************************************************
177 *************************************************************/
log_rn(double x)178 double log_rn(double x){
179 db_number xdb;
180 double y, ed, ri, logih, logim, yrih, yril, th, zh, zl;
181 double polyHorner, zhSquareh, zhSquarel, polyUpper, zhSquareHalfh, zhSquareHalfl;
182 double t1h, t1l, t2h, t2l, ph, pl, log2edh, log2edl, logTabPolyh, logTabPolyl, logh, logm, logl, roundcst;
183 int E, index;
184
185 E=0;
186 xdb.d=x;
187
188 /* Filter cases */
189 if (xdb.i[HI] < 0x00100000){ /* x < 2^(-1022) */
190 if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0){
191 return -1.0/0.0;
192 } /* log(+/-0) = -Inf */
193 if (xdb.i[HI] < 0){
194 return (x-x)/0; /* log(-x) = Nan */
195 }
196 /* Subnormal number */
197 E = -52;
198 xdb.d *= two52; /* make x a normal number */
199 }
200
201 if (xdb.i[HI] >= 0x7ff00000){
202 return x+x; /* Inf or Nan */
203 }
204
205
206 /* Extract exponent and mantissa
207 Do range reduction,
208 yielding to E holding the exponent and
209 y the mantissa between sqrt(2)/2 and sqrt(2)
210 */
211 E += (xdb.i[HI]>>20)-1023; /* extract the exponent */
212 index = (xdb.i[HI] & 0x000fffff);
213 xdb.i[HI] = index | 0x3ff00000; /* do exponent = 0 */
214 index = (index + (1<<(20-L-1))) >> (20-L);
215
216 /* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
217 if (index >= MAXINDEX){ /* corresponds to xdb>sqrt(2)*/
218 xdb.i[HI] -= 0x00100000;
219 E++;
220 }
221 y = xdb.d;
222 index = index & INDEXMASK;
223 /* Cast integer E into double ed for multiplication later */
224 ed = (double) E;
225
226 /*
227 Read tables:
228 Read one float for ri
229 Read the first two doubles for -log(r_i) (out of three)
230
231 Organization of the table:
232
233 one struct entry per index, the struct entry containing
234 r, logih, logim and logil in this order
235 */
236
237
238 ri = argredtable[index].ri;
239 /*
240 Actually we don't need the logarithm entries now
241 Move the following two lines to the eventual reconstruction
242 As long as we don't have any if in the following code, we can overlap
243 memory access with calculations
244 */
245 logih = argredtable[index].logih;
246 logim = argredtable[index].logim;
247
248 /* Do range reduction:
249
250 zh + zl = y * ri - 1.0 correctly
251
252 Correctness is assured by use of Mul12 and Add12
253 even if we don't force ri to have its' LSBs set to zero
254
255 Discard zl for higher monome degrees
256 */
257
258 Mul12(&yrih, &yril, y, ri);
259 th = yrih - 1.0;
260 Add12Cond(zh, zl, th, yril);
261
262 /*
263 Polynomial evaluation
264
265 Use a 7 degree polynomial
266 Evaluate the higher 5 terms in double precision (-7 * 3 = -21) using Horner's scheme
267 Evaluate the lower 3 terms (the last is 0) in double double precision accounting also for zl
268 using an ad hoc method
269
270 */
271
272
273
274 #if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
275 polyHorner = FMA(FMA(FMA(FMA(c7,zh,c6),zh,c5),zh,c4),zh,c3);
276 #else
277 polyHorner = c3 + zh * (c4 + zh * (c5 + zh * (c6 + zh * c7)));
278 #endif
279
280 Mul12(&zhSquareh, &zhSquarel, zh, zh);
281 polyUpper = polyHorner * (zh * zhSquareh);
282 zhSquareHalfh = zhSquareh * -0.5;
283 zhSquareHalfl = zhSquarel * -0.5;
284 Add12(t1h, t1l, polyUpper, -1 * (zh * zl));
285 Add22(&t2h, &t2l, zh, zl, zhSquareHalfh, zhSquareHalfl);
286 Add22(&ph, &pl, t2h, t2l, t1h, t1l);
287
288 /* Reconstruction
289
290 Read logih and logim in the tables (already done)
291
292 Compute log(x) = E * log(2) + log(1+z) - log(ri)
293 i.e. log(x) = ed * (log2h + log2m) + (ph + pl) + (logih + logim) + delta
294
295 Carry out everything in double double precision
296
297 */
298
299 /*
300 We store log2 as log2h + log2m + log2l where log2h and log2m have 12 trailing zeros
301 Multiplication of ed (double E) and log2h is thus correct
302 The overall accuracy of log2h + log2m + log2l is 53 * 3 - 24 = 135 which
303 is enough for the accurate phase
304 The accuracy suffices also for the quick phase: 53 * 2 - 24 = 82
305 Nevertheless the storage with trailing zeros implies an overlap of the tabulated
306 triple double values. We have to take it into account for the accurate phase
307 basic procedures for addition and multiplication
308 The condition on the next Add12 is verified as log2m is smaller than log2h
309 and both are scaled by ed
310 */
311
312 Add12(log2edh, log2edl, log2h * ed, log2m * ed);
313
314 /* Add logih and logim to ph and pl
315
316 We must use conditioned Add22 as logih can move over ph
317 */
318
319 Add22Cond(&logTabPolyh, &logTabPolyl, logih, logim, ph, pl);
320
321 /* Add log2edh + log2edl to logTabPolyh + logTabPolyl */
322
323 Add22Cond(&logh, &logm, log2edh, log2edl, logTabPolyh, logTabPolyl);
324
325 /* Rounding test and eventual return or call to the accurate function */
326
327 if(E==0)
328 roundcst = ROUNDCST1;
329 else
330 roundcst = ROUNDCST2;
331
332
333 if(logh == (logh + (logm * roundcst)))
334 return logh;
335 else
336 {
337
338 #if DEBUG
339 printf("Going for Accurate Phase for x=%1.50e\n",x);
340 #endif
341
342 log_td_accurate(&logh, &logm, &logl, E, ed, index, zh, zl, logih, logim);
343
344 ReturnRoundToNearest3(logh, logm, logl);
345
346 } /* Accurate phase launched */
347 }
348
349
350 /*************************************************************
351 *************************************************************
352 * ROUNDED UPWARDS *
353 *************************************************************
354 *************************************************************/
log_ru(double x)355 double log_ru(double x) {
356 db_number xdb;
357 double y, ed, ri, logih, logim, yrih, yril, th, zh, zl;
358 double polyHorner, zhSquareh, zhSquarel, polyUpper, zhSquareHalfh, zhSquareHalfl;
359 double t1h, t1l, t2h, t2l, ph, pl, log2edh, log2edl, logTabPolyh, logTabPolyl, logh, logm, logl, roundcst;
360 int E, index;
361
362 if (x == 1.0) return 0.0; /* This the only case in which the image under log of a double is a double. */
363
364 E=0;
365 xdb.d=x;
366
367 /* Filter cases */
368 if (xdb.i[HI] < 0x00100000){ /* x < 2^(-1022) */
369 if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0){
370 return -1.0/0.0;
371 } /* log(+/-0) = -Inf */
372 if (xdb.i[HI] < 0){
373 return (x-x)/0; /* log(-x) = Nan */
374 }
375 /* Subnormal number */
376 E = -52;
377 xdb.d *= two52; /* make x a normal number */
378 }
379
380 if (xdb.i[HI] >= 0x7ff00000){
381 return x+x; /* Inf or Nan */
382 }
383
384
385 /* Extract exponent and mantissa
386 Do range reduction,
387 yielding to E holding the exponent and
388 y the mantissa between sqrt(2)/2 and sqrt(2)
389 */
390 E += (xdb.i[HI]>>20)-1023; /* extract the exponent */
391 index = (xdb.i[HI] & 0x000fffff);
392 xdb.i[HI] = index | 0x3ff00000; /* do exponent = 0 */
393 index = (index + (1<<(20-L-1))) >> (20-L);
394
395 /* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
396 if (index >= MAXINDEX){ /* corresponds to xdb>sqrt(2)*/
397 xdb.i[HI] -= 0x00100000;
398 E++;
399 }
400 y = xdb.d;
401 index = index & INDEXMASK;
402 /* Cast integer E into double ed for multiplication later */
403 ed = (double) E;
404
405 /*
406 Read tables:
407 Read one float for ri
408 Read the first two doubles for -log(r_i) (out of three)
409
410 Organization of the table:
411
412 one struct entry per index, the struct entry containing
413 r, logih, logim and logil in this order
414 */
415
416
417 ri = argredtable[index].ri;
418 /*
419 Actually we don't need the logarithm entries now
420 Move the following two lines to the eventual reconstruction
421 As long as we don't have any if in the following code, we can overlap
422 memory access with calculations
423 */
424 logih = argredtable[index].logih;
425 logim = argredtable[index].logim;
426
427 /* Do range reduction:
428
429 zh + zl = y * ri - 1.0 correctly
430
431 Correctness is assured by use of Mul12 and Add12
432 even if we don't force ri to have its' LSBs set to zero
433
434 Discard zl for higher monome degrees
435 */
436
437 Mul12(&yrih, &yril, y, ri);
438 th = yrih - 1.0;
439 Add12Cond(zh, zl, th, yril);
440
441 /*
442 Polynomial evaluation
443
444 Use a 7 degree polynomial
445 Evaluate the higher 5 terms in double precision (-7 * 3 = -21) using Horner's scheme
446 Evaluate the lower 3 terms (the last is 0) in double double precision accounting also for zl
447 using an ad hoc method
448
449 */
450
451
452
453 #if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
454 polyHorner = FMA(FMA(FMA(FMA(c7,zh,c6),zh,c5),zh,c4),zh,c3);
455 #else
456 polyHorner = c3 + zh * (c4 + zh * (c5 + zh * (c6 + zh * c7)));
457 #endif
458
459 Mul12(&zhSquareh, &zhSquarel, zh, zh);
460 polyUpper = polyHorner * (zh * zhSquareh);
461 zhSquareHalfh = zhSquareh * -0.5;
462 zhSquareHalfl = zhSquarel * -0.5;
463 Add12(t1h, t1l, polyUpper, -1 * (zh * zl));
464 Add22(&t2h, &t2l, zh, zl, zhSquareHalfh, zhSquareHalfl);
465 Add22(&ph, &pl, t2h, t2l, t1h, t1l);
466
467 /* Reconstruction
468
469 Read logih and logim in the tables (already done)
470
471 Compute log(x) = E * log(2) + log(1+z) - log(ri)
472 i.e. log(x) = ed * (log2h + log2m) + (ph + pl) + (logih + logim) + delta
473
474 Carry out everything in double double precision
475
476 */
477
478 /*
479 We store log2 as log2h + log2m + log2l where log2h and log2m have 12 trailing zeros
480 Multiplication of ed (double E) and log2h is thus correct
481 The overall accuracy of log2h + log2m + log2l is 53 * 3 - 24 = 135 which
482 is enough for the accurate phase
483 The accuracy suffices also for the quick phase: 53 * 2 - 24 = 82
484 Nevertheless the storage with trailing zeros implies an overlap of the tabulated
485 triple double values. We have to take it into account for the accurate phase
486 basic procedures for addition and multiplication
487 The condition on the next Add12 is verified as log2m is smaller than log2h
488 and both are scaled by ed
489 */
490
491 Add12(log2edh, log2edl, log2h * ed, log2m * ed);
492
493 /* Add logih and logim to ph and pl
494
495 We must use conditioned Add22 as logih can move over ph
496 */
497
498 Add22Cond(&logTabPolyh, &logTabPolyl, logih, logim, ph, pl);
499
500 /* Add log2edh + log2edl to logTabPolyh + logTabPolyl */
501
502 Add22Cond(&logh, &logm, log2edh, log2edl, logTabPolyh, logTabPolyl);
503
504 /* Rounding test and eventual return or call to the accurate function */
505
506 if(E==0)
507 roundcst = RDROUNDCST1;
508 else
509 roundcst = RDROUNDCST2;
510
511 TEST_AND_RETURN_RU(logh, logm, roundcst);
512
513 #if DEBUG
514 printf("Going for Accurate Phase for x=%1.50e\n",x);
515 #endif
516
517 log_td_accurate(&logh, &logm, &logl, E, ed, index, zh, zl, logih, logim);
518
519 ReturnRoundUpwards3(logh, logm, logl);
520
521 }
522
523
524 /*************************************************************
525 *************************************************************
526 * ROUNDED DOWNWARDS *
527 *************************************************************
528 *************************************************************/
log_rd(double x)529 double log_rd(double x) {
530 db_number xdb;
531 double y, ed, ri, logih, logim, yrih, yril, th, zh, zl;
532 double polyHorner, zhSquareh, zhSquarel, polyUpper, zhSquareHalfh, zhSquareHalfl;
533 double t1h, t1l, t2h, t2l, ph, pl, log2edh, log2edl, logTabPolyh, logTabPolyl, logh, logm, logl, roundcst;
534 int E, index;
535
536 if (x == 1.0) return 0.0; /* This the only case in which the image under log of a double is a double. */
537
538 E=0;
539 xdb.d=x;
540
541 /* Filter cases */
542 if (xdb.i[HI] < 0x00100000){ /* x < 2^(-1022) */
543 if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0){
544 return -1.0/0.0;
545 } /* log(+/-0) = -Inf */
546 if (xdb.i[HI] < 0){
547 return (x-x)/0; /* log(-x) = Nan */
548 }
549 /* Subnormal number */
550 E = -52;
551 xdb.d *= two52; /* make x a normal number */
552 }
553
554 if (xdb.i[HI] >= 0x7ff00000){
555 return x+x; /* Inf or Nan */
556 }
557
558
559 /* Extract exponent and mantissa
560 Do range reduction,
561 yielding to E holding the exponent and
562 y the mantissa between sqrt(2)/2 and sqrt(2)
563 */
564 E += (xdb.i[HI]>>20)-1023; /* extract the exponent */
565 index = (xdb.i[HI] & 0x000fffff);
566 xdb.i[HI] = index | 0x3ff00000; /* do exponent = 0 */
567 index = (index + (1<<(20-L-1))) >> (20-L);
568
569 /* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
570 if (index >= MAXINDEX){ /* corresponds to xdb>sqrt(2)*/
571 xdb.i[HI] -= 0x00100000;
572 E++;
573 }
574 y = xdb.d;
575 index = index & INDEXMASK;
576 /* Cast integer E into double ed for multiplication later */
577 ed = (double) E;
578
579 /*
580 Read tables:
581 Read one float for ri
582 Read the first two doubles for -log(r_i) (out of three)
583
584 Organization of the table:
585
586 one struct entry per index, the struct entry containing
587 r, logih, logim and logil in this order
588 */
589
590
591 ri = argredtable[index].ri;
592 /*
593 Actually we don't need the logarithm entries now
594 Move the following two lines to the eventual reconstruction
595 As long as we don't have any if in the following code, we can overlap
596 memory access with calculations
597 */
598 logih = argredtable[index].logih;
599 logim = argredtable[index].logim;
600
601 /* Do range reduction:
602
603 zh + zl = y * ri - 1.0 correctly
604
605 Correctness is assured by use of Mul12 and Add12
606 even if we don't force ri to have its' LSBs set to zero
607
608 Discard zl for higher monome degrees
609 */
610
611 Mul12(&yrih, &yril, y, ri);
612 th = yrih - 1.0;
613 Add12Cond(zh, zl, th, yril);
614
615 /*
616 Polynomial evaluation
617
618 Use a 7 degree polynomial
619 Evaluate the higher 5 terms in double precision (-7 * 3 = -21) using Horner's scheme
620 Evaluate the lower 3 terms (the last is 0) in double double precision accounting also for zl
621 using an ad hoc method
622
623 */
624
625
626
627 #if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
628 polyHorner = FMA(FMA(FMA(FMA(c7,zh,c6),zh,c5),zh,c4),zh,c3);
629 #else
630 polyHorner = c3 + zh * (c4 + zh * (c5 + zh * (c6 + zh * c7)));
631 #endif
632
633 Mul12(&zhSquareh, &zhSquarel, zh, zh);
634 polyUpper = polyHorner * (zh * zhSquareh);
635 zhSquareHalfh = zhSquareh * -0.5;
636 zhSquareHalfl = zhSquarel * -0.5;
637 Add12(t1h, t1l, polyUpper, -1 * (zh * zl));
638 Add22(&t2h, &t2l, zh, zl, zhSquareHalfh, zhSquareHalfl);
639 Add22(&ph, &pl, t2h, t2l, t1h, t1l);
640
641 /* Reconstruction
642
643 Read logih and logim in the tables (already done)
644
645 Compute log(x) = E * log(2) + log(1+z) - log(ri)
646 i.e. log(x) = ed * (log2h + log2m) + (ph + pl) + (logih + logim) + delta
647
648 Carry out everything in double double precision
649
650 */
651
652 /*
653 We store log2 as log2h + log2m + log2l where log2h and log2m have 12 trailing zeros
654 Multiplication of ed (double E) and log2h is thus correct
655 The overall accuracy of log2h + log2m + log2l is 53 * 3 - 24 = 135 which
656 is enough for the accurate phase
657 The accuracy suffices also for the quick phase: 53 * 2 - 24 = 82
658 Nevertheless the storage with trailing zeros implies an overlap of the tabulated
659 triple double values. We have to take it into account for the accurate phase
660 basic procedures for addition and multiplication
661 The condition on the next Add12 is verified as log2m is smaller than log2h
662 and both are scaled by ed
663 */
664
665 Add12(log2edh, log2edl, log2h * ed, log2m * ed);
666
667 /* Add logih and logim to ph and pl
668
669 We must use conditioned Add22 as logih can move over ph
670 */
671
672 Add22Cond(&logTabPolyh, &logTabPolyl, logih, logim, ph, pl);
673
674 /* Add log2edh + log2edl to logTabPolyh + logTabPolyl */
675
676 Add22Cond(&logh, &logm, log2edh, log2edl, logTabPolyh, logTabPolyl);
677
678 /* Rounding test and eventual return or call to the accurate function */
679
680 if(E==0)
681 roundcst = RDROUNDCST1;
682 else
683 roundcst = RDROUNDCST2;
684
685 TEST_AND_RETURN_RD(logh, logm, roundcst);
686
687 #if DEBUG
688 printf("Going for Accurate Phase for x=%1.50e\n",x);
689 #endif
690
691 log_td_accurate(&logh, &logm, &logl, E, ed, index, zh, zl, logih, logim);
692
693 ReturnRoundDownwards3(logh, logm, logl);
694 }
695
696 /*************************************************************
697 *************************************************************
698 * ROUNDED TOWARDS ZERO *
699 *************************************************************
700 *************************************************************/
log_rz(double x)701 double log_rz(double x) {
702 db_number xdb;
703 double y, ed, ri, logih, logim, yrih, yril, th, zh, zl;
704 double polyHorner, zhSquareh, zhSquarel, polyUpper, zhSquareHalfh, zhSquareHalfl;
705 double t1h, t1l, t2h, t2l, ph, pl, log2edh, log2edl, logTabPolyh, logTabPolyl, logh, logm, logl, roundcst;
706 int E, index;
707
708 if (x == 1.0) return 0.0; /* This the only case in which the image under log of a double is a double. */
709
710 E=0;
711 xdb.d=x;
712
713 /* Filter cases */
714 if (xdb.i[HI] < 0x00100000){ /* x < 2^(-1022) */
715 if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0){
716 return -1.0/0.0;
717 } /* log(+/-0) = -Inf */
718 if (xdb.i[HI] < 0){
719 return (x-x)/0; /* log(-x) = Nan */
720 }
721 /* Subnormal number */
722 E = -52;
723 xdb.d *= two52; /* make x a normal number */
724 }
725
726 if (xdb.i[HI] >= 0x7ff00000){
727 return x+x; /* Inf or Nan */
728 }
729
730
731 /* Extract exponent and mantissa
732 Do range reduction,
733 yielding to E holding the exponent and
734 y the mantissa between sqrt(2)/2 and sqrt(2)
735 */
736 E += (xdb.i[HI]>>20)-1023; /* extract the exponent */
737 index = (xdb.i[HI] & 0x000fffff);
738 xdb.i[HI] = index | 0x3ff00000; /* do exponent = 0 */
739 index = (index + (1<<(20-L-1))) >> (20-L);
740
741 /* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
742 if (index >= MAXINDEX){ /* corresponds to xdb>sqrt(2)*/
743 xdb.i[HI] -= 0x00100000;
744 E++;
745 }
746 y = xdb.d;
747 index = index & INDEXMASK;
748 /* Cast integer E into double ed for multiplication later */
749 ed = (double) E;
750
751 /*
752 Read tables:
753 Read one float for ri
754 Read the first two doubles for -log(r_i) (out of three)
755
756 Organization of the table:
757
758 one struct entry per index, the struct entry containing
759 r, logih, logim and logil in this order
760 */
761
762
763 ri = argredtable[index].ri;
764 /*
765 Actually we don't need the logarithm entries now
766 Move the following two lines to the eventual reconstruction
767 As long as we don't have any if in the following code, we can overlap
768 memory access with calculations
769 */
770 logih = argredtable[index].logih;
771 logim = argredtable[index].logim;
772
773 /* Do range reduction:
774
775 zh + zl = y * ri - 1.0 correctly
776
777 Correctness is assured by use of Mul12 and Add12
778 even if we don't force ri to have its' LSBs set to zero
779
780 Discard zl for higher monome degrees
781 */
782
783 Mul12(&yrih, &yril, y, ri);
784 th = yrih - 1.0;
785 Add12Cond(zh, zl, th, yril);
786
787 /*
788 Polynomial evaluation
789
790 Use a 7 degree polynomial
791 Evaluate the higher 5 terms in double precision (-7 * 3 = -21) using Horner's scheme
792 Evaluate the lower 3 terms (the last is 0) in double double precision accounting also for zl
793 using an ad hoc method
794
795 */
796
797
798
799 #if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
800 polyHorner = FMA(FMA(FMA(FMA(c7,zh,c6),zh,c5),zh,c4),zh,c3);
801 #else
802 polyHorner = c3 + zh * (c4 + zh * (c5 + zh * (c6 + zh * c7)));
803 #endif
804
805 Mul12(&zhSquareh, &zhSquarel, zh, zh);
806 polyUpper = polyHorner * (zh * zhSquareh);
807 zhSquareHalfh = zhSquareh * -0.5;
808 zhSquareHalfl = zhSquarel * -0.5;
809 Add12(t1h, t1l, polyUpper, -1 * (zh * zl));
810 Add22(&t2h, &t2l, zh, zl, zhSquareHalfh, zhSquareHalfl);
811 Add22(&ph, &pl, t2h, t2l, t1h, t1l);
812
813 /* Reconstruction
814
815 Read logih and logim in the tables (already done)
816
817 Compute log(x) = E * log(2) + log(1+z) - log(ri)
818 i.e. log(x) = ed * (log2h + log2m) + (ph + pl) + (logih + logim) + delta
819
820 Carry out everything in double double precision
821
822 */
823
824 /*
825 We store log2 as log2h + log2m + log2l where log2h and log2m have 12 trailing zeros
826 Multiplication of ed (double E) and log2h is thus correct
827 The overall accuracy of log2h + log2m + log2l is 53 * 3 - 24 = 135 which
828 is enough for the accurate phase
829 The accuracy suffices also for the quick phase: 53 * 2 - 24 = 82
830 Nevertheless the storage with trailing zeros implies an overlap of the tabulated
831 triple double values. We have to take it into account for the accurate phase
832 basic procedures for addition and multiplication
833 The condition on the next Add12 is verified as log2m is smaller than log2h
834 and both are scaled by ed
835 */
836
837 Add12(log2edh, log2edl, log2h * ed, log2m * ed);
838
839 /* Add logih and logim to ph and pl
840
841 We must use conditioned Add22 as logih can move over ph
842 */
843
844 Add22Cond(&logTabPolyh, &logTabPolyl, logih, logim, ph, pl);
845
846 /* Add log2edh + log2edl to logTabPolyh + logTabPolyl */
847
848 Add22Cond(&logh, &logm, log2edh, log2edl, logTabPolyh, logTabPolyl);
849
850 /* Rounding test and eventual return or call to the accurate function */
851
852 if(E==0)
853 roundcst = RDROUNDCST1;
854 else
855 roundcst = RDROUNDCST2;
856
857 TEST_AND_RETURN_RZ(logh, logm, roundcst);
858
859 #if DEBUG
860 printf("Going for Accurate Phase for x=%1.50e\n",x);
861 #endif
862
863 log_td_accurate(&logh, &logm, &logl, E, ed, index, zh, zl, logih, logim);
864
865 ReturnRoundTowardsZero3(logh, logm, logl);
866 }
867
868 #ifdef BUILD_INTERVAL_FUNCTIONS
j_log(interval x)869 interval j_log(interval x)
870 {
871 interval res;
872 int roundable;
873 int cs_inf=0; int cs_sup=0;
874 double x_inf,x_sup;
875 x_inf=LOW(x);
876 x_sup=UP(x);
877 double res_inf, res_sup, res_simple_inf, res_simple_sup;
878
879 db_number xdb_sup;
880 double y_sup, ed_sup, ri_sup, logih_sup, logim_sup, yrih_sup, yril_sup, th_sup, zh_sup, zl_sup;
881 double polyHorner_sup, zhSquareh_sup, zhSquarel_sup, polyUpper_sup, zhSquareHalfh_sup, zhSquareHalfl_sup;
882 double t1h_sup, t1l_sup, t2h_sup, t2l_sup, ph_sup, pl_sup, log2edh_sup, log2edl_sup, logTabPolyh_sup, logTabPolyl_sup, logh_sup, logm_sup, logl_sup, roundcst;
883 int E_sup, index_sup;
884
885
886 db_number xdb_inf;
887 double y_inf, ed_inf, ri_inf, logih_inf, logim_inf, yrih_inf, yril_inf, th_inf, zh_inf, zl_inf;
888 double polyHorner_inf, zhSquareh_inf, zhSquarel_inf, polyUpper_inf, zhSquareHalfh_inf, zhSquareHalfl_inf;
889 double t1h_inf, t1l_inf, t2h_inf, t2l_inf, ph_inf, pl_inf, log2edh_inf, log2edl_inf, logTabPolyh_inf, logTabPolyl_inf, logh_inf, logm_inf, logl_inf;
890 int E_inf, index_inf;
891
892 E_inf=0;
893 xdb_inf.d=x_inf;
894 E_sup=0;
895 xdb_sup.d=x_sup;
896
897 if (__builtin_expect(
898 (x_inf == 1.0) || (!(x_inf<=x_sup)) || (xdb_sup.i[HI] < 0)
899 || (xdb_inf.i[HI] < 0x00100000) || (((xdb_inf.i[HI] & 0x7fffffff)|xdb_inf.i[LO])==0)
900 || (xdb_inf.i[HI] < 0)
901 || (xdb_inf.i[HI] >= 0x7ff00000)
902 || (x_sup == 1.0)
903 || (xdb_sup.i[HI] < 0x00100000)
904 || (((xdb_sup.i[HI] & 0x7fffffff)|xdb_sup.i[LO])==0)
905 || (xdb_sup.i[HI] < 0)
906 || (xdb_sup.i[HI] >= 0x7ff00000)
907 || ((xdb_inf.d<00) && (xdb_sup.d>0) )
908 ,FALSE))
909 {
910 if (!(x_inf<=x_sup)) RETURN_EMPTY_INTERVAL;
911 if (xdb_sup.i[HI] < 0) RETURN_EMPTY_INTERVAL;
912 if ((xdb_inf.d<00) && (xdb_sup.d>0) )
913 {
914 ASSIGN_LOW(res,-1.0/0.0);
915 ASSIGN_UP(res,log_ru(UP(x)));
916 return res;
917 }
918 ASSIGN_LOW(res,log_rd(LOW(x)));
919 ASSIGN_UP(res,log_ru(UP(x)));
920 return res;
921 }
922 /* Extract exponent and mantissa
923 Do range reduction,
924 yielding to E holding the exponent and
925 y the mantissa between sqrt(2)/2 and sqrt(2)
926 */
927 E_inf += (xdb_inf.i[HI]>>20)-1023; /* extract the exponent */
928 E_sup += (xdb_sup.i[HI]>>20)-1023; /* extract the exponent */
929 index_inf = (xdb_inf.i[HI] & 0x000fffff);
930 index_sup = (xdb_sup.i[HI] & 0x000fffff);
931 xdb_inf.i[HI] = index_inf | 0x3ff00000; /* do exponent = 0 */
932 xdb_sup.i[HI] = index_sup | 0x3ff00000; /* do exponent = 0 */
933 index_inf = (index_inf + (1<<(20-L-1))) >> (20-L);
934 index_sup = (index_sup + (1<<(20-L-1))) >> (20-L);
935 /* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
936 if (index_inf >= MAXINDEX){ /* corresponds to xdb>sqrt(2)*/
937 xdb_inf.i[HI] -= 0x00100000;
938 E_inf++;
939 }
940 /* reduce such that sqrt(2)/2 < xdb.d < sqrt(2) */
941 if (index_sup >= MAXINDEX){ /* corresponds to xdb>sqrt(2)*/
942 xdb_sup.i[HI] -= 0x00100000;
943 E_sup++;
944 }
945 y_inf = xdb_inf.d;
946 y_sup = xdb_sup.d;
947 index_inf = index_inf & INDEXMASK;
948 index_sup = index_sup & INDEXMASK;
949 /* Cast integer E into double ed for multiplication later */
950 ed_inf = (double) E_inf;
951 ed_sup = (double) E_sup;
952 /*
953 Read tables:
954 Read one float for ri
955 Read the first two doubles for -log(r_i) (out of three)
956 Organization of the table:
957 one struct entry per index, the struct entry containing
958 r, logih, logim and logil in this order
959 */
960 ri_inf = argredtable[index_inf].ri;
961 ri_sup = argredtable[index_sup].ri;
962 /*
963 Actually we don't need the logarithm entries now
964 Move the following two lines to the eventual reconstruction
965 As long as we don't have any if in the following code, we can overlap
966 memory access with calculations
967 */
968 logih_inf = argredtable[index_inf].logih;
969 logih_sup = argredtable[index_sup].logih;
970 logim_inf = argredtable[index_inf].logim;
971 logim_sup = argredtable[index_sup].logim;
972 /* Do range reduction:
973 zh + zl = y * ri - 1.0 correctly
974 Correctness is assured by use of Mul12 and Add12
975 even if we don't force ri to have its' LSBs set to zero
976 Discard zl for higher monome degrees
977 */
978
979 Mul12(&yrih_inf, &yril_inf, y_inf, ri_inf);
980 Mul12(&yrih_sup, &yril_sup, y_sup, ri_sup);
981 th_inf = yrih_inf - 1.0;
982 th_sup = yrih_sup - 1.0;
983 /* Do range reduction:
984 zh + zl = y * ri - 1.0 correctly
985 Correctness is assured by use of Mul12 and Add12
986 even if we don't force ri to have its' LSBs set to zero
987 Discard zl for higher monome degrees
988 */
989
990
991 Add12Cond(zh_inf, zl_inf, th_inf, yril_inf);
992 Add12Cond(zh_sup, zl_sup, th_sup, yril_sup);
993 /*
994 Polynomial evaluation
995 Use a 7 degree polynomial
996 Evaluate the higher 5 terms in double precision (-7 * 3 = -21) using Horner's scheme
997 Evaluate the lower 3 terms (the last is 0) in double double precision accounting also for zl
998 using an ad hoc method
999 */
1000
1001
1002
1003 #if defined(PROCESSOR_HAS_FMA) && !defined(AVOID_FMA)
1004 polyHorner_inf = FMA(FMA(FMA(FMA(c7,zh_inf,c6),zh_inf,c5),zh_inf,c4),zh_inf,c3);
1005 polyHorner_sup = FMA(FMA(FMA(FMA(c7,zh_sup,c6),zh_sup,c5),zh_sup,c4),zh_sup,c3);
1006 #else
1007 polyHorner_inf = c3 + zh_inf * (c4 + zh_inf * (c5 + zh_inf * (c6 + zh_inf * c7)));
1008 polyHorner_sup = c3 + zh_sup * (c4 + zh_sup * (c5 + zh_sup * (c6 + zh_sup * c7)));
1009 #endif
1010
1011 Mul12(&zhSquareh_inf, &zhSquarel_inf, zh_inf, zh_inf);
1012 Mul12(&zhSquareh_sup, &zhSquarel_sup, zh_sup, zh_sup);
1013 polyUpper_inf = polyHorner_inf * (zh_inf * zhSquareh_inf);
1014 polyUpper_sup = polyHorner_sup * (zh_sup * zhSquareh_sup);
1015 zhSquareHalfh_inf = zhSquareh_inf * -0.5;
1016 zhSquareHalfh_sup = zhSquareh_sup * -0.5;
1017 zhSquareHalfl_inf = zhSquarel_inf * -0.5;
1018 zhSquareHalfl_sup = zhSquarel_sup * -0.5;
1019 Add12(t1h_inf, t1l_inf, polyUpper_inf, -1 * (zh_inf * zl_inf));
1020 Add12(t1h_sup, t1l_sup, polyUpper_sup, -1 * (zh_sup * zl_sup));
1021 Add22(&t2h_inf, &t2l_inf, zh_inf, zl_inf, zhSquareHalfh_inf, zhSquareHalfl_inf);
1022 Add22(&t2h_sup, &t2l_sup, zh_sup, zl_sup, zhSquareHalfh_sup, zhSquareHalfl_sup);
1023 Add22(&ph_inf, &pl_inf, t2h_inf, t2l_inf, t1h_inf, t1l_inf);
1024 Add22(&ph_sup, &pl_sup, t2h_sup, t2l_sup, t1h_sup, t1l_sup);
1025 /* Reconstruction
1026
1027 Read logih and logim in the tables (already done)
1028
1029 Compute log(x) = E * log(2) + log(1+z) - log(ri)
1030 i.e. log(x) = ed * (log2h + log2m) + (ph + pl) + (logih + logim) + delta
1031
1032 Carry out everything in double double precision
1033
1034 */
1035
1036 /*
1037 We store log2 as log2h + log2m + log2l where log2h and log2m have 12 trailing zeros
1038 Multiplication of ed (double E) and log2h is thus correct
1039 The overall accuracy of log2h + log2m + log2l is 53 * 3 - 24 = 135 which
1040 is enough for the accurate phase
1041 The accuracy suffices also for the quick phase: 53 * 2 - 24 = 82
1042 Nevertheless the storage with trailing zeros implies an overlap of the tabulated
1043 triple double values. We have to take it into account for the accurate phase
1044 basic procedures for addition and multiplication
1045 The condition on the next Add12 is verified as log2m is smaller than log2h
1046 and both are scaled by ed
1047 */
1048
1049 Add12(log2edh_inf, log2edl_inf, log2h * ed_inf, log2m * ed_inf);
1050
1051 /*
1052 We store log2 as log2h + log2m + log2l where log2h and log2m have 12 trailing zeros
1053 Multiplication of ed (double E) and log2h is thus correct
1054 The overall accuracy of log2h + log2m + log2l is 53 * 3 - 24 = 135 which
1055 is enough for the accurate phase
1056 The accuracy suffices also for the quick phase: 53 * 2 - 24 = 82
1057 Nevertheless the storage with trailing zeros implies an overlap of the tabulated
1058 triple double values. We have to take it into account for the accurate phase
1059 basic procedures for addition and multiplication
1060 The condition on the next Add12 is verified as log2m is smaller than log2h
1061 and both are scaled by ed
1062 */
1063
1064 Add12(log2edh_sup, log2edl_sup, log2h * ed_sup, log2m * ed_sup);
1065
1066
1067 /* Add logih and logim to ph and pl
1068
1069 We must use conditioned Add22 as logih can move over ph
1070 */
1071
1072 Add22Cond(&logTabPolyh_inf, &logTabPolyl_inf, logih_inf, logim_inf, ph_inf, pl_inf);
1073
1074 /* Add log2edh + log2edl to logTabPolyh + logTabPolyl */
1075
1076 Add22Cond(&logh_inf, &logm_inf, log2edh_inf, log2edl_inf, logTabPolyh_inf, logTabPolyl_inf);
1077
1078 /* Add logih and logim to ph and pl
1079
1080 We must use conditioned Add22 as logih can move over ph
1081 */
1082
1083 Add22Cond(&logTabPolyh_sup, &logTabPolyl_sup, logih_sup, logim_sup, ph_sup, pl_sup);
1084
1085 /* Add log2edh + log2edl to logTabPolyh + logTabPolyl */
1086
1087 Add22Cond(&logh_sup, &logm_sup, log2edh_sup, log2edl_sup, logTabPolyh_sup, logTabPolyl_sup);
1088
1089
1090 /* Rounding test and eventual return or call to the accurate function */
1091
1092 roundcst = RDROUNDCST1;
1093
1094 if(cs_inf) { res_inf=res_simple_inf; }
1095 if(cs_sup) { res_sup=res_simple_sup; }
1096 //TEST_AND_COPY_RDRU_LOG(roundable,res_inf,logh_inf,logm_inf,res_sup,logh_sup,logm_sup,roundcst);
1097
1098 //#define TEST_AND_COPY_RDRU_LOG(__cond__, __res_inf__, __yh_inf__, __yl_inf__, __res_sup__, __yh_sup__, __yl_sup__, __eps__)
1099
1100 db_number yh_inf, yl_inf, u53_inf, yh_sup, yl_sup, u53_sup;
1101 int yh_inf_neg, yl_inf_neg, yh_sup_neg, yl_sup_neg;
1102 int rd_ok, ru_ok;
1103 double save_res_inf=res_inf;
1104 yh_inf.d = logh_inf; yl_inf.d = logm_inf;
1105 yh_inf_neg = (yh_inf.i[HI] & 0x80000000);
1106 yl_inf_neg = (yl_inf.i[HI] & 0x80000000);
1107 yh_inf.l = yh_inf.l & 0x7fffffffffffffffLL; /* compute the absolute value*/
1108 yl_inf.l = yl_inf.l & 0x7fffffffffffffffLL; /* compute the absolute value*/
1109 u53_inf.l = (yh_inf.l & ULL(7ff0000000000000)) + ULL(0010000000000000);
1110 yh_sup.d = logh_sup; yl_sup.d = logm_sup;
1111 yh_sup_neg = (yh_sup.i[HI] & 0x80000000);
1112 yl_sup_neg = (yl_sup.i[HI] & 0x80000000);
1113 yh_sup.l = yh_sup.l & 0x7fffffffffffffffLL; /* compute the absolute value*/
1114 yl_sup.l = yl_sup.l & 0x7fffffffffffffffLL; /* compute the absolute value*/
1115 u53_sup.l = (yh_sup.l & ULL(7ff0000000000000)) + ULL(0010000000000000);
1116 roundable = 0;
1117 rd_ok=(yl_inf.d > roundcst * u53_inf.d);
1118 ru_ok=(yl_sup.d > roundcst * u53_sup.d);
1119 if(yl_inf_neg) { /* The case yl==0 is filtered by the above test*/
1120 /* return next down */
1121 yh_inf.d = logh_inf;
1122 if(yh_inf_neg) yh_inf.l++; else yh_inf.l--; /* Beware: fails for zero */
1123 res_inf = yh_inf.d ;
1124 }
1125 else {
1126 res_inf = logh_inf;
1127 }
1128 if(!yl_sup_neg) { /* The case yl==0 is filtered by the above test*/
1129 /* return next up */
1130 yh_sup.d = logh_sup;
1131 if(yh_sup_neg) yh_sup.l--; else yh_sup.l++; /* Beware: fails for zero */
1132 res_sup = yh_sup.d ;
1133 }
1134 else {
1135 res_sup = logh_sup;
1136 }
1137 if(save_res_inf==-1.0/0.0) res_inf=-1.0/0.0;
1138 if(rd_ok && ru_ok){
1139 ASSIGN_LOW(res,res_inf);
1140 ASSIGN_UP(res,res_sup);
1141 return(res);
1142 }
1143 else if (rd_ok){
1144 roundable=1;
1145 }
1146 else if (ru_ok){
1147 roundable=2;
1148 }
1149
1150
1151 #if DEBUG
1152 printf("Going for Accurate Phase for x=%1.50e\n",x);
1153 #endif
1154 if (roundable==1)
1155 {
1156 log_td_accurate(&logh_sup, &logm_sup, &logl_sup, E_sup, ed_sup, index_sup, zh_sup, zl_sup, logih_sup, logim_sup);
1157 RoundUpwards3(&res_sup, logh_sup, logm_sup, logl_sup);
1158 }
1159
1160 if (roundable==2)
1161 {
1162 log_td_accurate(&logh_inf, &logm_inf, &logl_inf, E_inf, ed_inf, index_inf, zh_inf, zl_inf, logih_inf, logim_inf);
1163 RoundDownwards3(&res_inf, logh_inf, logm_inf, logl_inf);
1164 }
1165 if (roundable==0)
1166 {
1167 log_td_accurate(&logh_inf, &logm_inf, &logl_inf, E_inf, ed_inf, index_inf, zh_inf, zl_inf, logih_inf, logim_inf);
1168 RoundDownwards3(&res_inf, logh_inf, logm_inf, logl_inf);
1169 log_td_accurate(&logh_sup, &logm_sup, &logl_sup, E_sup, ed_sup, index_sup, zh_sup, zl_sup, logih_sup, logim_sup);
1170 RoundUpwards3(&res_sup, logh_sup, logm_sup, logl_sup);
1171 }
1172 ASSIGN_LOW(res,res_inf);
1173 ASSIGN_UP(res,res_sup);
1174 return res;
1175 }
1176 #endif
1177
1178