lib/expm1l.c

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