11573Srgrimes/* 21573Srgrimes * Copyright (c) 1992, 1993 31573Srgrimes * The Regents of the University of California. All rights reserved. 41573Srgrimes * 51573Srgrimes * Redistribution and use in source and binary forms, with or without 61573Srgrimes * modification, are permitted provided that the following conditions 71573Srgrimes * are met: 81573Srgrimes * 1. Redistributions of source code must retain the above copyright 91573Srgrimes * notice, this list of conditions and the following disclaimer. 101573Srgrimes * 2. Redistributions in binary form must reproduce the above copyright 111573Srgrimes * notice, this list of conditions and the following disclaimer in the 121573Srgrimes * documentation and/or other materials provided with the distribution. 131573Srgrimes * 3. All advertising materials mentioning features or use of this software 141573Srgrimes * must display the following acknowledgement: 151573Srgrimes * This product includes software developed by the University of 161573Srgrimes * California, Berkeley and its contributors. 171573Srgrimes * 4. Neither the name of the University nor the names of its contributors 181573Srgrimes * may be used to endorse or promote products derived from this software 191573Srgrimes * without specific prior written permission. 201573Srgrimes * 211573Srgrimes * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 221573Srgrimes * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 231573Srgrimes * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 241573Srgrimes * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 251573Srgrimes * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 261573Srgrimes * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 271573Srgrimes * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 281573Srgrimes * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 291573Srgrimes * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 301573Srgrimes * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 311573Srgrimes * SUCH DAMAGE. 321573Srgrimes */ 331573Srgrimes 34176449Sdas/* @(#)log.c 8.2 (Berkeley) 11/30/93 */ 3592887Sobrien#include <sys/cdefs.h> 3692887Sobrien__FBSDID("$FreeBSD$"); 371573Srgrimes 381573Srgrimes#include <math.h> 391573Srgrimes#include <errno.h> 401573Srgrimes 411573Srgrimes#include "mathimpl.h" 421573Srgrimes 431573Srgrimes/* Table-driven natural logarithm. 441573Srgrimes * 451573Srgrimes * This code was derived, with minor modifications, from: 461573Srgrimes * Peter Tang, "Table-Driven Implementation of the 471573Srgrimes * Logarithm in IEEE Floating-Point arithmetic." ACM Trans. 481573Srgrimes * Math Software, vol 16. no 4, pp 378-400, Dec 1990). 491573Srgrimes * 501573Srgrimes * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256, 511573Srgrimes * where F = j/128 for j an integer in [0, 128]. 521573Srgrimes * 531573Srgrimes * log(2^m) = log2_hi*m + log2_tail*m 541573Srgrimes * since m is an integer, the dominant term is exact. 551573Srgrimes * m has at most 10 digits (for subnormal numbers), 561573Srgrimes * and log2_hi has 11 trailing zero bits. 571573Srgrimes * 581573Srgrimes * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h 591573Srgrimes * logF_hi[] + 512 is exact. 601573Srgrimes * 611573Srgrimes * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ... 621573Srgrimes * the leading term is calculated to extra precision in two 631573Srgrimes * parts, the larger of which adds exactly to the dominant 641573Srgrimes * m and F terms. 651573Srgrimes * There are two cases: 661573Srgrimes * 1. when m, j are non-zero (m | j), use absolute 671573Srgrimes * precision for the leading term. 681573Srgrimes * 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1). 691573Srgrimes * In this case, use a relative precision of 24 bits. 701573Srgrimes * (This is done differently in the original paper) 711573Srgrimes * 721573Srgrimes * Special cases: 731573Srgrimes * 0 return signalling -Inf 741573Srgrimes * neg return signalling NaN 751573Srgrimes * +Inf return +Inf 761573Srgrimes*/ 771573Srgrimes 781573Srgrimes#define N 128 791573Srgrimes 801573Srgrimes/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128. 811573Srgrimes * Used for generation of extend precision logarithms. 821573Srgrimes * The constant 35184372088832 is 2^45, so the divide is exact. 831573Srgrimes * It ensures correct reading of logF_head, even for inaccurate 841573Srgrimes * decimal-to-binary conversion routines. (Everybody gets the 851573Srgrimes * right answer for integers less than 2^53.) 861573Srgrimes * Values for log(F) were generated using error < 10^-57 absolute 871573Srgrimes * with the bc -l package. 881573Srgrimes*/ 891573Srgrimesstatic double A1 = .08333333333333178827; 901573Srgrimesstatic double A2 = .01250000000377174923; 911573Srgrimesstatic double A3 = .002232139987919447809; 921573Srgrimesstatic double A4 = .0004348877777076145742; 931573Srgrimes 941573Srgrimesstatic double logF_head[N+1] = { 951573Srgrimes 0., 961573Srgrimes .007782140442060381246, 971573Srgrimes .015504186535963526694, 981573Srgrimes .023167059281547608406, 991573Srgrimes .030771658666765233647, 1001573Srgrimes .038318864302141264488, 1011573Srgrimes .045809536031242714670, 1021573Srgrimes .053244514518837604555, 1031573Srgrimes .060624621816486978786, 1041573Srgrimes .067950661908525944454, 1051573Srgrimes .075223421237524235039, 1061573Srgrimes .082443669210988446138, 1071573Srgrimes .089612158689760690322, 1081573Srgrimes .096729626458454731618, 1091573Srgrimes .103796793681567578460, 1101573Srgrimes .110814366340264314203, 1111573Srgrimes .117783035656430001836, 1121573Srgrimes .124703478501032805070, 1131573Srgrimes .131576357788617315236, 1141573Srgrimes .138402322859292326029, 1151573Srgrimes .145182009844575077295, 1161573Srgrimes .151916042025732167530, 1171573Srgrimes .158605030176659056451, 1181573Srgrimes .165249572895390883786, 1191573Srgrimes .171850256926518341060, 1201573Srgrimes .178407657472689606947, 1211573Srgrimes .184922338493834104156, 1221573Srgrimes .191394852999565046047, 1231573Srgrimes .197825743329758552135, 1241573Srgrimes .204215541428766300668, 1251573Srgrimes .210564769107350002741, 1261573Srgrimes .216873938300523150246, 1271573Srgrimes .223143551314024080056, 1281573Srgrimes .229374101064877322642, 1291573Srgrimes .235566071312860003672, 1301573Srgrimes .241719936886966024758, 1311573Srgrimes .247836163904594286577, 1321573Srgrimes .253915209980732470285, 1331573Srgrimes .259957524436686071567, 1341573Srgrimes .265963548496984003577, 1351573Srgrimes .271933715484010463114, 1361573Srgrimes .277868451003087102435, 1371573Srgrimes .283768173130738432519, 1381573Srgrimes .289633292582948342896, 1391573Srgrimes .295464212893421063199, 1401573Srgrimes .301261330578199704177, 1411573Srgrimes .307025035294827830512, 1421573Srgrimes .312755710004239517729, 1431573Srgrimes .318453731118097493890, 1441573Srgrimes .324119468654316733591, 1451573Srgrimes .329753286372579168528, 1461573Srgrimes .335355541920762334484, 1471573Srgrimes .340926586970454081892, 1481573Srgrimes .346466767346100823488, 1491573Srgrimes .351976423156884266063, 1501573Srgrimes .357455888922231679316, 1511573Srgrimes .362905493689140712376, 1521573Srgrimes .368325561158599157352, 1531573Srgrimes .373716409793814818840, 1541573Srgrimes .379078352934811846353, 1551573Srgrimes .384411698910298582632, 1561573Srgrimes .389716751140440464951, 1571573Srgrimes .394993808240542421117, 1581573Srgrimes .400243164127459749579, 1591573Srgrimes .405465108107819105498, 1601573Srgrimes .410659924985338875558, 1611573Srgrimes .415827895143593195825, 1621573Srgrimes .420969294644237379543, 1631573Srgrimes .426084395310681429691, 1641573Srgrimes .431173464818130014464, 1651573Srgrimes .436236766774527495726, 1661573Srgrimes .441274560805140936281, 1671573Srgrimes .446287102628048160113, 1681573Srgrimes .451274644139630254358, 1691573Srgrimes .456237433481874177232, 1701573Srgrimes .461175715122408291790, 1711573Srgrimes .466089729924533457960, 1721573Srgrimes .470979715219073113985, 1731573Srgrimes .475845904869856894947, 1741573Srgrimes .480688529345570714212, 1751573Srgrimes .485507815781602403149, 1761573Srgrimes .490303988045525329653, 1771573Srgrimes .495077266798034543171, 1781573Srgrimes .499827869556611403822, 1791573Srgrimes .504556010751912253908, 1801573Srgrimes .509261901790523552335, 1811573Srgrimes .513945751101346104405, 1821573Srgrimes .518607764208354637958, 1831573Srgrimes .523248143765158602036, 1841573Srgrimes .527867089620485785417, 1851573Srgrimes .532464798869114019908, 1861573Srgrimes .537041465897345915436, 1871573Srgrimes .541597282432121573947, 1881573Srgrimes .546132437597407260909, 1891573Srgrimes .550647117952394182793, 1901573Srgrimes .555141507540611200965, 1911573Srgrimes .559615787935399566777, 1921573Srgrimes .564070138285387656651, 1931573Srgrimes .568504735352689749561, 1941573Srgrimes .572919753562018740922, 1951573Srgrimes .577315365035246941260, 1961573Srgrimes .581691739635061821900, 1971573Srgrimes .586049045003164792433, 1981573Srgrimes .590387446602107957005, 1991573Srgrimes .594707107746216934174, 2001573Srgrimes .599008189645246602594, 2011573Srgrimes .603290851438941899687, 2021573Srgrimes .607555250224322662688, 2031573Srgrimes .611801541106615331955, 2041573Srgrimes .616029877215623855590, 2051573Srgrimes .620240409751204424537, 2061573Srgrimes .624433288012369303032, 2071573Srgrimes .628608659422752680256, 2081573Srgrimes .632766669570628437213, 2091573Srgrimes .636907462236194987781, 2101573Srgrimes .641031179420679109171, 2111573Srgrimes .645137961373620782978, 2121573Srgrimes .649227946625615004450, 2131573Srgrimes .653301272011958644725, 2141573Srgrimes .657358072709030238911, 2151573Srgrimes .661398482245203922502, 2161573Srgrimes .665422632544505177065, 2171573Srgrimes .669430653942981734871, 2181573Srgrimes .673422675212350441142, 2191573Srgrimes .677398823590920073911, 2201573Srgrimes .681359224807238206267, 2211573Srgrimes .685304003098281100392, 2221573Srgrimes .689233281238557538017, 2231573Srgrimes .693147180560117703862 2241573Srgrimes}; 2251573Srgrimes 2261573Srgrimesstatic double logF_tail[N+1] = { 2271573Srgrimes 0., 2281573Srgrimes -.00000000000000543229938420049, 2291573Srgrimes .00000000000000172745674997061, 2301573Srgrimes -.00000000000001323017818229233, 2311573Srgrimes -.00000000000001154527628289872, 2321573Srgrimes -.00000000000000466529469958300, 2331573Srgrimes .00000000000005148849572685810, 2341573Srgrimes -.00000000000002532168943117445, 2351573Srgrimes -.00000000000005213620639136504, 2361573Srgrimes -.00000000000001819506003016881, 2371573Srgrimes .00000000000006329065958724544, 2381573Srgrimes .00000000000008614512936087814, 2391573Srgrimes -.00000000000007355770219435028, 2401573Srgrimes .00000000000009638067658552277, 2411573Srgrimes .00000000000007598636597194141, 2421573Srgrimes .00000000000002579999128306990, 2431573Srgrimes -.00000000000004654729747598444, 2441573Srgrimes -.00000000000007556920687451336, 2451573Srgrimes .00000000000010195735223708472, 2461573Srgrimes -.00000000000017319034406422306, 2471573Srgrimes -.00000000000007718001336828098, 2481573Srgrimes .00000000000010980754099855238, 2491573Srgrimes -.00000000000002047235780046195, 2501573Srgrimes -.00000000000008372091099235912, 2511573Srgrimes .00000000000014088127937111135, 2521573Srgrimes .00000000000012869017157588257, 2531573Srgrimes .00000000000017788850778198106, 2541573Srgrimes .00000000000006440856150696891, 2551573Srgrimes .00000000000016132822667240822, 2561573Srgrimes -.00000000000007540916511956188, 2571573Srgrimes -.00000000000000036507188831790, 2581573Srgrimes .00000000000009120937249914984, 2591573Srgrimes .00000000000018567570959796010, 2601573Srgrimes -.00000000000003149265065191483, 2611573Srgrimes -.00000000000009309459495196889, 2621573Srgrimes .00000000000017914338601329117, 2631573Srgrimes -.00000000000001302979717330866, 2641573Srgrimes .00000000000023097385217586939, 2651573Srgrimes .00000000000023999540484211737, 2661573Srgrimes .00000000000015393776174455408, 2671573Srgrimes -.00000000000036870428315837678, 2681573Srgrimes .00000000000036920375082080089, 2691573Srgrimes -.00000000000009383417223663699, 2701573Srgrimes .00000000000009433398189512690, 2711573Srgrimes .00000000000041481318704258568, 2721573Srgrimes -.00000000000003792316480209314, 2731573Srgrimes .00000000000008403156304792424, 2741573Srgrimes -.00000000000034262934348285429, 2751573Srgrimes .00000000000043712191957429145, 2761573Srgrimes -.00000000000010475750058776541, 2771573Srgrimes -.00000000000011118671389559323, 2781573Srgrimes .00000000000037549577257259853, 2791573Srgrimes .00000000000013912841212197565, 2801573Srgrimes .00000000000010775743037572640, 2811573Srgrimes .00000000000029391859187648000, 2821573Srgrimes -.00000000000042790509060060774, 2831573Srgrimes .00000000000022774076114039555, 2841573Srgrimes .00000000000010849569622967912, 2851573Srgrimes -.00000000000023073801945705758, 2861573Srgrimes .00000000000015761203773969435, 2871573Srgrimes .00000000000003345710269544082, 2881573Srgrimes -.00000000000041525158063436123, 2891573Srgrimes .00000000000032655698896907146, 2901573Srgrimes -.00000000000044704265010452446, 2911573Srgrimes .00000000000034527647952039772, 2921573Srgrimes -.00000000000007048962392109746, 2931573Srgrimes .00000000000011776978751369214, 2941573Srgrimes -.00000000000010774341461609578, 2951573Srgrimes .00000000000021863343293215910, 2961573Srgrimes .00000000000024132639491333131, 2971573Srgrimes .00000000000039057462209830700, 2981573Srgrimes -.00000000000026570679203560751, 2991573Srgrimes .00000000000037135141919592021, 3001573Srgrimes -.00000000000017166921336082431, 3011573Srgrimes -.00000000000028658285157914353, 3021573Srgrimes -.00000000000023812542263446809, 3031573Srgrimes .00000000000006576659768580062, 3041573Srgrimes -.00000000000028210143846181267, 3051573Srgrimes .00000000000010701931762114254, 3061573Srgrimes .00000000000018119346366441110, 3071573Srgrimes .00000000000009840465278232627, 3081573Srgrimes -.00000000000033149150282752542, 3091573Srgrimes -.00000000000018302857356041668, 3101573Srgrimes -.00000000000016207400156744949, 3111573Srgrimes .00000000000048303314949553201, 3121573Srgrimes -.00000000000071560553172382115, 3131573Srgrimes .00000000000088821239518571855, 3141573Srgrimes -.00000000000030900580513238244, 3151573Srgrimes -.00000000000061076551972851496, 3161573Srgrimes .00000000000035659969663347830, 3171573Srgrimes .00000000000035782396591276383, 3181573Srgrimes -.00000000000046226087001544578, 3191573Srgrimes .00000000000062279762917225156, 3201573Srgrimes .00000000000072838947272065741, 3211573Srgrimes .00000000000026809646615211673, 3221573Srgrimes -.00000000000010960825046059278, 3231573Srgrimes .00000000000002311949383800537, 3241573Srgrimes -.00000000000058469058005299247, 3251573Srgrimes -.00000000000002103748251144494, 3261573Srgrimes -.00000000000023323182945587408, 3271573Srgrimes -.00000000000042333694288141916, 3281573Srgrimes -.00000000000043933937969737844, 3291573Srgrimes .00000000000041341647073835565, 3301573Srgrimes .00000000000006841763641591466, 3311573Srgrimes .00000000000047585534004430641, 3321573Srgrimes .00000000000083679678674757695, 3331573Srgrimes -.00000000000085763734646658640, 3341573Srgrimes .00000000000021913281229340092, 3351573Srgrimes -.00000000000062242842536431148, 3361573Srgrimes -.00000000000010983594325438430, 3371573Srgrimes .00000000000065310431377633651, 3381573Srgrimes -.00000000000047580199021710769, 3391573Srgrimes -.00000000000037854251265457040, 3401573Srgrimes .00000000000040939233218678664, 3411573Srgrimes .00000000000087424383914858291, 3421573Srgrimes .00000000000025218188456842882, 3431573Srgrimes -.00000000000003608131360422557, 3441573Srgrimes -.00000000000050518555924280902, 3451573Srgrimes .00000000000078699403323355317, 3461573Srgrimes -.00000000000067020876961949060, 3471573Srgrimes .00000000000016108575753932458, 3481573Srgrimes .00000000000058527188436251509, 3491573Srgrimes -.00000000000035246757297904791, 3501573Srgrimes -.00000000000018372084495629058, 3511573Srgrimes .00000000000088606689813494916, 3521573Srgrimes .00000000000066486268071468700, 3531573Srgrimes .00000000000063831615170646519, 3541573Srgrimes .00000000000025144230728376072, 3551573Srgrimes -.00000000000017239444525614834 3561573Srgrimes}; 3571573Srgrimes 35893211Sbde#if 0 3591573Srgrimesdouble 3601573Srgrimes#ifdef _ANSI_SOURCE 3611573Srgrimeslog(double x) 3621573Srgrimes#else 3631573Srgrimeslog(x) double x; 3641573Srgrimes#endif 3651573Srgrimes{ 3661573Srgrimes int m, j; 3671573Srgrimes double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0; 3681573Srgrimes volatile double u1; 3691573Srgrimes 3701573Srgrimes /* Catch special cases */ 3711573Srgrimes if (x <= 0) 372138924Sdas if (x == zero) /* log(0) = -Inf */ 3731573Srgrimes return (-one/zero); 374138924Sdas else /* log(neg) = NaN */ 3751573Srgrimes return (zero/zero); 3761573Srgrimes else if (!finite(x)) 377138924Sdas return (x+x); /* x = NaN, Inf */ 3788870Srgrimes 3791573Srgrimes /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 3801573Srgrimes /* y = F*(1 + f/F) for |f| <= 2^-8 */ 3811573Srgrimes 3821573Srgrimes m = logb(x); 3831573Srgrimes g = ldexp(x, -m); 384138924Sdas if (m == -1022) { 3851573Srgrimes j = logb(g), m += j; 3861573Srgrimes g = ldexp(g, -j); 3871573Srgrimes } 3881573Srgrimes j = N*(g-1) + .5; 3891573Srgrimes F = (1.0/N) * j + 1; /* F*128 is an integer in [128, 512] */ 3901573Srgrimes f = g - F; 3911573Srgrimes 3921573Srgrimes /* Approximate expansion for log(1+f/F) ~= u + q */ 3931573Srgrimes g = 1/(2*F+f); 3941573Srgrimes u = 2*f*g; 3951573Srgrimes v = u*u; 3961573Srgrimes q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 3971573Srgrimes 3981573Srgrimes /* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8, 3991573Srgrimes * u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits. 4001573Srgrimes * It also adds exactly to |m*log2_hi + log_F_head[j] | < 750 4011573Srgrimes */ 4021573Srgrimes if (m | j) 4031573Srgrimes u1 = u + 513, u1 -= 513; 4041573Srgrimes 4051573Srgrimes /* case 2: |1-x| < 1/256. The m- and j- dependent terms are zero; 4061573Srgrimes * u1 = u to 24 bits. 4071573Srgrimes */ 4081573Srgrimes else 4091573Srgrimes u1 = u, TRUNC(u1); 4101573Srgrimes u2 = (2.0*(f - F*u1) - u1*f) * g; 4111573Srgrimes /* u1 + u2 = 2f/(2F+f) to extra precision. */ 4121573Srgrimes 4131573Srgrimes /* log(x) = log(2^m*F*(1+f/F)) = */ 4141573Srgrimes /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); */ 4151573Srgrimes /* (exact) + (tiny) */ 4161573Srgrimes 4171573Srgrimes u1 += m*logF_head[N] + logF_head[j]; /* exact */ 4181573Srgrimes u2 = (u2 + logF_tail[j]) + q; /* tiny */ 4191573Srgrimes u2 += logF_tail[N]*m; 4201573Srgrimes return (u1 + u2); 4211573Srgrimes} 42293211Sbde#endif 4231573Srgrimes 4241573Srgrimes/* 4251573Srgrimes * Extra precision variant, returning struct {double a, b;}; 426108533Sschweikh * log(x) = a+b to 63 bits, with a rounded to 26 bits. 4271573Srgrimes */ 4281573Srgrimesstruct Double 4291573Srgrimes#ifdef _ANSI_SOURCE 4301573Srgrimes__log__D(double x) 4311573Srgrimes#else 4321573Srgrimes__log__D(x) double x; 4331573Srgrimes#endif 4341573Srgrimes{ 4351573Srgrimes int m, j; 436150318Sbde double F, f, g, q, u, v, u2; 4371573Srgrimes volatile double u1; 4381573Srgrimes struct Double r; 4391573Srgrimes 4401573Srgrimes /* Argument reduction: 1 <= g < 2; x/2^m = g; */ 4411573Srgrimes /* y = F*(1 + f/F) for |f| <= 2^-8 */ 4421573Srgrimes 4431573Srgrimes m = logb(x); 4441573Srgrimes g = ldexp(x, -m); 445138924Sdas if (m == -1022) { 4461573Srgrimes j = logb(g), m += j; 4471573Srgrimes g = ldexp(g, -j); 4481573Srgrimes } 4491573Srgrimes j = N*(g-1) + .5; 4501573Srgrimes F = (1.0/N) * j + 1; 4511573Srgrimes f = g - F; 4521573Srgrimes 4531573Srgrimes g = 1/(2*F+f); 4541573Srgrimes u = 2*f*g; 4551573Srgrimes v = u*u; 4561573Srgrimes q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); 4571573Srgrimes if (m | j) 4581573Srgrimes u1 = u + 513, u1 -= 513; 4591573Srgrimes else 4601573Srgrimes u1 = u, TRUNC(u1); 4611573Srgrimes u2 = (2.0*(f - F*u1) - u1*f) * g; 4621573Srgrimes 4631573Srgrimes u1 += m*logF_head[N] + logF_head[j]; 4641573Srgrimes 4651573Srgrimes u2 += logF_tail[j]; u2 += q; 4661573Srgrimes u2 += logF_tail[N]*m; 4671573Srgrimes r.a = u1 + u2; /* Only difference is here */ 4681573Srgrimes TRUNC(r.a); 4691573Srgrimes r.b = (u1 - r.a) + u2; 4701573Srgrimes return (r); 4711573Srgrimes} 472