lib/expm1.c

   1 /* Exponential function minus one.
   2    Copyright (C) 2012-2014 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 #include <float.h>
  23
  24 /* A value slightly larger than log(2).  */
  25 #define LOG2_PLUS_EPSILON 0.6931471805599454
  26
  27 /* Best possible approximation of log(2) as a 'double'.  */
  28 #define LOG2 0.693147180559945309417232121458176568075
  29
  30 /* Best possible approximation of 1/log(2) as a 'double'.  */
  31 #define LOG2_INVERSE 1.44269504088896340735992468100189213743
  32
  33 /* Best possible approximation of log(2)/256 as a 'double'.  */
  34 #define LOG2_BY_256 0.00270760617406228636491106297444600221904
  35
  36 /* Best possible approximation of 256/log(2) as a 'double'.  */
  37 #define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181
  38
  39 /* The upper 32 bits of log(2)/256.  */
  40 #define LOG2_BY_256_HI_PART 0.0027076061733168899081647396087646484375
  41 /* log(2)/256 - LOG2_HI_PART.  */
  42 #define LOG2_BY_256_LO_PART \
  43   0.000000000000745396456746323365681353781544922399845
  44
  45 double
  46 expm1 (double x)
  47 {
  48   if (isnand (x))
  49     return x;
  50
  51   if (x >= (double) DBL_MAX_EXP * LOG2_PLUS_EPSILON)
  52     /* x > DBL_MAX_EXP * log(2)
  53        hence exp(x) > 2^DBL_MAX_EXP, overflows to Infinity.  */
  54     return HUGE_VAL;
  55
  56   if (x <= (double) (- DBL_MANT_DIG) * LOG2_PLUS_EPSILON)
  57     /* x < (- DBL_MANT_DIG) * log(2)
  58        hence 0 < exp(x) < 2^-DBL_MANT_DIG,
  59        hence -1 < exp(x)-1 < -1 + 2^-DBL_MANT_DIG
  60        rounds to -1.  */
  61     return -1.0;
  62
  63   if (x <= - LOG2_PLUS_EPSILON)
  64     /* 0 < exp(x) < 1/2.
  65        Just compute exp(x), then subtract 1.  */
  66     return exp (x) - 1.0;
  67
  68   if (x == 0.0)
  69     /* Return a zero with the same sign as x.  */
  70     return x;
  71
  72   /* Decompose x into
  73        x = n * log(2) + m * log(2)/256 + y
  74      where
  75        n is an integer, n >= -1,
  76        m is an integer, -128 <= m <= 128,
  77        y is a number, |y| <= log(2)/512 + epsilon = 0.00135...
  78      Then
  79        exp(x) = 2^n * exp(m * log(2)/256) * exp(y)
  80      Compute each factor minus one, then combine them through the
  81      formula (1+a)*(1+b) = 1 + (a+b*(1+a)),
  82      that is (1+a)*(1+b) - 1 = a + b*(1+a).
  83      The first factor is an ldexpl() call.
  84      The second factor is a table lookup.
  85      The third factor minus one is computed
  86      - either as sinh(y) + sinh(y)^2 / (cosh(y) + 1)
  87        where sinh(y) is computed through the power series:
  88          sinh(y) = y + y^3/3! + y^5/5! + ...
  89        and cosh(y) is computed as hypot(1, sinh(y)),
  90      - or as exp(2*z) - 1 = 2 * tanh(z) / (1 - tanh(z))
  91        where z = y/2
  92        and tanh(z) is computed through its power series:
  93          tanh(z) = z
  94                    - 1/3 * z^3
  95                    + 2/15 * z^5
  96                    - 17/315 * z^7
  97                    + 62/2835 * z^9
  98                    - 1382/155925 * z^11
  99                    + 21844/6081075 * z^13
 100                    - 929569/638512875 * z^15
 101                    + ...
 102        Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the
 103        z^7 term is < 0.0007^6 < 2^-60 <= 2^-DBL_MANT_DIG, therefore we can
 104        truncate the series after the z^5 term.
 105
 106      Given the usual bounds DBL_MAX_EXP <= 16384, DBL_MANT_DIG <= 120, we
 107      can estimate x:  -84 <= x <= 11357.
 108      This means, when dividing x by log(2), where we want x mod log(2)
 109      to be precise to DBL_MANT_DIG bits, we have to use an approximation
 110      to log(2) that has 14+DBL_MANT_DIG bits.  */
 111
 112   {
 113     double nm = round (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */
 114     /* n has at most 15 bits, nm therefore has at most 23 bits, therefore
 115        n * LOG2_HI_PART is computed exactly, and n * LOG2_LO_PART is computed
 116        with an absolute error < 2^15 * 2e-10 * 2^-DBL_MANT_DIG.  */
 117     double y_tmp = x - nm * LOG2_BY_256_HI_PART;
 118     double y = y_tmp - nm * LOG2_BY_256_LO_PART;
 119     double z = 0.5L * y;
 120
 121 /* Coefficients of the power series for tanh(z).  */
 122 #define TANH_COEFF_1   1.0
 123 #define TANH_COEFF_3  -0.333333333333333333333333333333333333334
 124 #define TANH_COEFF_5   0.133333333333333333333333333333333333334
 125 #define TANH_COEFF_7  -0.053968253968253968253968253968253968254
 126 #define TANH_COEFF_9   0.0218694885361552028218694885361552028218
 127 #define TANH_COEFF_11 -0.00886323552990219656886323552990219656886
 128 #define TANH_COEFF_13  0.00359212803657248101692546136990581435026
 129 #define TANH_COEFF_15 -0.00145583438705131826824948518070211191904
 130
 131     double z2 = z * z;
 132     double tanh_z =
 133       ((TANH_COEFF_5
 134         * z2 + TANH_COEFF_3)
 135        * z2 + TANH_COEFF_1)
 136       * z;
 137
 138     double exp_y_minus_1 = 2.0 * tanh_z / (1.0 - tanh_z);
 139
 140     int n = (int) round (nm * (1.0 / 256.0));
 141     int m = (int) nm - 256 * n;
 142
 143     /* expm1_table[i] = exp((i - 128) * log(2)/256) - 1.
 144        Computed in GNU clisp through
 145          (setf (long-float-digits) 128)
 146          (setq a 0L0)
 147          (setf (long-float-digits) 256)
 148          (dotimes (i 257)
 149            (format t "        ~D,~%"
 150                    (float (- (exp (* (/ (- i 128) 256) (log 2L0))) 1) a)))  */
 151     static const double expm1_table[257] =
 152       {
 153         -0.292893218813452475599155637895150960716,
 154         -0.290976057839792401079436677742323809165,
 155         -0.289053698915417220095325702647879950038,
 156         -0.287126127947252846596498423285616993819,
 157         -0.285193330804014994382467110862430046956,
 158         -0.283255293316105578740250215722626632811,
 159         -0.281312001275508837198386957752147486471,
 160         -0.279363440435687168635744042695052413926,
 161         -0.277409596511476689981496879264164547161,
 162         -0.275450455178982509740597294512888729286,
 163         -0.273486002075473717576963754157712706214,
 164         -0.271516222799278089184548475181393238264,
 165         -0.269541102909676505674348554844689233423,
 166         -0.267560627926797086703335317887720824384,
 167         -0.265574783331509036569177486867109287348,
 168         -0.263583554565316202492529493866889713058,
 169         -0.261586927030250344306546259812975038038,
 170         -0.259584886088764114771170054844048746036,
 171         -0.257577417063623749727613604135596844722,
 172         -0.255564505237801467306336402685726757248,
 173         -0.253546135854367575399678234256663229163,
 174         -0.251522294116382286608175138287279137577,
 175         -0.2494929651867872398674385184702356751864,
 176         -0.247458134188296727960327722100283867508,
 177         -0.24541778620328863011699022448340323429,
 178         -0.243371906273695048903181511842366886387,
 179         -0.24132047940089265059510885341281062657,
 180         -0.239263490545592708236869372901757573532,
 181         -0.237200924627730846574373155241529522695,
 182         -0.23513276652635648805745654063657412692,
 183         -0.233059001079521999099699248246140670544,
 184         -0.230979613084171535783261520405692115669,
 185         -0.228894587296029588193854068954632579346,
 186         -0.226803908429489222568744221853864674729,
 187         -0.224707561157500020438486294646580877171,
 188         -0.222605530111455713940842831198332609562,
 189         -0.2204977998810815164831359552625710592544,
 190         -0.218384355014321147927034632426122058645,
 191         -0.2162651800172235534675441445217774245016,
 192         -0.214140259353829315375718509234297186439,
 193         -0.212009577446056756772364919909047495547,
 194         -0.209873118673587736597751517992039478005,
 195         -0.2077308673737531349400659265343210916196,
 196         -0.205582807841418027883101951185666435317,
 197         -0.2034289243288665510313756784404656320656,
 198         -0.201269201045686450868589852895683430425,
 199         -0.199103622158653323103076879204523186316,
 200         -0.196932171791614537151556053482436428417,
 201         -0.19475483402537284591023966632129970827,
 202         -0.192571592897569679960015418424270885733,
 203         -0.190382432402568125350119133273631796029,
 204         -0.188187336491335584102392022226559177731,
 205         -0.185986289071326116575890738992992661386,
 206         -0.183779274006362464829286135533230759947,
 207         -0.181566275116517756116147982921992768975,
 208         -0.17934727617799688564586793151548689933,
 209         -0.1771222609230175777406216376370887771665,
 210         -0.1748912130396911245164132617275148983224,
 211         -0.1726541161719028012138814282020908791644,
 212         -0.170410953919191957302175212789218768074,
 213         -0.168161709836631782476831771511804777363,
 214         -0.165906367434708746670203829291463807099,
 215         -0.1636449101792017131905953879307692887046,
 216         -0.161377321491060724103867675441291294819,
 217         -0.15910358474628545696887452376678510496,
 218         -0.15682368327580335203567701228614769857,
 219         -0.154537600365347409013071332406381692911,
 220         -0.152245319255333652509541396360635796882,
 221         -0.149946823140738265249318713251248832456,
 222         -0.147642095170974388162796469615281683674,
 223         -0.145331118449768586448102562484668501975,
 224         -0.143013876035036980698187522160833990549,
 225         -0.140690350938761042185327811771843747742,
 226         -0.138360526126863051392482883127641270248,
 227         -0.136024384519081218878475585385633792948,
 228         -0.133681908988844467561490046485836530346,
 229         -0.131333082363146875502898959063916619876,
 230         -0.128977887422421778270943284404535317759,
 231         -0.126616306900415529961291721709773157771,
 232         -0.1242483234840609219490048572320697039866,
 233         -0.121873919813350258443919690312343389353,
 234         -0.1194930784812080879189542126763637438278,
 235         -0.11710578203336358947830887503073906297,
 236         -0.1147120129682226132300120925687579825894,
 237         -0.1123117537367393737247203999003383961205,
 238         -0.1099049867422877955201404475637647649574,
 239         -0.1074916943405325099278897180135900838485,
 240         -0.1050718588392995019970556101123417014993,
 241         -0.102645462498446406786148378936109092823,
 242         -0.1002124875297324539725723033374854302454,
 243         -0.097772916096688059846161368344495155786,
 244         -0.0953267303144840657307406742107731280055,
 245         -0.092873912249800621875082699818829828767,
 246         -0.0904144439206957158520284361718212536293,
 247         -0.0879483072964733445019372468353990225585,
 248         -0.0854754842975513284540160873038416459095,
 249         -0.0829959567953287682564584052058555719614,
 250         -0.080509706612053141143695628825336081184,
 251         -0.078016715520687037466429613329061550362,
 252         -0.075516965244774535807472733052603963221,
 253         -0.073010437458307215803773464831151680239,
 254         -0.070497113785589807692349282254427317595,
 255         -0.067976975801105477595185454402763710658,
 256         -0.0654500050293807475554878955602008567352,
 257         -0.06291618294485004933500052502277673278,
 258         -0.0603754909717199109794126487955155117284,
 259         -0.0578279104838327751561896480162548451191,
 260         -0.055273422804530448266460732621318468453,
 261         -0.0527120092065171793298906732865376926237,
 262         -0.0501436509117223676387482401930039000769,
 263         -0.0475683290911628981746625337821392744829,
 264         -0.044986024864805103778829470427200864833,
 265         -0.0423967193014263530636943648520845560749,
 266         -0.0398003934184762630513928111129293882558,
 267         -0.0371970281819375355214808849088086316225,
 268         -0.0345866045061864160477270517354652168038,
 269         -0.0319691032538527747009720477166542375817,
 270         -0.0293445052356798073922893825624102948152,
 271         -0.0267127912103833568278979766786970786276,
 272         -0.0240739418845108520444897665995250062307,
 273         -0.0214279379122998654908388741865642544049,
 274         -0.018774759895536286618755114942929674984,
 275         -0.016114388383412110943633198761985316073,
 276         -0.01344680387238284353202993186779328685225,
 277         -0.0107719868060245158708750409344163322253,
 278         -0.00808991757489031507008688867384418356197,
 279         -0.00540057651636682434752231377783368554176,
 280         -0.00270394391452987374234008615207739887604,
 281         0.0,
 282         0.00271127505020248543074558845036204047301,
 283         0.0054299011128028213513839559347998147001,
 284         0.00815589811841751578309489081720103927357,
 285         0.0108892860517004600204097905618605243881,
 286         0.01363008495148943884025892906393992959584,
 287         0.0163783149109530379404931137862940627635,
 288         0.0191339960777379496848780958207928793998,
 289         0.0218971486541166782344801347832994397821,
 290         0.0246677928971356451482890762708149276281,
 291         0.0274459491187636965388611939222137814994,
 292         0.0302316376860410128717079024539045670944,
 293         0.0330248790212284225001082839704609180866,
 294         0.0358256936019571200299832090180813718441,
 295         0.0386341019613787906124366979546397325796,
 296         0.0414501246883161412645460790118931264803,
 297         0.0442737824274138403219664787399290087847,
 298         0.0471050958792898661299072502271122405627,
 299         0.049944085800687266082038126515907909062,
 300         0.0527907730046263271198912029807463031904,
 301         0.05564517836055715880834132515293865216,
 302         0.0585073227945126901057721096837166450754,
 303         0.0613772272892620809505676780038837262945,
 304         0.0642549128844645497886112570015802206798,
 305         0.0671404006768236181695211209928091626068,
 306         0.070033711820241773542411936757623568504,
 307         0.0729348675259755513850354508738275853402,
 308         0.0758438890627910378032286484760570740623,
 309         0.0787607977571197937406800374384829584908,
 310         0.081685614993215201942115594422531125645,
 311         0.0846183622133092378161051719066143416095,
 312         0.0875590609177696653467978309440397078697,
 313         0.090507732665257659207010655760707978993,
 314         0.0934643990728858542282201462504471620805,
 315         0.096429081816376823386138295859248481766,
 316         0.099401802630221985463696968238829904039,
 317         0.1023825833078409435564142094256468575113,
 318         0.1053714457017412555882746962569503110404,
 319         0.1083684117236786380094236494266198501387,
 320         0.111373503344817603850149254228916637444,
 321         0.1143867425958925363088129569196030678004,
 322         0.1174081515673691990545799630857802666544,
 323         0.120437752409606684429003879866313012766,
 324         0.1234755673330198007337297397753214319548,
 325         0.1265216186082418997947986437870347776336,
 326         0.12957592856628814599726498884024982591,
 327         0.1326385195987192279870737236776230843835,
 328         0.135709414157805514240390330676117013429,
 329         0.1387886347566916537038302838415112547204,
 330         0.14187620396956162271229760828788093894,
 331         0.144972144431804219394413888222915895793,
 332         0.148076478840179006778799662697342680031,
 333         0.15118922995298270581775963520198253612,
 334         0.154310420590216039548221528724806960684,
 335         0.157440073633751029613085766293796821108,
 336         0.160578212027498746369459472576090986253,
 337         0.163724858777577513813573599092185312343,
 338         0.166880036952481570555516298414089287832,
 339         0.1700437696832501880802590357927385730016,
 340         0.1732160801636372475348043545132453888896,
 341         0.176396991650281276284645728483848641053,
 342         0.1795865274628759454861005667694405189764,
 343         0.182784710984341029924457204693850757963,
 344         0.185991565660993831371265649534215563735,
 345         0.189207115002721066717499970560475915293,
 346         0.192431382583151222142727558145431011481,
 347         0.1956643920398273745838370498654519757025,
 348         0.1989061670743804817703025579763002069494,
 349         0.202156731452703142096396957497765876,
 350         0.205416109005123825604211432558411335666,
 351         0.208684323626581577354792255889216998483,
 352         0.211961399276801194468168917732493045449,
 353         0.2152473599804688781165202513387984576236,
 354         0.218542229827408361758207148117394510722,
 355         0.221846032972757516903891841911570785834,
 356         0.225158793637145437709464594384845353705,
 357         0.2284805361068700056940089577927818403626,
 358         0.231811284734075935884556653212794816605,
 359         0.235151063936933305692912507415415760296,
 360         0.238499898199816567833368865859612431546,
 361         0.241857812073484048593677468726595605511,
 362         0.245224830175257932775204967486152674173,
 363         0.248600977189204736621766097302495545187,
 364         0.251986277866316270060206031789203597321,
 365         0.255380757024691089579390657442301194598,
 366         0.258784439549716443077860441815162618762,
 367         0.262197350394250708014010258518416459672,
 368         0.265619514578806324196273999873453036297,
 369         0.269050957191733222554419081032338004715,
 370         0.272491703389402751236692044184602176772,
 371         0.27594177839639210038120243475928938891,
 372         0.279401207505669226913587970027852545961,
 373         0.282870016078778280726669781021514051111,
 374         0.286348229546025533601482208069738348358,
 375         0.289835873406665812232747295491552189677,
 376         0.293332973229089436725559789048704304684,
 377         0.296839554651009665933754117792451159835,
 378         0.300355643379650651014140567070917791291,
 379         0.303881265191935898574523648951997368331,
 380         0.30741644593467724479715157747196172848,
 381         0.310961211524764341922991786330755849366,
 382         0.314515587949354658485983613383997794966,
 383         0.318079601266063994690185647066116617661,
 384         0.321653277603157514326511812330609226158,
 385         0.325236643159741294629537095498721674113,
 386         0.32882972420595439547865089632866510792,
 387         0.33243254708316144935164337949073577407,
 388         0.336045138204145773442627904371869759286,
 389         0.339667524053303005360030669724352576023,
 390         0.343299731186835263824217146181630875424,
 391         0.346941786232945835788173713229537282073,
 392         0.350593715892034391408522196060133960038,
 393         0.354255546936892728298014740140702804344,
 394         0.357927306212901046494536695671766697444,
 395         0.361609020638224755585535938831941474643,
 396         0.365300717204011815430698360337542855432,
 397         0.369002422974590611929601132982192832168,
 398         0.372714165087668369284997857144717215791,
 399         0.376435970754530100216322805518686960261,
 400         0.380167867260238095581945274358283464698,
 401         0.383909881963831954872659527265192818003,
 402         0.387662042298529159042861017950775988895,
 403         0.391424375771926187149835529566243446678,
 404         0.395196909966200178275574599249220994717,
 405         0.398979672538311140209528136715194969206,
 406         0.402772691220204706374713524333378817108,
 407         0.40657599381901544248361973255451684411,
 408         0.410389608217270704414375128268675481146,
 409         0.414213562373095048801688724209698078569
 410       };
 411
 412     double t = expm1_table[128 + m];
 413
 414     /* (1+t) * (1+exp_y_minus_1) - 1 = t + (1+t)*exp_y_minus_1 */
 415     double p_minus_1 = t + (1.0 + t) * exp_y_minus_1;
 416
 417     double s = ldexp (1.0, n) - 1.0;
 418
 419     /* (1+s) * (1+p_minus_1) - 1 = s + (1+s)*p_minus_1 */
 420     return s + (1.0 + s) * p_minus_1;
 421   }
 422 }