1 /* Exponential function.
2 Copyright (C) 2011-2012 Free Software Foundation, Inc.
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.
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.
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/>. */
22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
34 /* A value slightly larger than log(2). */
35 #define LOG2_PLUS_EPSILON 0.6931471805599454L
37 /* Best possible approximation of log(2) as a 'long double'. */
38 #define LOG2 0.693147180559945309417232121458176568075L
40 /* Best possible approximation of 1/log(2) as a 'long double'. */
41 #define LOG2_INVERSE 1.44269504088896340735992468100189213743L
43 /* Best possible approximation of log(2)/256 as a 'long double'. */
44 #define LOG2_BY_256 0.00270760617406228636491106297444600221904L
46 /* Best possible approximation of 256/log(2) as a 'long double'. */
47 #define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L
49 /* The upper 32 bits of log(2)/256. */
50 #define LOG2_BY_256_HI_PART 0.0027076061733168899081647396087646484375L
51 /* log(2)/256 - LOG2_HI_PART. */
52 #define LOG2_BY_256_LO_PART \
53 0.000000000000745396456746323365681353781544922399845L
61 if (x >= (long double) LDBL_MAX_EXP * LOG2_PLUS_EPSILON)
62 /* x > LDBL_MAX_EXP * log(2)
63 hence exp(x) > 2^LDBL_MAX_EXP, overflows to Infinity. */
66 if (x <= (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * LOG2_PLUS_EPSILON)
67 /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * log(2)
68 hence exp(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG),
69 underflows to zero. */
73 x = n * log(2) + m * log(2)/256 + y
76 m is an integer, -128 <= m <= 128,
77 y is a number, |y| <= log(2)/512 + epsilon = 0.00135...
79 exp(x) = 2^n * exp(m * log(2)/256) * exp(y)
80 The first factor is an ldexpl() call.
81 The second factor is a table lookup.
82 The third factor is computed
83 - either as sinh(y) + cosh(y)
84 where sinh(y) is computed through the power series:
85 sinh(y) = y + y^3/3! + y^5/5! + ...
86 and cosh(y) is computed as hypot(1, sinh(y)),
87 - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z))
89 and tanh(z) is computed through its power series:
96 + 21844/6081075 * z^13
97 - 929569/638512875 * z^15
99 Since |z| <= log(2)/1024 < 0.0007, the relative error of the z^13 term
100 is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we can truncate
101 the series after the z^11 term.
103 Given the usual bounds LDBL_MAX_EXP <= 16384, LDBL_MIN_EXP >= -16381,
104 LDBL_MANT_DIG <= 120, we can estimate x: -11440 <= x <= 11357.
105 This means, when dividing x by log(2), where we want x mod log(2)
106 to be precise to LDBL_MANT_DIG bits, we have to use an approximation
107 to log(2) that has 14+LDBL_MANT_DIG bits. */
110 long double nm = roundl (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */
111 /* n has at most 15 bits, nm therefore has at most 23 bits, therefore
112 n * LOG2_HI_PART is computed exactly, and n * LOG2_LO_PART is computed
113 with an absolute error < 2^15 * 2e-10 * 2^-LDBL_MANT_DIG. */
114 long double y_tmp = x - nm * LOG2_BY_256_HI_PART;
115 long double y = y_tmp - nm * LOG2_BY_256_LO_PART;
116 long double z = 0.5L * y;
118 /* Coefficients of the power series for tanh(z). */
119 #define TANH_COEFF_1 1.0L
120 #define TANH_COEFF_3 -0.333333333333333333333333333333333333334L
121 #define TANH_COEFF_5 0.133333333333333333333333333333333333334L
122 #define TANH_COEFF_7 -0.053968253968253968253968253968253968254L
123 #define TANH_COEFF_9 0.0218694885361552028218694885361552028218L
124 #define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L
125 #define TANH_COEFF_13 0.00359212803657248101692546136990581435026L
126 #define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L
128 long double z2 = z * z;
138 long double exp_y = (1.0L + tanh_z) / (1.0L - tanh_z);
140 int n = (int) roundl (nm * (1.0L / 256.0L));
141 int m = (int) nm - 256 * n;
143 /* expl_table[i] = exp((i - 128) * log(2)/256).
144 Computed in GNU clisp through
145 (setf (long-float-digits) 128)
147 (setf (long-float-digits) 256)
150 (float (exp (* (/ (- i 128) 256) (log 2L0))) a))) */
151 static const long double expl_table[257] =
153 0.707106781186547524400844362104849039284L,
154 0.709023942160207598920563322257676190836L,
155 0.710946301084582779904674297352120049962L,
156 0.71287387205274715340350157671438300618L,
157 0.714806669195985005617532889137569953044L,
158 0.71674470668389442125974978427737336719L,
159 0.71868799872449116280161304224785251353L,
160 0.720636559564312831364255957304947586072L,
161 0.72259040348852331001850312073583545284L,
162 0.724549544821017490259402705487111270714L,
163 0.726513997924526282423036245842287293786L,
164 0.728483777200721910815451524818606761737L,
165 0.730458897090323494325651445155310766577L,
166 0.732439372073202913296664682112279175616L,
167 0.734425216668490963430822513132890712652L,
168 0.736416445434683797507470506133110286942L,
169 0.738413072969749655693453740187024961962L,
170 0.740415113911235885228829945155951253966L,
171 0.742422582936376250272386395864403155277L,
172 0.744435494762198532693663597314273242753L,
173 0.746453864145632424600321765743336770838L,
174 0.748477705883617713391824861712720862423L,
175 0.750507034813212760132561481529764324813L,
176 0.752541865811703272039672277899716132493L,
177 0.75458221379671136988300977551659676571L,
178 0.756628093726304951096818488157633113612L,
179 0.75867952059910734940489114658718937343L,
180 0.760736509454407291763130627098242426467L,
181 0.762799075372269153425626844758470477304L,
182 0.76486723347364351194254345936342587308L,
183 0.766940998920478000900300751753859329456L,
184 0.769020386915828464216738479594307884331L,
185 0.771105412703970411806145931045367420652L,
186 0.773196091570510777431255778146135325272L,
187 0.77529243884249997956151370535341912283L,
188 0.777394469888544286059157168801667390437L,
189 0.779502200118918483516864044737428940745L,
190 0.781615644985678852072965367573877941354L,
191 0.783734819982776446532455855478222575498L,
192 0.78585974064617068462428149076570281356L,
193 0.787990422553943243227635080090952504452L,
194 0.790126881326412263402248482007960521995L,
195 0.79226913262624686505993407346567890838L,
196 0.794417192158581972116898048814333564685L,
197 0.796571075671133448968624321559534367934L,
198 0.798730798954313549131410147104316569576L,
199 0.800896377841346676896923120795476813684L,
200 0.803067828208385462848443946517563571584L,
201 0.805245165974627154089760333678700291728L,
202 0.807428407102430320039984581575729114268L,
203 0.809617567597431874649880866726368203972L,
204 0.81181266350866441589760797777344082227L,
205 0.814013710928673883424109261007007338614L,
206 0.816220725993637535170713864466769240053L,
207 0.818433724883482243883852017078007231025L,
208 0.82065272382200311435413206848451310067L,
209 0.822877739076982422259378362362911222833L,
210 0.825108786960308875483586738272485101678L,
211 0.827345883828097198786118571797909120834L,
212 0.829589046080808042697824787210781231927L,
213 0.831838290163368217523168228488195222638L,
214 0.834093632565291253329796170708536192903L,
215 0.836355089820798286809404612069230711295L,
216 0.83862267850893927589613232455870870518L,
217 0.84089641525371454303112547623321489504L,
218 0.84317631672419664796432298771385230143L,
219 0.84546239963465259098692866759361830709L,
220 0.84775468074466634749045860363936420312L,
221 0.850053176859261734750681286748751167545L,
222 0.852357904829025611837203530384718316326L,
223 0.854668881550231413551897437515331498025L,
224 0.856986123964963019301812477839166009452L,
225 0.859309649061238957814672188228156252257L,
226 0.861639473873136948607517116872358729753L,
227 0.863975615480918781121524414614366207052L,
228 0.866318091011155532438509953514163469652L,
229 0.868666917636853124497101040936083380124L,
230 0.871022112577578221729056715595464682243L,
231 0.873383693099584470038708278290226842228L,
232 0.875751676515939078050995142767930296012L,
233 0.878126080186649741556080309687656610647L,
234 0.880506921518791912081045787323636256171L,
235 0.882894217966636410521691124969260937028L,
236 0.885287987031777386769987907431242017412L,
237 0.88768824626326062627527960009966160388L,
238 0.89009501325771220447985955243623523504L,
239 0.892508305659467490072110281986409916153L,
240 0.8949281411607004980029443898876582985L,
241 0.897354537501553593213851621063890907178L,
242 0.899787512470267546027427696662514569756L,
243 0.902227083903311940153838631655504844215L,
244 0.904673269685515934269259325789226871994L,
245 0.907126087750199378124917300181170171233L,
246 0.909585556079304284147971563828178746372L,
247 0.91205169270352665549806275316460097744L,
248 0.914524515702448671545983912696158354092L,
249 0.91700404320467123174354159479414442804L,
250 0.919490293387946858856304371174663918816L,
251 0.921983284479312962533570386670938449637L,
252 0.92448303475522546419252726694739603678L,
253 0.92698956254169278419622653516884831976L,
254 0.929502886214410192307650717745572682403L,
255 0.932023024198894522404814545597236289343L,
256 0.934549994970619252444512104439799143264L,
257 0.93708381705514995066499947497722326722L,
258 0.93962450902828008902058735120448448827L,
259 0.942172089516167224843810351983745154882L,
260 0.944726577195469551733539267378681531548L,
261 0.947287990793482820670109326713462307376L,
262 0.949856349088277632361251759806996099924L,
263 0.952431670908837101825337466217860725517L,
264 0.955013975135194896221170529572799135168L,
265 0.957603280698573646936305635147915443924L,
266 0.960199606581523736948607188887070611744L,
267 0.962802971818062464478519115091191368377L,
268 0.965413395493813583952272948264534783197L,
269 0.968030896746147225299027952283345762418L,
270 0.970655494764320192607710617437589705184L,
271 0.973287208789616643172102023321302921373L,
272 0.97592605811548914795551023340047499377L,
273 0.978572062087700134509161125813435745597L,
274 0.981225240104463713381244885057070325016L,
275 0.983885611616587889056366801238014683926L,
276 0.98655319612761715646797006813220671315L,
277 0.989228013193975484129124959065583667775L,
278 0.99191008242510968492991311132615581644L,
279 0.994599423483633175652477686222166314457L,
280 0.997296056085470126257659913847922601123L,
282 1.00271127505020248543074558845036204047L,
283 1.0054299011128028213513839559347998147L,
284 1.008155898118417515783094890817201039276L,
285 1.01088928605170046002040979056186052439L,
286 1.013630084951489438840258929063939929597L,
287 1.01637831491095303794049311378629406276L,
288 1.0191339960777379496848780958207928794L,
289 1.02189714865411667823448013478329943978L,
290 1.02466779289713564514828907627081492763L,
291 1.0274459491187636965388611939222137815L,
292 1.030231637686041012871707902453904567093L,
293 1.033024879021228422500108283970460918086L,
294 1.035825693601957120029983209018081371844L,
295 1.03863410196137879061243669795463973258L,
296 1.04145012468831614126454607901189312648L,
297 1.044273782427413840321966478739929008784L,
298 1.04710509587928986612990725022711224056L,
299 1.04994408580068726608203812651590790906L,
300 1.05279077300462632711989120298074630319L,
301 1.05564517836055715880834132515293865216L,
302 1.058507322794512690105772109683716645074L,
303 1.061377227289262080950567678003883726294L,
304 1.06425491288446454978861125700158022068L,
305 1.06714040067682361816952112099280916261L,
306 1.0700337118202417735424119367576235685L,
307 1.072934867525975551385035450873827585343L,
308 1.075843889062791037803228648476057074063L,
309 1.07876079775711979374068003743848295849L,
310 1.081685614993215201942115594422531125643L,
311 1.08461836221330923781610517190661434161L,
312 1.087559060917769665346797830944039707867L,
313 1.09050773266525765920701065576070797899L,
314 1.09346439907288585422822014625044716208L,
315 1.096429081816376823386138295859248481766L,
316 1.09940180263022198546369696823882990404L,
317 1.10238258330784094355641420942564685751L,
318 1.10537144570174125558827469625695031104L,
319 1.108368411723678638009423649426619850137L,
320 1.111373503344817603850149254228916637444L,
321 1.1143867425958925363088129569196030678L,
322 1.11740815156736919905457996308578026665L,
323 1.12043775240960668442900387986631301277L,
324 1.123475567333019800733729739775321431954L,
325 1.12652161860824189979479864378703477763L,
326 1.129575928566288145997264988840249825907L,
327 1.13263851959871922798707372367762308438L,
328 1.13570941415780551424039033067611701343L,
329 1.13878863475669165370383028384151125472L,
330 1.14187620396956162271229760828788093894L,
331 1.14497214443180421939441388822291589579L,
332 1.14807647884017900677879966269734268003L,
333 1.15118922995298270581775963520198253612L,
334 1.154310420590216039548221528724806960684L,
335 1.157440073633751029613085766293796821106L,
336 1.16057821202749874636945947257609098625L,
337 1.16372485877757751381357359909218531234L,
338 1.166880036952481570555516298414089287834L,
339 1.170043769683250188080259035792738573L,
340 1.17321608016363724753480435451324538889L,
341 1.176396991650281276284645728483848641054L,
342 1.17958652746287594548610056676944051898L,
343 1.182784710984341029924457204693850757966L,
344 1.18599156566099383137126564953421556374L,
345 1.18920711500272106671749997056047591529L,
346 1.19243138258315122214272755814543101148L,
347 1.195664392039827374583837049865451975705L,
348 1.19890616707438048177030255797630020695L,
349 1.202156731452703142096396957497765876003L,
350 1.205416109005123825604211432558411335666L,
351 1.208684323626581577354792255889216998484L,
352 1.21196139927680119446816891773249304545L,
353 1.215247359980468878116520251338798457624L,
354 1.218542229827408361758207148117394510724L,
355 1.221846032972757516903891841911570785836L,
356 1.225158793637145437709464594384845353707L,
357 1.22848053610687000569400895779278184036L,
358 1.2318112847340759358845566532127948166L,
359 1.235151063936933305692912507415415760294L,
360 1.238499898199816567833368865859612431545L,
361 1.24185781207348404859367746872659560551L,
362 1.24522483017525793277520496748615267417L,
363 1.24860097718920473662176609730249554519L,
364 1.25198627786631627006020603178920359732L,
365 1.255380757024691089579390657442301194595L,
366 1.25878443954971644307786044181516261876L,
367 1.26219735039425070801401025851841645967L,
368 1.265619514578806324196273999873453036296L,
369 1.26905095719173322255441908103233800472L,
370 1.27249170338940275123669204418460217677L,
371 1.27594177839639210038120243475928938891L,
372 1.27940120750566922691358797002785254596L,
373 1.28287001607877828072666978102151405111L,
374 1.286348229546025533601482208069738348355L,
375 1.28983587340666581223274729549155218968L,
376 1.293332973229089436725559789048704304684L,
377 1.296839554651009665933754117792451159835L,
378 1.30035564337965065101414056707091779129L,
379 1.30388126519193589857452364895199736833L,
380 1.30741644593467724479715157747196172848L,
381 1.310961211524764341922991786330755849366L,
382 1.314515587949354658485983613383997794965L,
383 1.318079601266063994690185647066116617664L,
384 1.32165327760315751432651181233060922616L,
385 1.32523664315974129462953709549872167411L,
386 1.32882972420595439547865089632866510792L,
387 1.33243254708316144935164337949073577407L,
388 1.33604513820414577344262790437186975929L,
389 1.33966752405330300536003066972435257602L,
390 1.34329973118683526382421714618163087542L,
391 1.346941786232945835788173713229537282075L,
392 1.35059371589203439140852219606013396004L,
393 1.35425554693689272829801474014070280434L,
394 1.357927306212901046494536695671766697446L,
395 1.36160902063822475558553593883194147464L,
396 1.36530071720401181543069836033754285543L,
397 1.36900242297459061192960113298219283217L,
398 1.37271416508766836928499785714471721579L,
399 1.37643597075453010021632280551868696026L,
400 1.380167867260238095581945274358283464697L,
401 1.383909881963831954872659527265192818L,
402 1.387662042298529159042861017950775988896L,
403 1.39142437577192618714983552956624344668L,
404 1.395196909966200178275574599249220994716L,
405 1.398979672538311140209528136715194969206L,
406 1.40277269122020470637471352433337881711L,
407 1.40657599381901544248361973255451684411L,
408 1.410389608217270704414375128268675481145L,
409 1.41421356237309504880168872420969807857L
412 return ldexpl (expl_table[128 + m] * exp_y, n);