xref: /netbsd-src/lib/libm/ld80/b_logl.c (revision 1bb3b9e15785f3a740f67c6b0ef40b49f73f0ab1)
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