1 /*
2 * Correctly rounded logarithm
3 *
4 * Copyright (C) 2004-2011 David Defour, Catherine Daramy-Loirat,
5 * Florent de Dinechin, Matthieu Gallet, Nicolas Gast, Christoph Quirin Lauter,
6 * and Jean-Michel Muller
7 *
8 * Author: David Defour
9 *
10 * This file is part of crlibm, the correctly rounded mathematical library,
11 * which has been developed by the Arénaire project at École normale supérieure
12 * de Lyon.
13 *
14 * This library is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU Lesser General Public
16 * License as published by the Free Software Foundation; either
17 * version 2.1 of the License, or (at your option) any later version.
18 *
19 * This library is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * Lesser General Public License for more details.
23 *
24 * You should have received a copy of the GNU Lesser General Public
25 * License along with this library; if not, write to the Free Software
26 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
27 */
28
29
30 #include <stdio.h>
31 #include "log_accurate.h"
32
33 /*
34 * 1) First reduction: exponent extraction
35 * E
36 * x = 2^ .(1+f) with 0 <= f < 1
37 *
38 * log(x) = E.log(2) + log(1+f) where:
39 * - log(2) is tabulated
40 * - log(1+f) need to be evaluated
41 *
42 *
43 * 2) Avoiding accuracy problem when E=-1 by testing
44 *
45 * if (1+f >= sqrt(2)) then
46 * 1+f = (1+f)/2; E = E+1;
47 * and,
48 * log(x) = (E+1).log(2) + log((1+f)/2)
49 *
50 * so now: sqrt(2)/2 <= (1+f) < sqrt(2)
51 *
52 *
53 * 3) Second reduction: tabular reduction
54 * -4
55 * wi = 1 + i. 2^
56 * 1
57 * log(1+f) = log(wi) + log ( 1 + --- . (1 + f - wi) )
58 * wi
59 *
60 * then |(1+f-wi)/wi| <= 2^-5 if we use rounded to nearest.
61 *
62 * 4) Computation:
63 * a) Table lookup of:
64 * - ti = log(wi)
65 * - inv_wi = 1/(wi)
66 * b) Polynomial evaluation of:
67 * - P(R) ~ log(1 + R), where R = (1+f-wi) * inv_wi
68 *
69 * -5
70 * with |R| < 2^
71 *
72 *
73 * 5) Reconstruction:
74 * log(x) = E.log(2) + t_i + P(R)
75 *
76 */
77
78
79
80
scs_log(scs_ptr res,db_number y,int E)81 void scs_log(scs_ptr res, db_number y, int E){
82 scs_t R, sc_ln2_times_E, res1, addi;
83 scs_ptr ti, inv_wi;
84 db_number z, wi;
85 int i;
86
87
88 #if EVAL_PERF
89 crlibm_second_step_taken++;
90 #endif
91
92
93 /* to normalize y.d and round to nearest */
94 /* + (1-trunc(sqrt(2.)/2 * 2^(4))*2^(-4) )+2.^(-(4+1))*/
95 z.d = y.d + norm_number.d;
96 i = (z.i[HI] & 0x000fffff);
97 i = i >> 16; /* 0<= i <=11 */
98
99
100 wi.d = ((double)(11+i))*0.0625;
101
102 /* (1+f-w_i) */
103 y.d -= wi.d;
104
105 /* Table reduction */
106 ti = table_ti_ptr[i];
107 inv_wi = table_inv_wi_ptr[i];
108
109 /* R = (1+f-w_i)/w_i */
110 scs_set_d(R, y.d);
111 scs_mul(R, R, inv_wi);
112
113
114 /*
115 * Polynomial evaluation of log(1 + R) with an error less than 2^(-130)
116 */
117
118 scs_mul(res1, constant_poly_ptr[0], R);
119 for(i=1; i<20; i++){
120 scs_add(addi, constant_poly_ptr[i], res1);
121 scs_mul(res1, addi, R);
122 }
123
124 if(E==0){
125 scs_add(res, res1, ti);
126 }else{
127 /* sc_ln2_times_E = E*log(2) */
128 scs_set(sc_ln2_times_E, sc_ln2_ptr);
129
130 if (E >= 0){
131 scs_mul_ui(sc_ln2_times_E, (unsigned int) E);
132 }else{
133 scs_mul_ui(sc_ln2_times_E, (unsigned int) -E);
134 sc_ln2_times_E->sign = -1;
135 }
136 scs_add(addi, res1, ti);
137 scs_add(res, addi, sc_ln2_times_E);
138 }
139 }
140
141