Mercurial > hg > octave-nkf > gnulib-hg

comparison lib/log1p.c @ 16655:e1c84de792a7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
log1p: Provide replacement for Minix and MSVC. * lib/math.in.h (log1p): New declaration. * lib/log1p.c: New file. * m4/log1p.m4: New file. * m4/math_h.m4 (gl_MATH_H): Test whether log1p is declared. (gl_MATH_H_DEFAULTS): Initialize GNULIB_LOG1P, HAVE_LOG1P. * modules/math (Makefile.am): Substitute GNULIB_LOG1P, HAVE_LOG1P. * modules/log1p (Files): Add lib/log1p.c, m4/log1p.m4. (Depends-on): Add math, isnand, log, round. (configure.ac): Invoke gl_FUNC_LOG1P. Arrange to compile replacement if HAVE_LOG1P is 0. * tests/test-math-c++.cc: Check the declaration of log1p. * doc/posix-functions/log1p.texi: Mention the replacement.
author Bruno Haible <bruno@clisp.org>
date Sat, 10 Mar 2012 19:54:08 +0100
parents
children e542fd46ad6f
comparison
equal deleted inserted replaced
16654:537282685581 16655:e1c84de792a7
1 /* Natural logarithm of 1 plus argument.
2 Copyright (C) 2012 Free Software Foundation, Inc.
3
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 3 of the License, or
7 (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18
19 /* Specification. */
20 #include <math.h>
21
22 double
23 log1p (double x)
24 {
25 if (isnand (x))
26 return x;
27
28 if (x <= -1.0)
29 {
30 if (x == -1.0)
31 /* Return -Infinity. */
32 return - HUGE_VAL;
33 else
34 {
35 /* Return NaN. */
36 #if defined _MSC_VER || (defined __sgi && !defined __GNUC__)
37 static double zero;
38 return zero / zero;
39 #else
40 return 0.0 / 0.0;
41 #endif
42 }
43 }
44
45 if (x < -0.5 || x > 1.0)
46 return log (1.0 + x);
47 /* Here -0.5 <= x <= 1.0. */
48
49 if (x == 0.0)
50 /* Return a zero with the same sign as x. */
51 return x;
52
53 /* Decompose x into
54 1 + x = (1 + m/256) * (1 + y)
55 where
56 m is an integer, -128 <= m <= 256,
57 y is a number, |y| <= 1/256.
58 y is computed as
59 y = (256 * x - m) / (256 + m).
60 Then
61 log(1+x) = log(m/256) + log(1+y)
62 The first summand is a table lookup.
63 The second summand is computed
64 - either through the power series
65 log(1+y) = y
66 - 1/2 * y^2
67 + 1/3 * y^3
68 - 1/4 * y^4
69 + 1/5 * y^5
70 - 1/6 * y^6
71 + 1/7 * y^7
72 - 1/8 * y^8
73 + 1/9 * y^9
74 - 1/10 * y^10
75 + 1/11 * y^11
76 - 1/12 * y^12
77 + 1/13 * y^13
78 - 1/14 * y^14
79 + 1/15 * y^15
80 - ...
81 - or as log(1+y) = log((1+z)/(1-z)) = 2 * atanh(z)
82 where z = y/(2+y)
83 and atanh(z) is computed through its power series:
84 atanh(z) = z
85 + 1/3 * z^3
86 + 1/5 * z^5
87 + 1/7 * z^7
88 + 1/9 * z^9
89 + 1/11 * z^11
90 + 1/13 * z^13
91 + 1/15 * z^15
92 + ...
93 Since |z| <= 1/511 < 0.002, the relative contribution of the z^9
94 term is < 1/9*0.002^8 < 2^-60 <= 2^-DBL_MANT_DIG, therefore we
95 can truncate the series after the z^7 term. */
96
97 {
98 double m = round (x * 256.0);
99 double y = ((x * 256.0) - m) / (m + 256.0);
100 double z = y / (2.0 + y);
101
102 /* Coefficients of the power series for atanh(z). */
103 #define ATANH_COEFF_1 1.0
104 #define ATANH_COEFF_3 0.333333333333333333333333333333333333334
105 #define ATANH_COEFF_5 0.2
106 #define ATANH_COEFF_7 0.142857142857142857142857142857142857143
107 #define ATANH_COEFF_9 0.1111111111111111111111111111111111111113
108 #define ATANH_COEFF_11 0.090909090909090909090909090909090909091
109 #define ATANH_COEFF_13 0.076923076923076923076923076923076923077
110 #define ATANH_COEFF_15 0.066666666666666666666666666666666666667
111
112 double z2 = z * z;
113 double atanh_z =
114 (((ATANH_COEFF_7
115 * z2 + ATANH_COEFF_5)
116 * z2 + ATANH_COEFF_3)
117 * z2 + ATANH_COEFF_1)
118 * z;
119
120 /* log_table[i] = log((i + 128) / 256).
121 Computed in GNU clisp through
122 (setf (long-float-digits) 128)
123 (setq a 0L0)
124 (setf (long-float-digits) 256)
125 (dotimes (i 385)
126 (format t " ~D,~%"
127 (float (log (* (/ (+ i 128) 256) 1L0)) a))) */
128 static const double log_table[385] =
129 {
130 -0.693147180559945309417232121458176568075,
131 -0.6853650401178903604697692213970398044,
132 -0.677642994023980055266378075415729732197,
133 -0.669980121278410931188432960495886651496,
134 -0.662375521893191621046203913861404403985,
135 -0.65482831625780871022347679633437927773,
136 -0.647337644528651106250552853843513225963,
137 -0.639902666041133026551361927671647791137,
138 -0.632522558743510466836625989417756304788,
139 -0.625196518651437560022666843685547154042,
140 -0.617923759322357783718626781474514153438,
141 -0.61070351134887071814907205278986876216,
142 -0.60353502187025817679728065207969203929,
143 -0.59641755410139419712166106497071313106,
144 -0.58935038687830174459117031769420187977,
145 -0.582332814219655195222425952134964639978,
146 -0.575364144903561854878438011987654863008,
147 -0.568443702058988073553825606077313299585,
148 -0.561570822771226036828515992768693405624,
149 -0.554744857700826173731906247856527380683,
150 -0.547965170715447412135297057717612244552,
151 -0.541231138534103334345428696561292056747,
152 -0.534542150383306725323860946832334992828,
153 -0.527897607664638146541620672180936254347,
154 -0.52129692363328608707713317540302930314,
155 -0.514739523087127012297831879947234599722,
156 -0.50822484206593331675332852879892694707,
157 -0.50175232756031585480793331389686769463,
158 -0.495321437230025429054660050261215099,
159 -0.488931639131254417913411735261937295862,
160 -0.482582411452595671747679308725825054355,
161 -0.476273242259330949798142595713829069596,
162 -0.470003629245735553650937031148342064701,
163 -0.463773079495099479425751396412036696525,
164 -0.457581109247178400339643902517133157939,
165 -0.451427243672800141272924605544662667972,
166 -0.445311016655364052636629355711651820077,
167 -0.43923197057898186527990882355156990061,
168 -0.4331896561230192424451526269158655235,
169 -0.427183632062807368078106194920633178807,
170 -0.421213465076303550585562626925177406092,
171 -0.415278729556489003230882088534775334993,
172 -0.409379007429300711070330899107921801414,
173 -0.403513887976902632538339065932507598071,
174 -0.397682967666109433030550215403212372894,
175 -0.391885849981783528404356583224421075418,
176 -0.386122145265033447342107580922798666387,
177 -0.380391470556048421030985561769857535915,
178 -0.374693449441410693606984907867576972481,
179 -0.369027711905733333326561361023189215893,
180 -0.363393894187477327602809309537386757124,
181 -0.357791638638807479160052541644010369001,
182 -0.352220593589352099112142921677820359633,
183 -0.346680413213736728498769933032403617363,
184 -0.341170757402767124761784665198737642087,
185 -0.33569129163814153519122263131727209364,
186 -0.330241686870576856279407775480686721935,
187 -0.324821619401237656369001967407777741178,
188 -0.31943077076636122859621528874235306143,
189 -0.314068827624975851026378775827156709194,
190 -0.308735481649613269682442058976885699557,
191 -0.303430429419920096046768517454655701024,
192 -0.298153372319076331310838085093194799765,
193 -0.292904016432932602487907019463045397996,
194 -0.287682072451780927439219005993827431504,
195 -0.282487255574676923482925918282353780414,
196 -0.277319285416234343803903228503274262719,
197 -0.272177885915815673288364959951380595626,
198 -0.267062785249045246292687241862699949179,
199 -0.261973715741573968558059642502581569596,
200 -0.256910413785027239068190798397055267412,
201 -0.251872619755070079927735679796875342712,
202 -0.2468600779315257978846419408385075613265,
203 -0.24187253642048672427253973837916408939,
204 -0.2369097470783577150364265832942468196375,
205 -0.2319714654377751430492321958603212094726,
206 -0.2270574506353460848586128739534071682175,
207 -0.222167465341154296870334265401817316702,
208 -0.2173012756899813951520225351537951559,
209 -0.212458651214193401740613666010165016867,
210 -0.2076393647782445016154410442673876674964,
211 -0.202843192514751471266885961812429707545,
212 -0.1980699137620937948192675366153429027185,
213 -0.193319311003495979595900706211132426563,
214 -0.188591169807550022358923589720001638093,
215 -0.183885278770137362613157202229852743197,
216 -0.179201429457710992616226033183958974965,
217 -0.174539416351899677264255125093377869519,
218 -0.169899036795397472900424896523305726435,
219 -0.165280090939102924303339903679875604517,
220 -0.160682381690473465543308397998034325468,
221 -0.156105714663061654850502877304344269052,
222 -0.1515498981272009378406898175577424691056,
223 -0.1470147429618096590348349122269674042104,
224 -0.142500062607283030157283942253263107981,
225 -0.1380056730194437167017517619422725179055,
226 -0.1335313926245226231463436209313499745895,
227 -0.129077042275142343345847831367985856258,
228 -0.124642445207276597338493356591214304499,
229 -0.1202274269981598003244753948319154994493,
230 -0.115831815525121705099120059938680166568,
231 -0.1114554409253228268966213677328042273655,
232 -0.1070981355563671005131126851708522185606,
233 -0.1027597339577689347753154133345778104976,
234 -0.098440072813252519902888574928971234883,
235 -0.094138990913861910035632096996525066015,
236 -0.0898563291218610470766469347968659624282,
237 -0.0855919303354035139161469686670511961825,
238 -0.0813456394539524058873423550293617843895,
239 -0.077117303344431289769666193261475917783,
240 -0.072906770808087780565737488890929711303,
241 -0.0687138925480518083746933774035034481663,
242 -0.064538521137571171672923915683992928129,
243 -0.0603805109889074798714456529545968095868,
244 -0.0562397183228760777967376942769773768851,
245 -0.0521160011390140183616307870527840213665,
246 -0.0480092191863606077520036253234446621373,
247 -0.0439192339348354905263921515528654458042,
248 -0.0398459085471996706586162402473026835046,
249 -0.0357891078515852792753420982122404025613,
250 -0.0317486983145803011569962827485256299276,
251 -0.0277245480148548604671395114515163869272,
252 -0.0237165266173160421183468505286730579517,
253 -0.0197245053477785891192717326571593033246,
254 -0.015748356968139168607549511460828269521,
255 -0.0117879557520422404691605618900871263399,
256 -0.0078431774610258928731840424909435816546,
257 -0.00391389932113632909231778364357266484272,
258 0.0,
259 0.00389864041565732301393734309584290701073,
260 0.00778214044205494894746290006113676367813,
261 0.01165061721997527413559144280921434893315,
262 0.0155041865359652541508540460424468358779,
263 0.01934296284313093463590553454155047018545,
264 0.0231670592815343782287991609622899165794,
265 0.0269765876982020757480692925396595457815,
266 0.0307716586667536883710282075967721640917,
267 0.0345523815066597334073715005898328652816,
268 0.038318864302136599193755325123797290346,
269 0.042071213920687054375203805926962379448,
270 0.045809536031294203166679267614663342114,
271 0.049533935122276630882096208829824573267,
272 0.0532445145188122828658701937865287769396,
273 0.0569413764001384247590131015404494943015,
274 0.0606246218164348425806061320404202632862,
275 0.0642943507053972572162284502656114944857,
276 0.0679506619085077493945652777726294140346,
277 0.071593653187008817925605272752092034269,
278 0.075223421237587525698605339983662414637,
279 0.078840061707776024531540577859198294559,
280 0.082443669211074591268160068668307805914,
281 0.086034337341803153381797826721996075141,
282 0.0896121586896871326199514693784845287854,
283 0.093177224854183289768781353027759396216,
284 0.096729626458551112295571056487463437015,
285 0.1002694531636751493081301751297276601964,
286 0.1037967936816435648260618037639746883066,
287 0.1073117357890880506671750303711543368066,
288 0.1108143663402901141948061693232119280986,
289 0.1143047712800586336342591448151747734094,
290 0.1177830356563834545387941094705217050686,
291 0.1212492436328696851612122640808405265723,
292 0.1247034785009572358634065153808632684918,
293 0.128145822691930038174109886961074873852,
294 0.1315763577887192725887161286894831624516,
295 0.134995164537504830601983291147085645626,
296 0.138402322859119135685325873601649187393,
297 0.1417979118602573498789527352804727189846,
298 0.1451820098444978972819350637405643235226,
299 0.1485546943231371429098223170672938691604,
300 0.151916042025841975071803424896884511328,
301 0.1552661289111239515223833017101021786436,
302 0.1586050301766385840933711746258415752456,
303 0.161932820269313253240338285123614220592,
304 0.165249572895307162875611449277240313729,
305 0.1685553610298066669415865321701023169345,
306 0.171850256926659222340098946055147264935,
307 0.1751343321278491480142914649863898412374,
308 0.1784076574728182971194002415109419683545,
309 0.181670303107634678260605595617079739242,
310 0.184922338494011992663903592659249621006,
311 0.1881638324181829868259905803105539806714,
312 0.191394852999629454609298807561308873447,
313 0.194615467699671658858138593767269731516,
314 0.1978257433299198803625720711969614690756,
315 0.201025746060590741340908337591797808969,
316 0.204215541428690891503820386196239272214,
317 0.2073951943460705871587455788490062338536,
318 0.210564769107349637669552812732351513721,
319 0.2137243293977181388619051976331987647734,
320 0.216873938300614359619089525744347498479,
321 0.220013658305282095907358638661628360712,
322 0.2231435513142097557662950903098345033745,
323 0.226263678650453389361787082280390161607,
324 0.229374101064845829991480725046139871551,
325 0.232474878743094064920705078095567528222,
326 0.235566071312766909077588218941043410137,
327 0.2386477378501750099171491363522813392526,
328 0.241719936887145168144307515913513900104,
329 0.244782726417690916434704717466314811104,
330 0.247836163904581256780602765746524747999,
331 0.25088030628580941658844644154994089393,
332 0.253915209980963444137323297906606667466,
333 0.256940930897500425446759867911224262093,
334 0.259957524436926066972079494542311044577,
335 0.26296504550088135182072917321108602859,
336 0.265963548497137941339125926537543389269,
337 0.268953087345503958932974357924497845489,
338 0.271933715483641758831669494532999161983,
339 0.274905485872799249167009582983018668293,
340 0.277868451003456306186350032923401233082,
341 0.280822662900887784639519758873134832073,
342 0.28376817313064459834690122235025476666,
343 0.286705032803954314653250930842073965668,
344 0.289633292583042676878893055525668970004,
345 0.292553002686377439978201258664126644308,
346 0.295464212893835876386681906054964195182,
347 0.298366972551797281464900430293496918012,
348 0.301261330578161781012875538233755492657,
349 0.304147335467296717015819874720446989991,
350 0.30702503529491186207512454053537790169,
351 0.309894477722864687861624550833227164546,
352 0.31275571000389688838624655968831903216,
353 0.315608778986303334901366180667483174144,
354 0.318453731118534615810247213590599595595,
355 0.321290612453734292057863145522557457887,
356 0.324119468654211976090670760434987352183,
357 0.326940344995853320592356894073809191681,
358 0.329753286372467981814422811920789810952,
359 0.332558337300076601412275626573419425269,
360 0.335355541921137830257179579814166199074,
361 0.338144944008716397710235913939267433111,
362 0.340926586970593210305089199780356208443,
363 0.34370051385331844468019789211029452987,
364 0.346466767346208580918462188425772950712,
365 0.349225389785288304181275421187371759687,
366 0.35197642315717818465544745625943892599,
367 0.354719909102929028355011218999317665826,
368 0.357455888921803774226009490140904474434,
369 0.360184403575007796281574967493016620926,
370 0.362905493689368453137824345977489846141,
371 0.365619199560964711319396875217046453067,
372 0.368325561158707653048230154050398826898,
373 0.371024618127872663911964910806824955394,
374 0.373716409793584080821016832715823506644,
375 0.376400975164253065997877633436251593315,
376 0.379078352934969458390853345631019858882,
377 0.38174858149084833985966626493567607862,
378 0.384411698910332039734790062481290868519,
379 0.387067742968448287898902502261817665695,
380 0.38971675114002521337046360400352086705,
381 0.392358760602863872479379611988215363485,
382 0.39499380824086897810639403636498176831,
383 0.397621930647138489104829072973405554918,
384 0.40024316412701270692932510199513117008,
385 0.402857544701083514655197565487057707577,
386 0.405465108108164381978013115464349136572,
387 0.408065889808221748430198682969084124381,
388 0.410659924985268385934306203175822787661,
389 0.41324724855021933092547601552548590025,
390 0.415827895143710965613328892954902305356,
391 0.418401899138883817510763261966760106515,
392 0.42096929464412963612886716150679597245,
393 0.423530115505803295718430478017910109426,
394 0.426084395310900063124544879595476618897,
395 0.428632167389698760206812276426639053152,
396 0.43117346481837134085917247895559499848,
397 0.433708320421559393435847903042186017095,
398 0.436236766774918070349041323061121300663,
399 0.438758836207627937745575058511446738878,
400 0.441274560804875229489496441661301225362,
401 0.443783972410300981171768440588146426918,
402 0.446287102628419511532590180619669006749,
403 0.448783982827006710512822115683937186274,
404 0.451274644139458585144692383079012478686,
405 0.453759117467120506644794794442263270651,
406 0.456237433481587594380805538163929748437,
407 0.458709622626976664843883309250877913511,
408 0.461175715122170166367999925597855358603,
409 0.463635740963032513092182277331163919118,
410 0.466089729924599224558619247504769399859,
411 0.468537711563239270375665237462973542708,
412 0.470979715218791012546897856056359251373,
413 0.473415770016672131372578393236978550606,
414 0.475845904869963914265209586304381412175,
415 0.478270148481470280383546145497464809096,
416 0.480688529345751907676618455448011551209,
417 0.48310107575113582273837458485214554795,
418 0.485507815781700807801791077190788900579,
419 0.487908777319238973246173184132656942487,
420 0.490303988045193838150346159645746860531,
421 0.492693475442575255695076950020077845328,
422 0.495077266797851514597964584842833665358,
423 0.497455389202818942250859256731684928918,
424 0.499827869556449329821331415247044141512,
425 0.502194734566715494273584171951812573586,
426 0.504556010752395287058308531738174929982,
427 0.506911724444854354113196312660089270034,
428 0.509261901789807946804074919228323824878,
429 0.51160656874906207851888487520338193135,
430 0.51394575110223431680100608827421759311,
431 0.51627947444845449617281928478756106467,
432 0.518607764208045632152976996364798698556,
433 0.520930645624185312409809834659637709188,
434 0.52324814376454783651680722493487084164,
435 0.525560283522927371382427602307131424923,
436 0.527867089620842385113892217778300963557,
437 0.530168586609121617841419630845212405063,
438 0.532464798869471843873923723460142242606,
439 0.534755750616027675477923292032637111077,
440 0.537041465896883654566729244153832299024,
441 0.539321968595608874655355158077341155752,
442 0.54159728243274437157654230390043409897,
443 0.543867430967283517663338989065998323965,
444 0.546132437598135650382397209231209163864,
445 0.548392325565573162748150286179863158565,
446 0.550647117952662279259948179204913460093,
447 0.552896837686677737580717902230624314327,
448 0.55514150754050159271548035951590405017,
449 0.557381150134006357049816540361233647898,
450 0.559615787935422686270888500526826593487,
451 0.561845443262691817915664819160697456814,
452 0.564070138284802966071384290090190711817,
453 0.566289895023115872590849979337124343595,
454 0.568504735352668712078738764866962263577,
455 0.5707146810034715448536245647415894503,
456 0.572919753561785509092756726626261068625,
457 0.575119974471387940421742546569273429365,
458 0.577315365034823604318112061519496401506,
459 0.579505946414642223855274409488070989814,
460 0.58169173963462248252061075372537234071,
461 0.583872765580982679097413356975291104927,
462 0.586049045003578208904119436287324349516,
463 0.588220598517086043034868221609113995052,
464 0.590387446602176374641916708123598757576,
465 0.59254960960667159874199020959329739696,
466 0.594707107746692789514343546529205333192,
467 0.59685996110779383658731192302565801002,
468 0.59900818964608339938160002446165150206,
469 0.601151813189334836191674317068856441547,
470 0.603290851438084262340585186661310605647,
471 0.6054253239667168894375677681414899356,
472 0.607555250224541795501085152791125371894,
473 0.609680649536855273481833501660588408785,
474 0.611801541105992903529889766428814783686,
475 0.613917944012370492196929119645563790777,
476 0.616029877215514019647565928196700650293,
477 0.618137359555078733872689126674816271683,
478 0.620240409751857528851494632567246856773,
479 0.62233904640877874159710264120869663505,
480 0.62443328801189350104253874405467311991,
481 0.626523152931352759778820859734204069282,
482 0.628608659422374137744308205774183639946,
483 0.6306898256261987050837261409313532241,
484 0.63276666957103782954578646850357975849,
485 0.634839209173010211969493840510489008123,
486 0.63690746223706923162049442718119919119,
487 0.63897144645792072137962398326473680873,
488 0.64103117942093129105560133440539254671,
489 0.643086678603027315392053859585132960477,
490 0.645137961373584701665228496134731905937,
491 0.647185044995309550122320631377863036675,
492 0.64922794662510981889083996990531112227,
493 0.651266683314958103396333353349672108398,
494 0.653301272012745638758615881210873884572,
495 0.65533172956312763209494967856962559648,
496 0.657358072708360030141890023245936165513,
497 0.659380318089127826115336413370955804038,
498 0.661398482245365008260235838709650938148,
499 0.66341258161706625109695030429080128179,
500 0.665422632545090448950092610006660181147,
501 0.667428651271956189947234166318980478403,
502 0.669430653942629267298885270929503510123,
503 0.67142865660530232331713904200189252584,
504 0.67342267521216672029796038880101726475,
505 0.67541272562017673108090414397019748722,
506 0.677398823591806140809682609997348298556,
507 0.67938098479579735014710062847376425181,
508 0.681359224807903068948071559568089441735,
509 0.683333559111620688164363148387750369654,
510 0.68530400309891941654404807896723298642,
511 0.687270572070960267497006884394346103924,
512 0.689233281238808980324914337814603903233,
513 0.691192145724141958859604629216309755938,
514 0.693147180559945309417232121458176568075
515 };
516 return log_table[128 + (int)m] + 2.0 * atanh_z;
517 }
518 }