Mercurial > hg > octave-lyh > gnulib-hg
changeset 16576:6c7c671909d3
expl: Fix precision of computed result.
* lib/expl.c: Completely rewritten.
* modules/expl (Depends-on): Add isnanl, roundl, ldexpl. Remove floorl.
(Maintainer): Add me.
* m4/expl.m4 (gl_FUNC_EXPL): Update computation of EXPL_LIBM.
author | Bruno Haible <bruno@clisp.org> |
---|---|
date | Tue, 06 Mar 2012 03:01:51 +0100 |
parents | 2faf2762788a |
children | 7bea32ad5560 |
files | ChangeLog lib/expl.c m4/expl.m4 modules/expl |
diffstat | 4 files changed, 404 insertions(+), 108 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2012-03-05 Bruno Haible <bruno@clisp.org> + + expl: Fix precision of computed result. + * lib/expl.c: Completely rewritten. + * modules/expl (Depends-on): Add isnanl, roundl, ldexpl. Remove floorl. + (Maintainer): Add me. + * m4/expl.m4 (gl_FUNC_EXPL): Update computation of EXPL_LIBM. + 2012-03-05 Bruno Haible <bruno@clisp.org> cbrt* tests: More tests.
--- a/lib/expl.c +++ b/lib/expl.c @@ -1,9 +1,5 @@ -/* Emulation for expl. - Contributed by Paolo Bonzini - - Copyright 2002-2003, 2007, 2009-2012 Free Software Foundation, Inc. - - This file is part of gnulib. +/* Exponential function. + Copyright (C) 2011-2012 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -33,117 +29,390 @@ #else -/* Code based on glibc/sysdeps/ieee754/ldbl-128/e_expl.c. */ - # include <float.h> -static const long double C[] = { -/* Chebyshev polynomial coefficients for (exp(x)-1)/x */ -# define P1 C[0] -# define P2 C[1] -# define P3 C[2] -# define P4 C[3] -# define P5 C[4] -# define P6 C[5] - 0.5L, - 1.66666666666666666666666666666666683E-01L, - 4.16666666666666666666654902320001674E-02L, - 8.33333333333333333333314659767198461E-03L, - 1.38888888889899438565058018857254025E-03L, - 1.98412698413981650382436541785404286E-04L, +/* A value slightly larger than log(2). */ +#define LOG2_PLUS_EPSILON 0.6931471805599454L + +/* Best possible approximation of log(2) as a 'long double'. */ +#define LOG2 0.693147180559945309417232121458176568075L + +/* Best possible approximation of 1/log(2) as a 'long double'. */ +#define LOG2_INVERSE 1.44269504088896340735992468100189213743L -/* Smallest integer x for which e^x overflows. */ -# define himark C[6] - 11356.523406294143949491931077970765L, +/* Best possible approximation of log(2)/256 as a 'long double'. */ +#define LOG2_BY_256 0.00270760617406228636491106297444600221904L -/* Largest integer x for which e^x underflows. */ -# define lomark C[7] --11433.4627433362978788372438434526231L, +/* Best possible approximation of 256/log(2) as a 'long double'. */ +#define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L -/* very small number */ -# define TINY C[8] - 1.0e-4900L, - -/* 2^16383 */ -# define TWO16383 C[9] - 5.94865747678615882542879663314003565E+4931L}; +/* The upper 32 bits of log(2)/256. */ +#define LOG2_BY_256_HI_PART 0.0027076061733168899081647396087646484375L +/* log(2)/256 - LOG2_HI_PART. */ +#define LOG2_BY_256_LO_PART \ + 0.000000000000745396456746323365681353781544922399845L long double expl (long double x) { - /* Check for usual case. */ - if (x < himark && x > lomark) - { - int exponent; - long double t, x22; - int k = 1; - long double result = 1.0; + if (isnanl (x)) + return x; - /* Compute an integer power of e with a granularity of 0.125. */ - exponent = (int) floorl (x * 8.0L); - x -= exponent / 8.0L; + if (x >= (long double) LDBL_MAX_EXP * LOG2_PLUS_EPSILON) + /* x > LDBL_MAX_EXP * log(2) + hence exp(x) > 2^LDBL_MAX_EXP, overflows to Infinity. */ + return HUGE_VALL; + + if (x <= (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * LOG2_PLUS_EPSILON) + /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * log(2) + hence exp(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG), + underflows to zero. */ + return 0.0L; - if (x > 0.0625) - { - exponent++; - x -= 0.125L; - } + /* Decompose x into + x = n * log(2) + m * log(2)/256 + y + where + n is an integer, + m is an integer, -128 <= m <= 128, + y is a number, |y| <= log(2)/512 + epsilon = 0.00135... + Then + exp(x) = 2^n * exp(m * log(2)/256) * exp(y) + The first factor is an ldexpl() call. + The second factor is a table lookup. + The third factor is computed + - either as sinh(y) + cosh(y) + where sinh(y) is computed through the power series: + sinh(y) = y + y^3/3! + y^5/5! + ... + and cosh(y) is computed as hypot(1, sinh(y)), + - or as exp(2*z) = (1 + tanh(z))^2 / (1 - tanh(z)^2) + where z = y/2 + and tanh(z) is computed through its power series: + tanh(z) = z + - 1/3 * z^3 + + 2/15 * z^5 + - 17/315 * z^7 + + 62/2835 * z^9 + - 1382/155925 * z^11 + + 21844/6081075 * z^13 + - 929569/638512875 * z^15 + + ... + Since |z| <= log(2)/1024 < 0.0007, the relative error of the z^13 term + is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we can truncate + the series after the z^11 term. - if (exponent < 0) - { - t = 0.8824969025845954028648921432290507362220L; /* e^-0.25 */ - exponent = -exponent; - } - else - t = 1.1331484530668263168290072278117938725655L; /* e^0.25 */ + Given the usual bounds LDBL_MAX_EXP <= 16384, LDBL_MIN_EXP >= -16381, + LDBL_MANT_DIG <= 120, we can estimate x: -11440 <= x <= 11357. + This means, when dividing x by log(2), where we want x mod log(2) + to be precise to LDBL_MANT_DIG bits, we have to use an approximation + to log(2) that has 14+LDBL_MANT_DIG bits. */ + + { + long double nm = roundl (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */ + /* n has at most 15 bits, nm therefore has at most 23 bits, therefore + n * LOG2_HI_PART is computed exactly, and n * LOG2_LO_PART is computed + with an absolute error < 2^15 * 2e-10 * 2^-LDBL_MANT_DIG. */ + long double y_tmp = x - nm * LOG2_BY_256_HI_PART; + long double y = y_tmp - nm * LOG2_BY_256_LO_PART; + long double z = 0.5L * y; + +/* Coefficients of the power series for tanh(z). */ +#define TANH_COEFF_1 1.0L +#define TANH_COEFF_3 -0.333333333333333333333333333333333333334L +#define TANH_COEFF_5 0.133333333333333333333333333333333333334L +#define TANH_COEFF_7 -0.053968253968253968253968253968253968254L +#define TANH_COEFF_9 0.0218694885361552028218694885361552028218L +#define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L +#define TANH_COEFF_13 0.00359212803657248101692546136990581435026L +#define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L + + long double z2 = z * z; + long double tanh_z = + (((((TANH_COEFF_11 + * z2 + TANH_COEFF_9) + * z2 + TANH_COEFF_7) + * z2 + TANH_COEFF_5) + * z2 + TANH_COEFF_3) + * z2 + TANH_COEFF_1) + * z; + + long double exp_y = + ((1.0L + tanh_z) * (1.0L + tanh_z)) / (1.0L - tanh_z * tanh_z); + + int n = (int) roundl (nm * (1.0L / 256.0L)); + int m = (int) nm - 256 * n; - while (exponent) - { - if (exponent & k) - { - result *= t; - exponent ^= k; - } - t *= t; - k <<= 1; - } - - /* Approximate (e^x - 1)/x, using a seventh-degree polynomial, - with maximum error in [-2^-16-2^-53,2^-16+2^-53] - less than 4.8e-39. */ - x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6))))); + /* expl_table[i] = exp((i - 128) * log(2)/256). + Computed in GNU clisp through + (progn + (setf (long-float-digits) 128) + (setq a 0L0) + (setf (long-float-digits) 256) + (dotimes (i 257) + (format t " ~D,~%" + (float (exp (* (/ (- i 128) 256) (log 2L0))) a)))) */ + static const long double expl_table[257] = + { + 0.707106781186547524400844362104849039284L, + 0.709023942160207598920563322257676190836L, + 0.710946301084582779904674297352120049962L, + 0.71287387205274715340350157671438300618L, + 0.714806669195985005617532889137569953044L, + 0.71674470668389442125974978427737336719L, + 0.71868799872449116280161304224785251353L, + 0.720636559564312831364255957304947586072L, + 0.72259040348852331001850312073583545284L, + 0.724549544821017490259402705487111270714L, + 0.726513997924526282423036245842287293786L, + 0.728483777200721910815451524818606761737L, + 0.730458897090323494325651445155310766577L, + 0.732439372073202913296664682112279175616L, + 0.734425216668490963430822513132890712652L, + 0.736416445434683797507470506133110286942L, + 0.738413072969749655693453740187024961962L, + 0.740415113911235885228829945155951253966L, + 0.742422582936376250272386395864403155277L, + 0.744435494762198532693663597314273242753L, + 0.746453864145632424600321765743336770838L, + 0.748477705883617713391824861712720862423L, + 0.750507034813212760132561481529764324813L, + 0.752541865811703272039672277899716132493L, + 0.75458221379671136988300977551659676571L, + 0.756628093726304951096818488157633113612L, + 0.75867952059910734940489114658718937343L, + 0.760736509454407291763130627098242426467L, + 0.762799075372269153425626844758470477304L, + 0.76486723347364351194254345936342587308L, + 0.766940998920478000900300751753859329456L, + 0.769020386915828464216738479594307884331L, + 0.771105412703970411806145931045367420652L, + 0.773196091570510777431255778146135325272L, + 0.77529243884249997956151370535341912283L, + 0.777394469888544286059157168801667390437L, + 0.779502200118918483516864044737428940745L, + 0.781615644985678852072965367573877941354L, + 0.783734819982776446532455855478222575498L, + 0.78585974064617068462428149076570281356L, + 0.787990422553943243227635080090952504452L, + 0.790126881326412263402248482007960521995L, + 0.79226913262624686505993407346567890838L, + 0.794417192158581972116898048814333564685L, + 0.796571075671133448968624321559534367934L, + 0.798730798954313549131410147104316569576L, + 0.800896377841346676896923120795476813684L, + 0.803067828208385462848443946517563571584L, + 0.805245165974627154089760333678700291728L, + 0.807428407102430320039984581575729114268L, + 0.809617567597431874649880866726368203972L, + 0.81181266350866441589760797777344082227L, + 0.814013710928673883424109261007007338614L, + 0.816220725993637535170713864466769240053L, + 0.818433724883482243883852017078007231025L, + 0.82065272382200311435413206848451310067L, + 0.822877739076982422259378362362911222833L, + 0.825108786960308875483586738272485101678L, + 0.827345883828097198786118571797909120834L, + 0.829589046080808042697824787210781231927L, + 0.831838290163368217523168228488195222638L, + 0.834093632565291253329796170708536192903L, + 0.836355089820798286809404612069230711295L, + 0.83862267850893927589613232455870870518L, + 0.84089641525371454303112547623321489504L, + 0.84317631672419664796432298771385230143L, + 0.84546239963465259098692866759361830709L, + 0.84775468074466634749045860363936420312L, + 0.850053176859261734750681286748751167545L, + 0.852357904829025611837203530384718316326L, + 0.854668881550231413551897437515331498025L, + 0.856986123964963019301812477839166009452L, + 0.859309649061238957814672188228156252257L, + 0.861639473873136948607517116872358729753L, + 0.863975615480918781121524414614366207052L, + 0.866318091011155532438509953514163469652L, + 0.868666917636853124497101040936083380124L, + 0.871022112577578221729056715595464682243L, + 0.873383693099584470038708278290226842228L, + 0.875751676515939078050995142767930296012L, + 0.878126080186649741556080309687656610647L, + 0.880506921518791912081045787323636256171L, + 0.882894217966636410521691124969260937028L, + 0.885287987031777386769987907431242017412L, + 0.88768824626326062627527960009966160388L, + 0.89009501325771220447985955243623523504L, + 0.892508305659467490072110281986409916153L, + 0.8949281411607004980029443898876582985L, + 0.897354537501553593213851621063890907178L, + 0.899787512470267546027427696662514569756L, + 0.902227083903311940153838631655504844215L, + 0.904673269685515934269259325789226871994L, + 0.907126087750199378124917300181170171233L, + 0.909585556079304284147971563828178746372L, + 0.91205169270352665549806275316460097744L, + 0.914524515702448671545983912696158354092L, + 0.91700404320467123174354159479414442804L, + 0.919490293387946858856304371174663918816L, + 0.921983284479312962533570386670938449637L, + 0.92448303475522546419252726694739603678L, + 0.92698956254169278419622653516884831976L, + 0.929502886214410192307650717745572682403L, + 0.932023024198894522404814545597236289343L, + 0.934549994970619252444512104439799143264L, + 0.93708381705514995066499947497722326722L, + 0.93962450902828008902058735120448448827L, + 0.942172089516167224843810351983745154882L, + 0.944726577195469551733539267378681531548L, + 0.947287990793482820670109326713462307376L, + 0.949856349088277632361251759806996099924L, + 0.952431670908837101825337466217860725517L, + 0.955013975135194896221170529572799135168L, + 0.957603280698573646936305635147915443924L, + 0.960199606581523736948607188887070611744L, + 0.962802971818062464478519115091191368377L, + 0.965413395493813583952272948264534783197L, + 0.968030896746147225299027952283345762418L, + 0.970655494764320192607710617437589705184L, + 0.973287208789616643172102023321302921373L, + 0.97592605811548914795551023340047499377L, + 0.978572062087700134509161125813435745597L, + 0.981225240104463713381244885057070325016L, + 0.983885611616587889056366801238014683926L, + 0.98655319612761715646797006813220671315L, + 0.989228013193975484129124959065583667775L, + 0.99191008242510968492991311132615581644L, + 0.994599423483633175652477686222166314457L, + 0.997296056085470126257659913847922601123L, + 1.0L, + 1.00271127505020248543074558845036204047L, + 1.0054299011128028213513839559347998147L, + 1.008155898118417515783094890817201039276L, + 1.01088928605170046002040979056186052439L, + 1.013630084951489438840258929063939929597L, + 1.01637831491095303794049311378629406276L, + 1.0191339960777379496848780958207928794L, + 1.02189714865411667823448013478329943978L, + 1.02466779289713564514828907627081492763L, + 1.0274459491187636965388611939222137815L, + 1.030231637686041012871707902453904567093L, + 1.033024879021228422500108283970460918086L, + 1.035825693601957120029983209018081371844L, + 1.03863410196137879061243669795463973258L, + 1.04145012468831614126454607901189312648L, + 1.044273782427413840321966478739929008784L, + 1.04710509587928986612990725022711224056L, + 1.04994408580068726608203812651590790906L, + 1.05279077300462632711989120298074630319L, + 1.05564517836055715880834132515293865216L, + 1.058507322794512690105772109683716645074L, + 1.061377227289262080950567678003883726294L, + 1.06425491288446454978861125700158022068L, + 1.06714040067682361816952112099280916261L, + 1.0700337118202417735424119367576235685L, + 1.072934867525975551385035450873827585343L, + 1.075843889062791037803228648476057074063L, + 1.07876079775711979374068003743848295849L, + 1.081685614993215201942115594422531125643L, + 1.08461836221330923781610517190661434161L, + 1.087559060917769665346797830944039707867L, + 1.09050773266525765920701065576070797899L, + 1.09346439907288585422822014625044716208L, + 1.096429081816376823386138295859248481766L, + 1.09940180263022198546369696823882990404L, + 1.10238258330784094355641420942564685751L, + 1.10537144570174125558827469625695031104L, + 1.108368411723678638009423649426619850137L, + 1.111373503344817603850149254228916637444L, + 1.1143867425958925363088129569196030678L, + 1.11740815156736919905457996308578026665L, + 1.12043775240960668442900387986631301277L, + 1.123475567333019800733729739775321431954L, + 1.12652161860824189979479864378703477763L, + 1.129575928566288145997264988840249825907L, + 1.13263851959871922798707372367762308438L, + 1.13570941415780551424039033067611701343L, + 1.13878863475669165370383028384151125472L, + 1.14187620396956162271229760828788093894L, + 1.14497214443180421939441388822291589579L, + 1.14807647884017900677879966269734268003L, + 1.15118922995298270581775963520198253612L, + 1.154310420590216039548221528724806960684L, + 1.157440073633751029613085766293796821106L, + 1.16057821202749874636945947257609098625L, + 1.16372485877757751381357359909218531234L, + 1.166880036952481570555516298414089287834L, + 1.170043769683250188080259035792738573L, + 1.17321608016363724753480435451324538889L, + 1.176396991650281276284645728483848641054L, + 1.17958652746287594548610056676944051898L, + 1.182784710984341029924457204693850757966L, + 1.18599156566099383137126564953421556374L, + 1.18920711500272106671749997056047591529L, + 1.19243138258315122214272755814543101148L, + 1.195664392039827374583837049865451975705L, + 1.19890616707438048177030255797630020695L, + 1.202156731452703142096396957497765876003L, + 1.205416109005123825604211432558411335666L, + 1.208684323626581577354792255889216998484L, + 1.21196139927680119446816891773249304545L, + 1.215247359980468878116520251338798457624L, + 1.218542229827408361758207148117394510724L, + 1.221846032972757516903891841911570785836L, + 1.225158793637145437709464594384845353707L, + 1.22848053610687000569400895779278184036L, + 1.2318112847340759358845566532127948166L, + 1.235151063936933305692912507415415760294L, + 1.238499898199816567833368865859612431545L, + 1.24185781207348404859367746872659560551L, + 1.24522483017525793277520496748615267417L, + 1.24860097718920473662176609730249554519L, + 1.25198627786631627006020603178920359732L, + 1.255380757024691089579390657442301194595L, + 1.25878443954971644307786044181516261876L, + 1.26219735039425070801401025851841645967L, + 1.265619514578806324196273999873453036296L, + 1.26905095719173322255441908103233800472L, + 1.27249170338940275123669204418460217677L, + 1.27594177839639210038120243475928938891L, + 1.27940120750566922691358797002785254596L, + 1.28287001607877828072666978102151405111L, + 1.286348229546025533601482208069738348355L, + 1.28983587340666581223274729549155218968L, + 1.293332973229089436725559789048704304684L, + 1.296839554651009665933754117792451159835L, + 1.30035564337965065101414056707091779129L, + 1.30388126519193589857452364895199736833L, + 1.30741644593467724479715157747196172848L, + 1.310961211524764341922991786330755849366L, + 1.314515587949354658485983613383997794965L, + 1.318079601266063994690185647066116617664L, + 1.32165327760315751432651181233060922616L, + 1.32523664315974129462953709549872167411L, + 1.32882972420595439547865089632866510792L, + 1.33243254708316144935164337949073577407L, + 1.33604513820414577344262790437186975929L, + 1.33966752405330300536003066972435257602L, + 1.34329973118683526382421714618163087542L, + 1.346941786232945835788173713229537282075L, + 1.35059371589203439140852219606013396004L, + 1.35425554693689272829801474014070280434L, + 1.357927306212901046494536695671766697446L, + 1.36160902063822475558553593883194147464L, + 1.36530071720401181543069836033754285543L, + 1.36900242297459061192960113298219283217L, + 1.37271416508766836928499785714471721579L, + 1.37643597075453010021632280551868696026L, + 1.380167867260238095581945274358283464697L, + 1.383909881963831954872659527265192818L, + 1.387662042298529159042861017950775988896L, + 1.39142437577192618714983552956624344668L, + 1.395196909966200178275574599249220994716L, + 1.398979672538311140209528136715194969206L, + 1.40277269122020470637471352433337881711L, + 1.40657599381901544248361973255451684411L, + 1.410389608217270704414375128268675481145L, + 1.41421356237309504880168872420969807857L + }; - return result + result * x22; - } - /* Exceptional cases: */ - else if (x < himark) - { - if (x + x == x) - /* e^-inf == 0, with no error. */ - return 0; - else - /* Underflow */ - return TINY * TINY; - } - else - /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ - return TWO16383*x; + return ldexpl (expl_table[128 + m] * exp_y, n); + } } #endif - -#if 0 -int -main (void) -{ - printf ("%.16Lg\n", expl (1.0L)); - printf ("%.16Lg\n", expl (-1.0L)); - printf ("%.16Lg\n", expl (2.0L)); - printf ("%.16Lg\n", expl (4.0L)); - printf ("%.16Lg\n", expl (-2.0L)); - printf ("%.16Lg\n", expl (-4.0L)); - printf ("%.16Lg\n", expl (0.0625L)); - printf ("%.16Lg\n", expl (0.3L)); - printf ("%.16Lg\n", expl (0.6L)); -} -#endif
--- a/m4/expl.m4 +++ b/m4/expl.m4 @@ -1,4 +1,4 @@ -# expl.m4 serial 6 +# expl.m4 serial 7 dnl Copyright (C) 2010-2012 Free Software Foundation, Inc. dnl This file is free software; the Free Software Foundation dnl gives unlimited permission to copy and/or distribute it, @@ -66,8 +66,25 @@ AC_REQUIRE([gl_FUNC_EXP]) EXPL_LIBM="$EXP_LIBM" else - AC_REQUIRE([gl_FUNC_FLOORL]) - EXPL_LIBM="$FLOORL_LIBM" + AC_REQUIRE([gl_FUNC_ISNANL]) + AC_REQUIRE([gl_FUNC_ROUNDL]) + AC_REQUIRE([gl_FUNC_LDEXPL]) + EXPL_LIBM= + dnl Append $ISNANL_LIBM to EXPL_LIBM, avoiding gratuitous duplicates. + case " $EXPL_LIBM " in + *" $ISNANL_LIBM "*) ;; + *) EXPL_LIBM="$EXPL_LIBM $ISNANL_LIBM" ;; + esac + dnl Append $ROUNDL_LIBM to EXPL_LIBM, avoiding gratuitous duplicates. + case " $EXPL_LIBM " in + *" $ROUNDL_LIBM "*) ;; + *) EXPL_LIBM="$EXPL_LIBM $ROUNDL_LIBM" ;; + esac + dnl Append $LDEXPL_LIBM to EXPL_LIBM, avoiding gratuitous duplicates. + case " $EXPL_LIBM " in + *" $LDEXPL_LIBM "*) ;; + *) EXPL_LIBM="$EXPL_LIBM $LDEXPL_LIBM" ;; + esac fi fi AC_SUBST([EXPL_LIBM])
--- a/modules/expl +++ b/modules/expl @@ -10,7 +10,9 @@ extensions exp [test $HAVE_EXPL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 1] float [test $HAVE_EXPL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0] -floorl [test $HAVE_EXPL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0] +isnanl [test $HAVE_EXPL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0] +roundl [test $HAVE_EXPL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0] +ldexpl [test $HAVE_EXPL = 0 && test $HAVE_SAME_LONG_DOUBLE_AS_DOUBLE = 0] configure.ac: gl_FUNC_EXPL @@ -31,4 +33,4 @@ LGPL Maintainer: -Paolo Bonzini +Bruno Haible