1cfe182f3Schristos /*-
2cfe182f3Schristos * SPDX-License-Identifier: BSD-3-Clause
3cfe182f3Schristos *
4cfe182f3Schristos * Copyright (c) 1992, 1993
5cfe182f3Schristos * The Regents of the University of California. All rights reserved.
6cfe182f3Schristos *
7cfe182f3Schristos * Redistribution and use in source and binary forms, with or without
8cfe182f3Schristos * modification, are permitted provided that the following conditions
9cfe182f3Schristos * are met:
10cfe182f3Schristos * 1. Redistributions of source code must retain the above copyright
11cfe182f3Schristos * notice, this list of conditions and the following disclaimer.
12cfe182f3Schristos * 2. Redistributions in binary form must reproduce the above copyright
13cfe182f3Schristos * notice, this list of conditions and the following disclaimer in the
14cfe182f3Schristos * documentation and/or other materials provided with the distribution.
15cfe182f3Schristos * 3. Neither the name of the University nor the names of its contributors
16cfe182f3Schristos * may be used to endorse or promote products derived from this software
17cfe182f3Schristos * without specific prior written permission.
18cfe182f3Schristos *
19cfe182f3Schristos * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
20cfe182f3Schristos * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21cfe182f3Schristos * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22cfe182f3Schristos * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
23cfe182f3Schristos * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24cfe182f3Schristos * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25cfe182f3Schristos * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
26cfe182f3Schristos * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27cfe182f3Schristos * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28cfe182f3Schristos * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
29cfe182f3Schristos * SUCH DAMAGE.
30cfe182f3Schristos */
31cfe182f3Schristos
32*1bb3b9e1Skre #include <sys/cdefs.h>
33*1bb3b9e1Skre
34cfe182f3Schristos /*
35cfe182f3Schristos * See bsdsrc/b_log.c for implementation details.
36cfe182f3Schristos *
37cfe182f3Schristos * bsdrc/b_log.c converted to long double by Steven G. Kargl.
38cfe182f3Schristos */
39cfe182f3Schristos
40*1bb3b9e1Skre #include "math_private.h"
41*1bb3b9e1Skre
42cfe182f3Schristos #define N 128
43cfe182f3Schristos
44cfe182f3Schristos /*
45cfe182f3Schristos * Coefficients in the polynomial approximation of log(1+f/F).
46cfe182f3Schristos * Domain of x is [0,1./256] with 2**(-84.48) precision.
47cfe182f3Schristos */
48cfe182f3Schristos static const union ieee_ext_u
49cfe182f3Schristos a1u = LD80C(0xaaaaaaaaaaaaaaab, -4, 8.33333333333333333356e-02L),
50cfe182f3Schristos a2u = LD80C(0xcccccccccccccd29, -7, 1.25000000000000000781e-02L),
51cfe182f3Schristos a3u = LD80C(0x9249249241ed3764, -9, 2.23214285711721994134e-03L),
52cfe182f3Schristos a4u = LD80C(0xe38e959e1e7e01cf, -12, 4.34030476540000360640e-04L);
53cfe182f3Schristos #define A1 (a1u.extu_ld)
54cfe182f3Schristos #define A2 (a2u.extu_ld)
55cfe182f3Schristos #define A3 (a3u.extu_ld)
56cfe182f3Schristos #define A4 (a4u.extu_ld)
57cfe182f3Schristos
58cfe182f3Schristos /*
59cfe182f3Schristos * Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
60cfe182f3Schristos * Used for generation of extend precision logarithms.
61cfe182f3Schristos * The constant 35184372088832 is 2^45, so the divide is exact.
62cfe182f3Schristos * It ensures correct reading of logF_head, even for inaccurate
63cfe182f3Schristos * decimal-to-binary conversion routines. (Everybody gets the
64cfe182f3Schristos * right answer for integers less than 2^53.)
65cfe182f3Schristos * Values for log(F) were generated using error < 10^-57 absolute
66cfe182f3Schristos * with the bc -l package.
67cfe182f3Schristos */
68cfe182f3Schristos
69cfe182f3Schristos static double logF_head[N+1] = {
70cfe182f3Schristos 0.,
71cfe182f3Schristos .007782140442060381246,
72cfe182f3Schristos .015504186535963526694,
73cfe182f3Schristos .023167059281547608406,
74cfe182f3Schristos .030771658666765233647,
75cfe182f3Schristos .038318864302141264488,
76cfe182f3Schristos .045809536031242714670,
77cfe182f3Schristos .053244514518837604555,
78cfe182f3Schristos .060624621816486978786,
79cfe182f3Schristos .067950661908525944454,
80cfe182f3Schristos .075223421237524235039,
81cfe182f3Schristos .082443669210988446138,
82cfe182f3Schristos .089612158689760690322,
83cfe182f3Schristos .096729626458454731618,
84cfe182f3Schristos .103796793681567578460,
85cfe182f3Schristos .110814366340264314203,
86cfe182f3Schristos .117783035656430001836,
87cfe182f3Schristos .124703478501032805070,
88cfe182f3Schristos .131576357788617315236,
89cfe182f3Schristos .138402322859292326029,
90cfe182f3Schristos .145182009844575077295,
91cfe182f3Schristos .151916042025732167530,
92cfe182f3Schristos .158605030176659056451,
93cfe182f3Schristos .165249572895390883786,
94cfe182f3Schristos .171850256926518341060,
95cfe182f3Schristos .178407657472689606947,
96cfe182f3Schristos .184922338493834104156,
97cfe182f3Schristos .191394852999565046047,
98cfe182f3Schristos .197825743329758552135,
99cfe182f3Schristos .204215541428766300668,
100cfe182f3Schristos .210564769107350002741,
101cfe182f3Schristos .216873938300523150246,
102cfe182f3Schristos .223143551314024080056,
103cfe182f3Schristos .229374101064877322642,
104cfe182f3Schristos .235566071312860003672,
105cfe182f3Schristos .241719936886966024758,
106cfe182f3Schristos .247836163904594286577,
107cfe182f3Schristos .253915209980732470285,
108cfe182f3Schristos .259957524436686071567,
109cfe182f3Schristos .265963548496984003577,
110cfe182f3Schristos .271933715484010463114,
111cfe182f3Schristos .277868451003087102435,
112cfe182f3Schristos .283768173130738432519,
113cfe182f3Schristos .289633292582948342896,
114cfe182f3Schristos .295464212893421063199,
115cfe182f3Schristos .301261330578199704177,
116cfe182f3Schristos .307025035294827830512,
117cfe182f3Schristos .312755710004239517729,
118cfe182f3Schristos .318453731118097493890,
119cfe182f3Schristos .324119468654316733591,
120cfe182f3Schristos .329753286372579168528,
121cfe182f3Schristos .335355541920762334484,
122cfe182f3Schristos .340926586970454081892,
123cfe182f3Schristos .346466767346100823488,
124cfe182f3Schristos .351976423156884266063,
125cfe182f3Schristos .357455888922231679316,
126cfe182f3Schristos .362905493689140712376,
127cfe182f3Schristos .368325561158599157352,
128cfe182f3Schristos .373716409793814818840,
129cfe182f3Schristos .379078352934811846353,
130cfe182f3Schristos .384411698910298582632,
131cfe182f3Schristos .389716751140440464951,
132cfe182f3Schristos .394993808240542421117,
133cfe182f3Schristos .400243164127459749579,
134cfe182f3Schristos .405465108107819105498,
135cfe182f3Schristos .410659924985338875558,
136cfe182f3Schristos .415827895143593195825,
137cfe182f3Schristos .420969294644237379543,
138cfe182f3Schristos .426084395310681429691,
139cfe182f3Schristos .431173464818130014464,
140cfe182f3Schristos .436236766774527495726,
141cfe182f3Schristos .441274560805140936281,
142cfe182f3Schristos .446287102628048160113,
143cfe182f3Schristos .451274644139630254358,
144cfe182f3Schristos .456237433481874177232,
145cfe182f3Schristos .461175715122408291790,
146cfe182f3Schristos .466089729924533457960,
147cfe182f3Schristos .470979715219073113985,
148cfe182f3Schristos .475845904869856894947,
149cfe182f3Schristos .480688529345570714212,
150cfe182f3Schristos .485507815781602403149,
151cfe182f3Schristos .490303988045525329653,
152cfe182f3Schristos .495077266798034543171,
153cfe182f3Schristos .499827869556611403822,
154cfe182f3Schristos .504556010751912253908,
155cfe182f3Schristos .509261901790523552335,
156cfe182f3Schristos .513945751101346104405,
157cfe182f3Schristos .518607764208354637958,
158cfe182f3Schristos .523248143765158602036,
159cfe182f3Schristos .527867089620485785417,
160cfe182f3Schristos .532464798869114019908,
161cfe182f3Schristos .537041465897345915436,
162cfe182f3Schristos .541597282432121573947,
163cfe182f3Schristos .546132437597407260909,
164cfe182f3Schristos .550647117952394182793,
165cfe182f3Schristos .555141507540611200965,
166cfe182f3Schristos .559615787935399566777,
167cfe182f3Schristos .564070138285387656651,
168cfe182f3Schristos .568504735352689749561,
169cfe182f3Schristos .572919753562018740922,
170cfe182f3Schristos .577315365035246941260,
171cfe182f3Schristos .581691739635061821900,
172cfe182f3Schristos .586049045003164792433,
173cfe182f3Schristos .590387446602107957005,
174cfe182f3Schristos .594707107746216934174,
175cfe182f3Schristos .599008189645246602594,
176cfe182f3Schristos .603290851438941899687,
177cfe182f3Schristos .607555250224322662688,
178cfe182f3Schristos .611801541106615331955,
179cfe182f3Schristos .616029877215623855590,
180cfe182f3Schristos .620240409751204424537,
181cfe182f3Schristos .624433288012369303032,
182cfe182f3Schristos .628608659422752680256,
183cfe182f3Schristos .632766669570628437213,
184cfe182f3Schristos .636907462236194987781,
185cfe182f3Schristos .641031179420679109171,
186cfe182f3Schristos .645137961373620782978,
187cfe182f3Schristos .649227946625615004450,
188cfe182f3Schristos .653301272011958644725,
189cfe182f3Schristos .657358072709030238911,
190cfe182f3Schristos .661398482245203922502,
191cfe182f3Schristos .665422632544505177065,
192cfe182f3Schristos .669430653942981734871,
193cfe182f3Schristos .673422675212350441142,
194cfe182f3Schristos .677398823590920073911,
195cfe182f3Schristos .681359224807238206267,
196cfe182f3Schristos .685304003098281100392,
197cfe182f3Schristos .689233281238557538017,
198cfe182f3Schristos .693147180560117703862
199cfe182f3Schristos };
200cfe182f3Schristos
201cfe182f3Schristos static double logF_tail[N+1] = {
202cfe182f3Schristos 0.,
203cfe182f3Schristos -.00000000000000543229938420049,
204cfe182f3Schristos .00000000000000172745674997061,
205cfe182f3Schristos -.00000000000001323017818229233,
206cfe182f3Schristos -.00000000000001154527628289872,
207cfe182f3Schristos -.00000000000000466529469958300,
208cfe182f3Schristos .00000000000005148849572685810,
209cfe182f3Schristos -.00000000000002532168943117445,
210cfe182f3Schristos -.00000000000005213620639136504,
211cfe182f3Schristos -.00000000000001819506003016881,
212cfe182f3Schristos .00000000000006329065958724544,
213cfe182f3Schristos .00000000000008614512936087814,
214cfe182f3Schristos -.00000000000007355770219435028,
215cfe182f3Schristos .00000000000009638067658552277,
216cfe182f3Schristos .00000000000007598636597194141,
217cfe182f3Schristos .00000000000002579999128306990,
218cfe182f3Schristos -.00000000000004654729747598444,
219cfe182f3Schristos -.00000000000007556920687451336,
220cfe182f3Schristos .00000000000010195735223708472,
221cfe182f3Schristos -.00000000000017319034406422306,
222cfe182f3Schristos -.00000000000007718001336828098,
223cfe182f3Schristos .00000000000010980754099855238,
224cfe182f3Schristos -.00000000000002047235780046195,
225cfe182f3Schristos -.00000000000008372091099235912,
226cfe182f3Schristos .00000000000014088127937111135,
227cfe182f3Schristos .00000000000012869017157588257,
228cfe182f3Schristos .00000000000017788850778198106,
229cfe182f3Schristos .00000000000006440856150696891,
230cfe182f3Schristos .00000000000016132822667240822,
231cfe182f3Schristos -.00000000000007540916511956188,
232cfe182f3Schristos -.00000000000000036507188831790,
233cfe182f3Schristos .00000000000009120937249914984,
234cfe182f3Schristos .00000000000018567570959796010,
235cfe182f3Schristos -.00000000000003149265065191483,
236cfe182f3Schristos -.00000000000009309459495196889,
237cfe182f3Schristos .00000000000017914338601329117,
238cfe182f3Schristos -.00000000000001302979717330866,
239cfe182f3Schristos .00000000000023097385217586939,
240cfe182f3Schristos .00000000000023999540484211737,
241cfe182f3Schristos .00000000000015393776174455408,
242cfe182f3Schristos -.00000000000036870428315837678,
243cfe182f3Schristos .00000000000036920375082080089,
244cfe182f3Schristos -.00000000000009383417223663699,
245cfe182f3Schristos .00000000000009433398189512690,
246cfe182f3Schristos .00000000000041481318704258568,
247cfe182f3Schristos -.00000000000003792316480209314,
248cfe182f3Schristos .00000000000008403156304792424,
249cfe182f3Schristos -.00000000000034262934348285429,
250cfe182f3Schristos .00000000000043712191957429145,
251cfe182f3Schristos -.00000000000010475750058776541,
252cfe182f3Schristos -.00000000000011118671389559323,
253cfe182f3Schristos .00000000000037549577257259853,
254cfe182f3Schristos .00000000000013912841212197565,
255cfe182f3Schristos .00000000000010775743037572640,
256cfe182f3Schristos .00000000000029391859187648000,
257cfe182f3Schristos -.00000000000042790509060060774,
258cfe182f3Schristos .00000000000022774076114039555,
259cfe182f3Schristos .00000000000010849569622967912,
260cfe182f3Schristos -.00000000000023073801945705758,
261cfe182f3Schristos .00000000000015761203773969435,
262cfe182f3Schristos .00000000000003345710269544082,
263cfe182f3Schristos -.00000000000041525158063436123,
264cfe182f3Schristos .00000000000032655698896907146,
265cfe182f3Schristos -.00000000000044704265010452446,
266cfe182f3Schristos .00000000000034527647952039772,
267cfe182f3Schristos -.00000000000007048962392109746,
268cfe182f3Schristos .00000000000011776978751369214,
269cfe182f3Schristos -.00000000000010774341461609578,
270cfe182f3Schristos .00000000000021863343293215910,
271cfe182f3Schristos .00000000000024132639491333131,
272cfe182f3Schristos .00000000000039057462209830700,
273cfe182f3Schristos -.00000000000026570679203560751,
274cfe182f3Schristos .00000000000037135141919592021,
275cfe182f3Schristos -.00000000000017166921336082431,
276cfe182f3Schristos -.00000000000028658285157914353,
277cfe182f3Schristos -.00000000000023812542263446809,
278cfe182f3Schristos .00000000000006576659768580062,
279cfe182f3Schristos -.00000000000028210143846181267,
280cfe182f3Schristos .00000000000010701931762114254,
281cfe182f3Schristos .00000000000018119346366441110,
282cfe182f3Schristos .00000000000009840465278232627,
283cfe182f3Schristos -.00000000000033149150282752542,
284cfe182f3Schristos -.00000000000018302857356041668,
285cfe182f3Schristos -.00000000000016207400156744949,
286cfe182f3Schristos .00000000000048303314949553201,
287cfe182f3Schristos -.00000000000071560553172382115,
288cfe182f3Schristos .00000000000088821239518571855,
289cfe182f3Schristos -.00000000000030900580513238244,
290cfe182f3Schristos -.00000000000061076551972851496,
291cfe182f3Schristos .00000000000035659969663347830,
292cfe182f3Schristos .00000000000035782396591276383,
293cfe182f3Schristos -.00000000000046226087001544578,
294cfe182f3Schristos .00000000000062279762917225156,
295cfe182f3Schristos .00000000000072838947272065741,
296cfe182f3Schristos .00000000000026809646615211673,
297cfe182f3Schristos -.00000000000010960825046059278,
298cfe182f3Schristos .00000000000002311949383800537,
299cfe182f3Schristos -.00000000000058469058005299247,
300cfe182f3Schristos -.00000000000002103748251144494,
301cfe182f3Schristos -.00000000000023323182945587408,
302cfe182f3Schristos -.00000000000042333694288141916,
303cfe182f3Schristos -.00000000000043933937969737844,
304cfe182f3Schristos .00000000000041341647073835565,
305cfe182f3Schristos .00000000000006841763641591466,
306cfe182f3Schristos .00000000000047585534004430641,
307cfe182f3Schristos .00000000000083679678674757695,
308cfe182f3Schristos -.00000000000085763734646658640,
309cfe182f3Schristos .00000000000021913281229340092,
310cfe182f3Schristos -.00000000000062242842536431148,
311cfe182f3Schristos -.00000000000010983594325438430,
312cfe182f3Schristos .00000000000065310431377633651,
313cfe182f3Schristos -.00000000000047580199021710769,
314cfe182f3Schristos -.00000000000037854251265457040,
315cfe182f3Schristos .00000000000040939233218678664,
316cfe182f3Schristos .00000000000087424383914858291,
317cfe182f3Schristos .00000000000025218188456842882,
318cfe182f3Schristos -.00000000000003608131360422557,
319cfe182f3Schristos -.00000000000050518555924280902,
320cfe182f3Schristos .00000000000078699403323355317,
321cfe182f3Schristos -.00000000000067020876961949060,
322cfe182f3Schristos .00000000000016108575753932458,
323cfe182f3Schristos .00000000000058527188436251509,
324cfe182f3Schristos -.00000000000035246757297904791,
325cfe182f3Schristos -.00000000000018372084495629058,
326cfe182f3Schristos .00000000000088606689813494916,
327cfe182f3Schristos .00000000000066486268071468700,
328cfe182f3Schristos .00000000000063831615170646519,
329cfe182f3Schristos .00000000000025144230728376072,
330cfe182f3Schristos -.00000000000017239444525614834
331cfe182f3Schristos };
332cfe182f3Schristos /*
333cfe182f3Schristos * Extra precision variant, returning struct {double a, b;};
334cfe182f3Schristos * log(x) = a + b to 63 bits, with 'a' rounded to 24 bits.
335cfe182f3Schristos */
336cfe182f3Schristos static struct LDouble
__log__LD(long double x)337cfe182f3Schristos __log__LD(long double x)
338cfe182f3Schristos {
339cfe182f3Schristos int m, j;
340cfe182f3Schristos long double F, f, g, q, u, v, u1, u2;
341cfe182f3Schristos struct LDouble r;
342cfe182f3Schristos
343cfe182f3Schristos /*
344cfe182f3Schristos * Argument reduction: 1 <= g < 2; x/2^m = g;
345cfe182f3Schristos * y = F*(1 + f/F) for |f| <= 2^-8
346cfe182f3Schristos */
347cfe182f3Schristos g = frexpl(x, &m);
348cfe182f3Schristos g *= 2;
349cfe182f3Schristos m--;
350cfe182f3Schristos if (m == DBL_MIN_EXP - 1) {
351cfe182f3Schristos j = ilogbl(g);
352cfe182f3Schristos m += j;
353cfe182f3Schristos g = ldexpl(g, -j);
354cfe182f3Schristos }
355cfe182f3Schristos j = N * (g - 1) + 0.5L;
356cfe182f3Schristos F = (1.L / N) * j + 1;
357cfe182f3Schristos f = g - F;
358cfe182f3Schristos
359cfe182f3Schristos g = 1 / (2 * F + f);
360cfe182f3Schristos u = 2 * f * g;
361cfe182f3Schristos v = u * u;
362cfe182f3Schristos q = u * v * (A1 + v * (A2 + v * (A3 + v * A4)));
363cfe182f3Schristos if (m | j) {
364cfe182f3Schristos u1 = u + 513;
365cfe182f3Schristos u1 -= 513;
366cfe182f3Schristos } else {
367cfe182f3Schristos u1 = (float)u;
368cfe182f3Schristos }
369cfe182f3Schristos u2 = (2 * (f - F * u1) - u1 * f) * g;
370cfe182f3Schristos
371cfe182f3Schristos u1 += m * (long double)logF_head[N] + logF_head[j];
372cfe182f3Schristos
373cfe182f3Schristos u2 += logF_tail[j];
374cfe182f3Schristos u2 += q;
375cfe182f3Schristos u2 += logF_tail[N] * m;
376cfe182f3Schristos r.a = (float)(u1 + u2); /* Only difference is here. */
377cfe182f3Schristos r.b = (u1 - r.a) + u2;
378cfe182f3Schristos return (r);
379cfe182f3Schristos }
380