lib/exp2.c

   1 /* Exponential base 2 function.
   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 /* Best possible approximation of log(2) as a 'double'.  */
  25 #define LOG2 0.693147180559945309417232121458176568075
  26
  27 /* Best possible approximation of 1/log(2) as a 'double'.  */
  28 #define LOG2_INVERSE 1.44269504088896340735992468100189213743
  29
  30 /* Best possible approximation of log(2)/256 as a 'double'.  */
  31 #define LOG2_BY_256 0.00270760617406228636491106297444600221904
  32
  33 /* Best possible approximation of 256/log(2) as a 'double'.  */
  34 #define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181
  35
  36 double
  37 exp2 (double x)
  38 {
  39   /* exp2(x) = exp(x*log(2)).
  40      If we would compute it like this, there would be rounding errors for
  41      integer or near-integer values of x.  To avoid these, we inline the
  42      algorithm for exp(), and the multiplication with log(2) cancels a
  43      division by log(2).  */
  44
  45   if (isnand (x))
  46     return x;
  47
  48   if (x > (double) DBL_MAX_EXP)
  49     /* x > DBL_MAX_EXP
  50        hence exp2(x) > 2^DBL_MAX_EXP, overflows to Infinity.  */
  51     return HUGE_VAL;
  52
  53   if (x < (double) (DBL_MIN_EXP - 1 - DBL_MANT_DIG))
  54     /* x < (DBL_MIN_EXP - 1 - DBL_MANT_DIG)
  55        hence exp2(x) < 2^(DBL_MIN_EXP-1-DBL_MANT_DIG),
  56        underflows to zero.  */
  57     return 0.0;
  58
  59   /* Decompose x into
  60        x = n + m/256 + y/log(2)
  61      where
  62        n is an integer,
  63        m is an integer, -128 <= m <= 128,
  64        y is a number, |y| <= log(2)/512 + epsilon = 0.00135...
  65      Then
  66        exp2(x) = 2^n * exp(m * log(2)/256) * exp(y)
  67      The first factor is an ldexpl() call.
  68      The second factor is a table lookup.
  69      The third factor is computed
  70      - either as sinh(y) + cosh(y)
  71        where sinh(y) is computed through the power series:
  72          sinh(y) = y + y^3/3! + y^5/5! + ...
  73        and cosh(y) is computed as hypot(1, sinh(y)),
  74      - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z))
  75        where z = y/2
  76        and tanh(z) is computed through its power series:
  77          tanh(z) = z
  78                    - 1/3 * z^3
  79                    + 2/15 * z^5
  80                    - 17/315 * z^7
  81                    + 62/2835 * z^9
  82                    - 1382/155925 * z^11
  83                    + 21844/6081075 * z^13
  84                    - 929569/638512875 * z^15
  85                    + ...
  86        Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the
  87        z^7 term is < 0.0007^6 < 2^-60 <= 2^-DBL_MANT_DIG, therefore we can
  88        truncate the series after the z^5 term.  */
  89
  90   {
  91     double nm = round (x * 256.0); /* = 256 * n + m */
  92     double z = (x * 256.0 - nm) * (LOG2_BY_256 * 0.5);
  93
  94 /* Coefficients of the power series for tanh(z).  */
  95 #define TANH_COEFF_1   1.0
  96 #define TANH_COEFF_3  -0.333333333333333333333333333333333333334
  97 #define TANH_COEFF_5   0.133333333333333333333333333333333333334
  98 #define TANH_COEFF_7  -0.053968253968253968253968253968253968254
  99 #define TANH_COEFF_9   0.0218694885361552028218694885361552028218
 100 #define TANH_COEFF_11 -0.00886323552990219656886323552990219656886
 101 #define TANH_COEFF_13  0.00359212803657248101692546136990581435026
 102 #define TANH_COEFF_15 -0.00145583438705131826824948518070211191904
 103
 104     double z2 = z * z;
 105     double tanh_z =
 106       ((TANH_COEFF_5
 107         * z2 + TANH_COEFF_3)
 108        * z2 + TANH_COEFF_1)
 109       * z;
 110
 111     double exp_y = (1.0 + tanh_z) / (1.0 - tanh_z);
 112
 113     int n = (int) round (nm * (1.0 / 256.0));
 114     int m = (int) nm - 256 * n;
 115
 116     /* exp_table[i] = exp((i - 128) * log(2)/256).
 117        Computed in GNU clisp through
 118          (setf (long-float-digits) 128)
 119          (setq a 0L0)
 120          (setf (long-float-digits) 256)
 121          (dotimes (i 257)
 122            (format t "        ~D,~%"
 123                    (float (exp (* (/ (- i 128) 256) (log 2L0))) a)))  */
 124     static const double exp_table[257] =
 125       {
 126         0.707106781186547524400844362104849039284,
 127         0.709023942160207598920563322257676190836,
 128         0.710946301084582779904674297352120049962,
 129         0.71287387205274715340350157671438300618,
 130         0.714806669195985005617532889137569953044,
 131         0.71674470668389442125974978427737336719,
 132         0.71868799872449116280161304224785251353,
 133         0.720636559564312831364255957304947586072,
 134         0.72259040348852331001850312073583545284,
 135         0.724549544821017490259402705487111270714,
 136         0.726513997924526282423036245842287293786,
 137         0.728483777200721910815451524818606761737,
 138         0.730458897090323494325651445155310766577,
 139         0.732439372073202913296664682112279175616,
 140         0.734425216668490963430822513132890712652,
 141         0.736416445434683797507470506133110286942,
 142         0.738413072969749655693453740187024961962,
 143         0.740415113911235885228829945155951253966,
 144         0.742422582936376250272386395864403155277,
 145         0.744435494762198532693663597314273242753,
 146         0.746453864145632424600321765743336770838,
 147         0.748477705883617713391824861712720862423,
 148         0.750507034813212760132561481529764324813,
 149         0.752541865811703272039672277899716132493,
 150         0.75458221379671136988300977551659676571,
 151         0.756628093726304951096818488157633113612,
 152         0.75867952059910734940489114658718937343,
 153         0.760736509454407291763130627098242426467,
 154         0.762799075372269153425626844758470477304,
 155         0.76486723347364351194254345936342587308,
 156         0.766940998920478000900300751753859329456,
 157         0.769020386915828464216738479594307884331,
 158         0.771105412703970411806145931045367420652,
 159         0.773196091570510777431255778146135325272,
 160         0.77529243884249997956151370535341912283,
 161         0.777394469888544286059157168801667390437,
 162         0.779502200118918483516864044737428940745,
 163         0.781615644985678852072965367573877941354,
 164         0.783734819982776446532455855478222575498,
 165         0.78585974064617068462428149076570281356,
 166         0.787990422553943243227635080090952504452,
 167         0.790126881326412263402248482007960521995,
 168         0.79226913262624686505993407346567890838,
 169         0.794417192158581972116898048814333564685,
 170         0.796571075671133448968624321559534367934,
 171         0.798730798954313549131410147104316569576,
 172         0.800896377841346676896923120795476813684,
 173         0.803067828208385462848443946517563571584,
 174         0.805245165974627154089760333678700291728,
 175         0.807428407102430320039984581575729114268,
 176         0.809617567597431874649880866726368203972,
 177         0.81181266350866441589760797777344082227,
 178         0.814013710928673883424109261007007338614,
 179         0.816220725993637535170713864466769240053,
 180         0.818433724883482243883852017078007231025,
 181         0.82065272382200311435413206848451310067,
 182         0.822877739076982422259378362362911222833,
 183         0.825108786960308875483586738272485101678,
 184         0.827345883828097198786118571797909120834,
 185         0.829589046080808042697824787210781231927,
 186         0.831838290163368217523168228488195222638,
 187         0.834093632565291253329796170708536192903,
 188         0.836355089820798286809404612069230711295,
 189         0.83862267850893927589613232455870870518,
 190         0.84089641525371454303112547623321489504,
 191         0.84317631672419664796432298771385230143,
 192         0.84546239963465259098692866759361830709,
 193         0.84775468074466634749045860363936420312,
 194         0.850053176859261734750681286748751167545,
 195         0.852357904829025611837203530384718316326,
 196         0.854668881550231413551897437515331498025,
 197         0.856986123964963019301812477839166009452,
 198         0.859309649061238957814672188228156252257,
 199         0.861639473873136948607517116872358729753,
 200         0.863975615480918781121524414614366207052,
 201         0.866318091011155532438509953514163469652,
 202         0.868666917636853124497101040936083380124,
 203         0.871022112577578221729056715595464682243,
 204         0.873383693099584470038708278290226842228,
 205         0.875751676515939078050995142767930296012,
 206         0.878126080186649741556080309687656610647,
 207         0.880506921518791912081045787323636256171,
 208         0.882894217966636410521691124969260937028,
 209         0.885287987031777386769987907431242017412,
 210         0.88768824626326062627527960009966160388,
 211         0.89009501325771220447985955243623523504,
 212         0.892508305659467490072110281986409916153,
 213         0.8949281411607004980029443898876582985,
 214         0.897354537501553593213851621063890907178,
 215         0.899787512470267546027427696662514569756,
 216         0.902227083903311940153838631655504844215,
 217         0.904673269685515934269259325789226871994,
 218         0.907126087750199378124917300181170171233,
 219         0.909585556079304284147971563828178746372,
 220         0.91205169270352665549806275316460097744,
 221         0.914524515702448671545983912696158354092,
 222         0.91700404320467123174354159479414442804,
 223         0.919490293387946858856304371174663918816,
 224         0.921983284479312962533570386670938449637,
 225         0.92448303475522546419252726694739603678,
 226         0.92698956254169278419622653516884831976,
 227         0.929502886214410192307650717745572682403,
 228         0.932023024198894522404814545597236289343,
 229         0.934549994970619252444512104439799143264,
 230         0.93708381705514995066499947497722326722,
 231         0.93962450902828008902058735120448448827,
 232         0.942172089516167224843810351983745154882,
 233         0.944726577195469551733539267378681531548,
 234         0.947287990793482820670109326713462307376,
 235         0.949856349088277632361251759806996099924,
 236         0.952431670908837101825337466217860725517,
 237         0.955013975135194896221170529572799135168,
 238         0.957603280698573646936305635147915443924,
 239         0.960199606581523736948607188887070611744,
 240         0.962802971818062464478519115091191368377,
 241         0.965413395493813583952272948264534783197,
 242         0.968030896746147225299027952283345762418,
 243         0.970655494764320192607710617437589705184,
 244         0.973287208789616643172102023321302921373,
 245         0.97592605811548914795551023340047499377,
 246         0.978572062087700134509161125813435745597,
 247         0.981225240104463713381244885057070325016,
 248         0.983885611616587889056366801238014683926,
 249         0.98655319612761715646797006813220671315,
 250         0.989228013193975484129124959065583667775,
 251         0.99191008242510968492991311132615581644,
 252         0.994599423483633175652477686222166314457,
 253         0.997296056085470126257659913847922601123,
 254         1.0,
 255         1.00271127505020248543074558845036204047,
 256         1.0054299011128028213513839559347998147,
 257         1.008155898118417515783094890817201039276,
 258         1.01088928605170046002040979056186052439,
 259         1.013630084951489438840258929063939929597,
 260         1.01637831491095303794049311378629406276,
 261         1.0191339960777379496848780958207928794,
 262         1.02189714865411667823448013478329943978,
 263         1.02466779289713564514828907627081492763,
 264         1.0274459491187636965388611939222137815,
 265         1.030231637686041012871707902453904567093,
 266         1.033024879021228422500108283970460918086,
 267         1.035825693601957120029983209018081371844,
 268         1.03863410196137879061243669795463973258,
 269         1.04145012468831614126454607901189312648,
 270         1.044273782427413840321966478739929008784,
 271         1.04710509587928986612990725022711224056,
 272         1.04994408580068726608203812651590790906,
 273         1.05279077300462632711989120298074630319,
 274         1.05564517836055715880834132515293865216,
 275         1.058507322794512690105772109683716645074,
 276         1.061377227289262080950567678003883726294,
 277         1.06425491288446454978861125700158022068,
 278         1.06714040067682361816952112099280916261,
 279         1.0700337118202417735424119367576235685,
 280         1.072934867525975551385035450873827585343,
 281         1.075843889062791037803228648476057074063,
 282         1.07876079775711979374068003743848295849,
 283         1.081685614993215201942115594422531125643,
 284         1.08461836221330923781610517190661434161,
 285         1.087559060917769665346797830944039707867,
 286         1.09050773266525765920701065576070797899,
 287         1.09346439907288585422822014625044716208,
 288         1.096429081816376823386138295859248481766,
 289         1.09940180263022198546369696823882990404,
 290         1.10238258330784094355641420942564685751,
 291         1.10537144570174125558827469625695031104,
 292         1.108368411723678638009423649426619850137,
 293         1.111373503344817603850149254228916637444,
 294         1.1143867425958925363088129569196030678,
 295         1.11740815156736919905457996308578026665,
 296         1.12043775240960668442900387986631301277,
 297         1.123475567333019800733729739775321431954,
 298         1.12652161860824189979479864378703477763,
 299         1.129575928566288145997264988840249825907,
 300         1.13263851959871922798707372367762308438,
 301         1.13570941415780551424039033067611701343,
 302         1.13878863475669165370383028384151125472,
 303         1.14187620396956162271229760828788093894,
 304         1.14497214443180421939441388822291589579,
 305         1.14807647884017900677879966269734268003,
 306         1.15118922995298270581775963520198253612,
 307         1.154310420590216039548221528724806960684,
 308         1.157440073633751029613085766293796821106,
 309         1.16057821202749874636945947257609098625,
 310         1.16372485877757751381357359909218531234,
 311         1.166880036952481570555516298414089287834,
 312         1.170043769683250188080259035792738573,
 313         1.17321608016363724753480435451324538889,
 314         1.176396991650281276284645728483848641054,
 315         1.17958652746287594548610056676944051898,
 316         1.182784710984341029924457204693850757966,
 317         1.18599156566099383137126564953421556374,
 318         1.18920711500272106671749997056047591529,
 319         1.19243138258315122214272755814543101148,
 320         1.195664392039827374583837049865451975705,
 321         1.19890616707438048177030255797630020695,
 322         1.202156731452703142096396957497765876003,
 323         1.205416109005123825604211432558411335666,
 324         1.208684323626581577354792255889216998484,
 325         1.21196139927680119446816891773249304545,
 326         1.215247359980468878116520251338798457624,
 327         1.218542229827408361758207148117394510724,
 328         1.221846032972757516903891841911570785836,
 329         1.225158793637145437709464594384845353707,
 330         1.22848053610687000569400895779278184036,
 331         1.2318112847340759358845566532127948166,
 332         1.235151063936933305692912507415415760294,
 333         1.238499898199816567833368865859612431545,
 334         1.24185781207348404859367746872659560551,
 335         1.24522483017525793277520496748615267417,
 336         1.24860097718920473662176609730249554519,
 337         1.25198627786631627006020603178920359732,
 338         1.255380757024691089579390657442301194595,
 339         1.25878443954971644307786044181516261876,
 340         1.26219735039425070801401025851841645967,
 341         1.265619514578806324196273999873453036296,
 342         1.26905095719173322255441908103233800472,
 343         1.27249170338940275123669204418460217677,
 344         1.27594177839639210038120243475928938891,
 345         1.27940120750566922691358797002785254596,
 346         1.28287001607877828072666978102151405111,
 347         1.286348229546025533601482208069738348355,
 348         1.28983587340666581223274729549155218968,
 349         1.293332973229089436725559789048704304684,
 350         1.296839554651009665933754117792451159835,
 351         1.30035564337965065101414056707091779129,
 352         1.30388126519193589857452364895199736833,
 353         1.30741644593467724479715157747196172848,
 354         1.310961211524764341922991786330755849366,
 355         1.314515587949354658485983613383997794965,
 356         1.318079601266063994690185647066116617664,
 357         1.32165327760315751432651181233060922616,
 358         1.32523664315974129462953709549872167411,
 359         1.32882972420595439547865089632866510792,
 360         1.33243254708316144935164337949073577407,
 361         1.33604513820414577344262790437186975929,
 362         1.33966752405330300536003066972435257602,
 363         1.34329973118683526382421714618163087542,
 364         1.346941786232945835788173713229537282075,
 365         1.35059371589203439140852219606013396004,
 366         1.35425554693689272829801474014070280434,
 367         1.357927306212901046494536695671766697446,
 368         1.36160902063822475558553593883194147464,
 369         1.36530071720401181543069836033754285543,
 370         1.36900242297459061192960113298219283217,
 371         1.37271416508766836928499785714471721579,
 372         1.37643597075453010021632280551868696026,
 373         1.380167867260238095581945274358283464697,
 374         1.383909881963831954872659527265192818,
 375         1.387662042298529159042861017950775988896,
 376         1.39142437577192618714983552956624344668,
 377         1.395196909966200178275574599249220994716,
 378         1.398979672538311140209528136715194969206,
 379         1.40277269122020470637471352433337881711,
 380         1.40657599381901544248361973255451684411,
 381         1.410389608217270704414375128268675481145,
 382         1.41421356237309504880168872420969807857
 383       };
 384
 385     return ldexp (exp_table[128 + m] * exp_y, n);
 386   }
 387 }