Mercurial > hg > CbC > CbC_gcc
annotate libdecnumber/decBasic.c @ 145:1830386684a0
gcc-9.2.0
| author | anatofuz |
|---|---|
| date | Thu, 13 Feb 2020 11:34:05 +0900 |
| parents | 84e7813d76e9 |
| children |
| rev | line source |
|---|---|
| 0 | 1 /* Common base code for the decNumber C Library. |
| 145 | 2 Copyright (C) 2007-2020 Free Software Foundation, Inc. |
| 0 | 3 Contributed by IBM Corporation. Author Mike Cowlishaw. |
| 4 | |
| 5 This file is part of GCC. | |
| 6 | |
| 7 GCC is free software; you can redistribute it and/or modify it under | |
| 8 the terms of the GNU General Public License as published by the Free | |
| 9 Software Foundation; either version 3, or (at your option) any later | |
| 10 version. | |
| 11 | |
| 12 GCC is distributed in the hope that it will be useful, but WITHOUT ANY | |
| 13 WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 15 for more details. | |
| 16 | |
| 17 Under Section 7 of GPL version 3, you are granted additional | |
| 18 permissions described in the GCC Runtime Library Exception, version | |
| 19 3.1, as published by the Free Software Foundation. | |
| 20 | |
| 21 You should have received a copy of the GNU General Public License and | |
| 22 a copy of the GCC Runtime Library Exception along with this program; | |
| 23 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see | |
| 24 <http://www.gnu.org/licenses/>. */ | |
| 25 | |
| 26 /* ------------------------------------------------------------------ */ | |
| 27 /* decBasic.c -- common base code for Basic decimal types */ | |
| 28 /* ------------------------------------------------------------------ */ | |
| 29 /* This module comprises code that is shared between decDouble and */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
30 /* decQuad (but not decSingle). The main arithmetic operations are */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
31 /* here (Add, Subtract, Multiply, FMA, and Division operators). */ |
| 0 | 32 /* */ |
| 33 /* Unlike decNumber, parameterization takes place at compile time */ | |
| 34 /* rather than at runtime. The parameters are set in the decDouble.c */ | |
| 35 /* (etc.) files, which then include this one to produce the compiled */ | |
| 36 /* code. The functions here, therefore, are code shared between */ | |
| 37 /* multiple formats. */ | |
| 38 /* */ | |
| 39 /* This must be included after decCommon.c. */ | |
| 40 /* ------------------------------------------------------------------ */ | |
| 41 /* Names here refer to decFloat rather than to decDouble, etc., and */ | |
| 42 /* the functions are in strict alphabetical order. */ | |
| 43 | |
| 44 /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */ | |
| 45 /* decCommon.c */ | |
| 46 #if !defined(QUAD) | |
| 47 #error decBasic.c must be included after decCommon.c | |
| 48 #endif | |
| 49 #if SINGLE | |
| 50 #error Routines in decBasic.c are for decDouble and decQuad only | |
| 51 #endif | |
| 52 | |
| 53 /* Private constants */ | |
| 54 #define DIVIDE 0x80000000 /* Divide operations [as flags] */ | |
| 55 #define REMAINDER 0x40000000 /* .. */ | |
| 56 #define DIVIDEINT 0x20000000 /* .. */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
57 #define REMNEAR 0x10000000 /* .. */ |
| 0 | 58 |
| 59 /* Private functions (local, used only by routines in this module) */ | |
| 60 static decFloat *decDivide(decFloat *, const decFloat *, | |
| 61 const decFloat *, decContext *, uInt); | |
| 62 static decFloat *decCanonical(decFloat *, const decFloat *); | |
| 63 static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *, | |
| 64 const decFloat *); | |
| 65 static decFloat *decInfinity(decFloat *, const decFloat *); | |
| 66 static decFloat *decInvalid(decFloat *, decContext *); | |
| 67 static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *, | |
| 68 decContext *); | |
| 69 static Int decNumCompare(const decFloat *, const decFloat *, Flag); | |
| 70 static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *, | |
| 71 enum rounding, Flag); | |
| 72 static uInt decToInt32(const decFloat *, decContext *, enum rounding, | |
| 73 Flag, Flag); | |
| 74 | |
| 75 /* ------------------------------------------------------------------ */ | |
| 76 /* decCanonical -- copy a decFloat, making canonical */ | |
| 77 /* */ | |
| 78 /* result gets the canonicalized df */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
79 /* df is the decFloat to copy and make canonical */ |
| 0 | 80 /* returns result */ |
| 81 /* */ | |
| 82 /* This is exposed via decFloatCanonical for Double and Quad only. */ | |
| 83 /* This works on specials, too; no error or exception is possible. */ | |
| 84 /* ------------------------------------------------------------------ */ | |
| 85 static decFloat * decCanonical(decFloat *result, const decFloat *df) { | |
| 86 uInt encode, precode, dpd; /* work */ | |
| 87 uInt inword, uoff, canon; /* .. */ | |
| 88 Int n; /* counter (down) */ | |
| 89 if (df!=result) *result=*df; /* effect copy if needed */ | |
| 90 if (DFISSPECIAL(result)) { | |
| 91 if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */ | |
| 92 /* is a NaN */ | |
| 93 DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */ | |
| 94 if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */ | |
| 95 /* drop through to check payload */ | |
| 96 } | |
| 97 /* return quickly if the coefficient continuation is canonical */ | |
| 98 { /* declare block */ | |
| 99 #if DOUBLE | |
| 100 uInt sourhi=DFWORD(df, 0); | |
| 101 uInt sourlo=DFWORD(df, 1); | |
| 102 if (CANONDPDOFF(sourhi, 8) | |
| 103 && CANONDPDTWO(sourhi, sourlo, 30) | |
| 104 && CANONDPDOFF(sourlo, 20) | |
| 105 && CANONDPDOFF(sourlo, 10) | |
| 106 && CANONDPDOFF(sourlo, 0)) return result; | |
| 107 #elif QUAD | |
| 108 uInt sourhi=DFWORD(df, 0); | |
| 109 uInt sourmh=DFWORD(df, 1); | |
| 110 uInt sourml=DFWORD(df, 2); | |
| 111 uInt sourlo=DFWORD(df, 3); | |
| 112 if (CANONDPDOFF(sourhi, 4) | |
| 113 && CANONDPDTWO(sourhi, sourmh, 26) | |
| 114 && CANONDPDOFF(sourmh, 16) | |
| 115 && CANONDPDOFF(sourmh, 6) | |
| 116 && CANONDPDTWO(sourmh, sourml, 28) | |
| 117 && CANONDPDOFF(sourml, 18) | |
| 118 && CANONDPDOFF(sourml, 8) | |
| 119 && CANONDPDTWO(sourml, sourlo, 30) | |
| 120 && CANONDPDOFF(sourlo, 20) | |
| 121 && CANONDPDOFF(sourlo, 10) | |
| 122 && CANONDPDOFF(sourlo, 0)) return result; | |
| 123 #endif | |
| 124 } /* block */ | |
| 125 | |
| 126 /* Loop to repair a non-canonical coefficent, as needed */ | |
| 127 inword=DECWORDS-1; /* current input word */ | |
| 128 uoff=0; /* bit offset of declet */ | |
| 129 encode=DFWORD(result, inword); | |
| 130 for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */ | |
| 131 dpd=encode>>uoff; | |
| 132 uoff+=10; | |
| 133 if (uoff>32) { /* crossed uInt boundary */ | |
| 134 inword--; | |
| 135 encode=DFWORD(result, inword); | |
| 136 uoff-=32; | |
| 137 dpd|=encode<<(10-uoff); /* get pending bits */ | |
| 138 } | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
139 dpd&=0x3ff; /* clear uninteresting bits */ |
| 0 | 140 if (dpd<0x16e) continue; /* must be canonical */ |
| 141 canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */ | |
| 142 if (canon==dpd) continue; /* have canonical declet */ | |
| 143 /* need to replace declet */ | |
| 144 if (uoff>=10) { /* all within current word */ | |
| 145 encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
146 encode|=canon<<(uoff-10); /* insert the canonical form */ |
| 0 | 147 DFWORD(result, inword)=encode; /* .. and save */ |
| 148 continue; | |
| 149 } | |
| 150 /* straddled words */ | |
| 151 precode=DFWORD(result, inword+1); /* get previous */ | |
| 152 precode&=0xffffffff>>(10-uoff); /* clear top bits */ | |
| 153 DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff))); | |
| 154 encode&=0xffffffff<<uoff; /* clear bottom bits */ | |
| 155 encode|=canon>>(10-uoff); /* insert canonical */ | |
| 156 DFWORD(result, inword)=encode; /* .. and save */ | |
| 157 } /* n */ | |
| 158 return result; | |
| 159 } /* decCanonical */ | |
| 160 | |
| 161 /* ------------------------------------------------------------------ */ | |
| 162 /* decDivide -- divide operations */ | |
| 163 /* */ | |
| 164 /* result gets the result of dividing dfl by dfr: */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
165 /* dfl is the first decFloat (lhs) */ |
| 0 | 166 /* dfr is the second decFloat (rhs) */ |
| 167 /* set is the context */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
168 /* op is the operation selector */ |
| 0 | 169 /* returns result */ |
| 170 /* */ | |
| 171 /* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */ | |
| 172 /* ------------------------------------------------------------------ */ | |
| 173 #define DIVCOUNT 0 /* 1 to instrument subtractions counter */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
174 #define DIVBASE ((uInt)BILLION) /* the base used for divide */ |
| 0 | 175 #define DIVOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */ |
| 176 #define DIVACCLEN (DIVOPLEN*3) /* accumulator length (ditto) */ | |
| 177 static decFloat * decDivide(decFloat *result, const decFloat *dfl, | |
| 178 const decFloat *dfr, decContext *set, uInt op) { | |
| 179 decFloat quotient; /* for remainders */ | |
| 180 bcdnum num; /* for final conversion */ | |
| 181 uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
182 uInt div[DIVOPLEN]; /* divisor in base-billion .. */ |
| 0 | 183 uInt quo[DIVOPLEN+1]; /* quotient in base-billion .. */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
184 uByte bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */ |
| 0 | 185 uInt *msua, *msud, *msuq; /* -> msu of acc, div, and quo */ |
| 186 Int divunits, accunits; /* lengths */ | |
| 187 Int quodigits; /* digits in quotient */ | |
| 188 uInt *lsua, *lsuq; /* -> current acc and quo lsus */ | |
| 189 Int length, multiplier; /* work */ | |
| 190 uInt carry, sign; /* .. */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
191 uInt *ua, *ud, *uq; /* .. */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
192 uByte *ub; /* .. */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
193 uInt uiwork; /* for macros */ |
| 0 | 194 uInt divtop; /* top unit of div adjusted for estimating */ |
| 195 #if DIVCOUNT | |
| 196 static uInt maxcount=0; /* worst-seen subtractions count */ | |
| 197 uInt divcount=0; /* subtractions count [this divide] */ | |
| 198 #endif | |
| 199 | |
| 200 /* calculate sign */ | |
| 201 num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign; | |
| 202 | |
| 203 if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */ | |
| 204 /* NaNs are handled as usual */ | |
| 205 if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set); | |
| 206 /* one or two infinities */ | |
| 207 if (DFISINF(dfl)) { | |
| 208 if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */ | |
| 209 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */ | |
| 210 /* Infinity/x is infinite and quiet, even if x=0 */ | |
| 211 DFWORD(result, 0)=num.sign; | |
| 212 return decInfinity(result, result); | |
| 213 } | |
| 214 /* must be x/Infinity -- remainders are lhs */ | |
| 215 if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl); | |
| 216 /* divides: return zero with correct sign and exponent depending */ | |
| 217 /* on op (Etiny for divide, 0 for divideInt) */ | |
| 218 decFloatZero(result); | |
| 219 if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */ | |
| 220 else DFWORD(result, 0)=num.sign; /* zeros the exponent, too */ | |
| 221 return result; | |
| 222 } | |
| 223 /* next, handle zero operands (x/0 and 0/x) */ | |
| 224 if (DFISZERO(dfr)) { /* x/0 */ | |
| 225 if (DFISZERO(dfl)) { /* 0/0 is undefined */ | |
| 226 decFloatZero(result); | |
| 227 DFWORD(result, 0)=DECFLOAT_qNaN; | |
| 228 set->status|=DEC_Division_undefined; | |
| 229 return result; | |
| 230 } | |
| 231 if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */ | |
| 232 set->status|=DEC_Division_by_zero; | |
| 233 DFWORD(result, 0)=num.sign; | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
234 return decInfinity(result, result); /* x/0 -> signed Infinity */ |
| 0 | 235 } |
| 236 num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */ | |
| 237 if (DFISZERO(dfl)) { /* 0/x (x!=0) */ | |
| 238 /* if divide, result is 0 with ideal exponent; divideInt has */ | |
| 239 /* exponent=0, remainders give zero with lower exponent */ | |
| 240 if (op&DIVIDEINT) { | |
| 241 decFloatZero(result); | |
| 242 DFWORD(result, 0)|=num.sign; /* add sign */ | |
| 243 return result; | |
| 244 } | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
245 if (!(op&DIVIDE)) { /* a remainder */ |
| 0 | 246 /* exponent is the minimum of the operands */ |
| 247 num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr)); | |
| 248 /* if the result is zero the sign shall be sign of dfl */ | |
| 249 num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign; | |
| 250 } | |
| 251 bcdacc[0]=0; | |
| 252 num.msd=bcdacc; /* -> 0 */ | |
| 253 num.lsd=bcdacc; /* .. */ | |
| 254 return decFinalize(result, &num, set); /* [divide may clamp exponent] */ | |
| 255 } /* 0/x */ | |
| 256 /* [here, both operands are known to be finite and non-zero] */ | |
| 257 | |
| 258 /* extract the operand coefficents into 'units' which are */ | |
| 259 /* base-billion; the lhs is high-aligned in acc and the msu of both */ | |
| 260 /* acc and div is at the right-hand end of array (offset length-1); */ | |
| 261 /* the quotient can need one more unit than the operands as digits */ | |
| 262 /* in it are not necessarily aligned neatly; further, the quotient */ | |
| 263 /* may not start accumulating until after the end of the initial */ | |
| 264 /* operand in acc if that is small (e.g., 1) so the accumulator */ | |
| 265 /* must have at least that number of units extra (at the ls end) */ | |
| 266 GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN); | |
| 267 GETCOEFFBILL(dfr, div); | |
| 268 /* zero the low uInts of acc */ | |
| 269 acc[0]=0; | |
| 270 acc[1]=0; | |
| 271 acc[2]=0; | |
| 272 acc[3]=0; | |
| 273 #if DOUBLE | |
| 274 #if DIVOPLEN!=2 | |
| 275 #error Unexpected Double DIVOPLEN | |
| 276 #endif | |
| 277 #elif QUAD | |
| 278 acc[4]=0; | |
| 279 acc[5]=0; | |
| 280 acc[6]=0; | |
| 281 acc[7]=0; | |
| 282 #if DIVOPLEN!=4 | |
| 283 #error Unexpected Quad DIVOPLEN | |
| 284 #endif | |
| 285 #endif | |
| 286 | |
| 287 /* set msu and lsu pointers */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
288 msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */ |
| 0 | 289 msuq=quo+DIVOPLEN; |
| 290 /*[loop for div will terminate because operands are non-zero] */ | |
| 291 for (msud=div+DIVOPLEN-1; *msud==0;) msud--; | |
| 292 /* the initial least-significant unit of acc is set so acc appears */ | |
| 293 /* to have the same length as div. */ | |
| 294 /* This moves one position towards the least possible for each */ | |
| 295 /* iteration */ | |
| 296 divunits=(Int)(msud-div+1); /* precalculate */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
297 lsua=msua-divunits+1; /* initial working lsu of acc */ |
| 0 | 298 lsuq=msuq; /* and of quo */ |
| 299 | |
| 300 /* set up the estimator for the multiplier; this is the msu of div, */ | |
| 301 /* plus two bits from the unit below (if any) rounded up by one if */ | |
| 302 /* there are any non-zero bits or units below that [the extra two */ | |
| 303 /* bits makes for a much better estimate when the top unit is small] */ | |
| 304 divtop=*msud<<2; | |
| 305 if (divunits>1) { | |
| 306 uInt *um=msud-1; | |
| 307 uInt d=*um; | |
| 308 if (d>=750000000) {divtop+=3; d-=750000000;} | |
| 309 else if (d>=500000000) {divtop+=2; d-=500000000;} | |
| 310 else if (d>=250000000) {divtop++; d-=250000000;} | |
| 311 if (d) divtop++; | |
| 312 else for (um--; um>=div; um--) if (*um) { | |
| 313 divtop++; | |
| 314 break; | |
| 315 } | |
| 316 } /* >1 unit */ | |
| 317 | |
| 318 #if DECTRACE | |
| 319 {Int i; | |
| 320 printf("----- div="); | |
| 321 for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]); | |
| 322 printf("\n");} | |
| 323 #endif | |
| 324 | |
| 325 /* now collect up to DECPMAX+1 digits in the quotient (this may */ | |
| 326 /* need OPLEN+1 uInts if unaligned) */ | |
| 327 quodigits=0; /* no digits yet */ | |
| 328 for (;; lsua--) { /* outer loop -- each input position */ | |
| 329 #if DECCHECK | |
| 330 if (lsua<acc) { | |
| 331 printf("Acc underrun...\n"); | |
| 332 break; | |
| 333 } | |
| 334 #endif | |
| 335 #if DECTRACE | |
| 336 printf("Outer: quodigits=%ld acc=", (LI)quodigits); | |
| 337 for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua); | |
| 338 printf("\n"); | |
| 339 #endif | |
| 340 *lsuq=0; /* default unit result is 0 */ | |
| 341 for (;;) { /* inner loop -- calculate quotient unit */ | |
| 342 /* strip leading zero units from acc (either there initially or */ | |
| 343 /* from subtraction below); this may strip all if exactly 0 */ | |
| 344 for (; *msua==0 && msua>=lsua;) msua--; | |
| 345 accunits=(Int)(msua-lsua+1); /* [maybe 0] */ | |
| 346 /* subtraction is only necessary and possible if there are as */ | |
| 347 /* least as many units remaining in acc for this iteration as */ | |
| 348 /* there are in div */ | |
| 349 if (accunits<divunits) { | |
| 350 if (accunits==0) msua++; /* restore */ | |
| 351 break; | |
| 352 } | |
| 353 | |
| 354 /* If acc is longer than div then subtraction is definitely */ | |
| 355 /* possible (as msu of both is non-zero), but if they are the */ | |
| 356 /* same length a comparison is needed. */ | |
| 357 /* If a subtraction is needed then a good estimate of the */ | |
| 358 /* multiplier for the subtraction is also needed in order to */ | |
| 359 /* minimise the iterations of this inner loop because the */ | |
| 360 /* subtractions needed dominate division performance. */ | |
| 361 if (accunits==divunits) { | |
| 362 /* compare the high divunits of acc and div: */ | |
| 363 /* acc<div: this quotient unit is unchanged; subtraction */ | |
| 364 /* will be possible on the next iteration */ | |
| 365 /* acc==div: quotient gains 1, set acc=0 */ | |
| 366 /* acc>div: subtraction necessary at this position */ | |
| 367 for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break; | |
| 368 /* [now at first mismatch or lsu] */ | |
| 369 if (*ud>*ua) break; /* next time... */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
370 if (*ud==*ua) { /* all compared equal */ |
| 0 | 371 *lsuq+=1; /* increment result */ |
| 372 msua=lsua; /* collapse acc units */ | |
| 373 *msua=0; /* .. to a zero */ | |
| 374 break; | |
| 375 } | |
| 376 | |
| 377 /* subtraction necessary; estimate multiplier [see above] */ | |
| 378 /* if both *msud and *msua are small it is cost-effective to */ | |
| 379 /* bring in part of the following units (if any) to get a */ | |
| 380 /* better estimate (assume some other non-zero in div) */ | |
| 381 #define DIVLO 1000000U | |
| 382 #define DIVHI (DIVBASE/DIVLO) | |
| 383 #if DECUSE64 | |
| 384 if (divunits>1) { | |
| 385 /* there cannot be a *(msud-2) for DECDOUBLE so next is */ | |
| 386 /* an exact calculation unless DECQUAD (which needs to */ | |
| 387 /* assume bits out there if divunits>2) */ | |
| 388 uLong mul=(uLong)*msua * DIVBASE + *(msua-1); | |
| 389 uLong div=(uLong)*msud * DIVBASE + *(msud-1); | |
| 390 #if QUAD | |
| 391 if (divunits>2) div++; | |
| 392 #endif | |
| 393 mul/=div; | |
| 394 multiplier=(Int)mul; | |
| 395 } | |
| 396 else multiplier=*msua/(*msud); | |
| 397 #else | |
| 398 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) { | |
| 399 multiplier=(*msua*DIVHI + *(msua-1)/DIVLO) | |
| 400 /(*msud*DIVHI + *(msud-1)/DIVLO +1); | |
| 401 } | |
| 402 else multiplier=(*msua<<2)/divtop; | |
| 403 #endif | |
| 404 } | |
| 405 else { /* accunits>divunits */ | |
| 406 /* msud is one unit 'lower' than msua, so estimate differently */ | |
| 407 #if DECUSE64 | |
| 408 uLong mul; | |
| 409 /* as before, bring in extra digits if possible */ | |
| 410 if (divunits>1 && *msua<DIVLO && *msud<DIVLO) { | |
| 411 mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI | |
| 412 + *(msua-2)/DIVLO; | |
| 413 mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1); | |
| 414 } | |
| 415 else if (divunits==1) { | |
| 416 mul=(uLong)*msua * DIVBASE + *(msua-1); | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
417 mul/=*msud; /* no more to the right */ |
| 0 | 418 } |
| 419 else { | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
420 mul=(uLong)(*msua) * (uInt)(DIVBASE<<2) |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
421 + (*(msua-1)<<2); |
| 0 | 422 mul/=divtop; /* [divtop already allows for sticky bits] */ |
| 423 } | |
| 424 multiplier=(Int)mul; | |
| 425 #else | |
| 426 multiplier=*msua * ((DIVBASE<<2)/divtop); | |
| 427 #endif | |
| 428 } | |
| 429 if (multiplier==0) multiplier=1; /* marginal case */ | |
| 430 *lsuq+=multiplier; | |
| 431 | |
| 432 #if DIVCOUNT | |
| 433 /* printf("Multiplier: %ld\n", (LI)multiplier); */ | |
| 434 divcount++; | |
| 435 #endif | |
| 436 | |
| 437 /* Carry out the subtraction acc-(div*multiplier); for each */ | |
| 438 /* unit in div, do the multiply, split to units (see */ | |
| 439 /* decFloatMultiply for the algorithm), and subtract from acc */ | |
| 440 #define DIVMAGIC 2305843009U /* 2**61/10**9 */ | |
| 441 #define DIVSHIFTA 29 | |
| 442 #define DIVSHIFTB 32 | |
| 443 carry=0; | |
| 444 for (ud=div, ua=lsua; ud<=msud; ud++, ua++) { | |
| 445 uInt lo, hop; | |
| 446 #if DECUSE64 | |
| 447 uLong sub=(uLong)multiplier*(*ud)+carry; | |
| 448 if (sub<DIVBASE) { | |
| 449 carry=0; | |
| 450 lo=(uInt)sub; | |
| 451 } | |
| 452 else { | |
| 453 hop=(uInt)(sub>>DIVSHIFTA); | |
| 454 carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB); | |
| 455 /* the estimate is now in hi; now calculate sub-hi*10**9 */ | |
| 456 /* to get the remainder (which will be <DIVBASE)) */ | |
| 457 lo=(uInt)sub; | |
| 458 lo-=carry*DIVBASE; /* low word of result */ | |
| 459 if (lo>=DIVBASE) { | |
| 460 lo-=DIVBASE; /* correct by +1 */ | |
| 461 carry++; | |
| 462 } | |
| 463 } | |
| 464 #else /* 32-bit */ | |
| 465 uInt hi; | |
| 466 /* calculate multiplier*(*ud) into hi and lo */ | |
| 467 LONGMUL32HI(hi, *ud, multiplier); /* get the high word */ | |
| 468 lo=multiplier*(*ud); /* .. and the low */ | |
| 469 lo+=carry; /* add the old hi */ | |
| 470 carry=hi+(lo<carry); /* .. with any carry */ | |
| 471 if (carry || lo>=DIVBASE) { /* split is needed */ | |
| 472 hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */ | |
| 473 LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */ | |
| 474 /* [DIVSHIFTB is 32, so carry can be used directly] */ | |
| 475 /* the estimate is now in carry; now calculate hi:lo-est*10**9; */ | |
| 476 /* happily the top word of the result is irrelevant because it */ | |
| 477 /* will always be zero so this needs only one multiplication */ | |
| 478 lo-=(carry*DIVBASE); | |
| 479 /* the correction here will be at most +1; do it */ | |
| 480 if (lo>=DIVBASE) { | |
| 481 lo-=DIVBASE; | |
| 482 carry++; | |
| 483 } | |
| 484 } | |
| 485 #endif | |
| 486 if (lo>*ua) { /* borrow needed */ | |
| 487 *ua+=DIVBASE; | |
| 488 carry++; | |
| 489 } | |
| 490 *ua-=lo; | |
| 491 } /* ud loop */ | |
| 492 if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */ | |
| 493 } /* inner loop */ | |
| 494 | |
| 495 /* the outer loop terminates when there is either an exact result */ | |
| 496 /* or enough digits; first update the quotient digit count and */ | |
| 497 /* pointer (if any significant digits) */ | |
| 498 #if DECTRACE | |
| 499 if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq); | |
| 500 #endif | |
| 501 if (quodigits) { | |
| 502 quodigits+=9; /* had leading unit earlier */ | |
| 503 lsuq--; | |
| 504 if (quodigits>DECPMAX+1) break; /* have enough */ | |
| 505 } | |
| 506 else if (*lsuq) { /* first quotient digits */ | |
| 507 const uInt *pow; | |
| 508 for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++; | |
| 509 lsuq--; | |
| 510 /* [cannot have >DECPMAX+1 on first unit] */ | |
| 511 } | |
| 512 | |
| 513 if (*msua!=0) continue; /* not an exact result */ | |
| 514 /* acc is zero iff used all of original units and zero down to lsua */ | |
| 515 /* (must also continue to original lsu for correct quotient length) */ | |
| 516 if (lsua>acc+DIVACCLEN-DIVOPLEN) continue; | |
| 517 for (; msua>lsua && *msua==0;) msua--; | |
| 518 if (*msua==0 && msua==lsua) break; | |
| 519 } /* outer loop */ | |
| 520 | |
| 521 /* all of the original operand in acc has been covered at this point */ | |
| 522 /* quotient now has at least DECPMAX+2 digits */ | |
| 523 /* *msua is now non-0 if inexact and sticky bits */ | |
| 524 /* lsuq is one below the last uint of the quotient */ | |
| 525 lsuq++; /* set -> true lsu of quo */ | |
| 526 if (*msua) *lsuq|=1; /* apply sticky bit */ | |
| 527 | |
| 528 /* quo now holds the (unrounded) quotient in base-billion; one */ | |
| 529 /* base-billion 'digit' per uInt. */ | |
| 530 #if DECTRACE | |
| 531 printf("DivQuo:"); | |
| 532 for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq); | |
| 533 printf("\n"); | |
| 534 #endif | |
| 535 | |
| 536 /* Now convert to BCD for rounding and cleanup, starting from the */ | |
| 537 /* most significant end [offset by one into bcdacc to leave room */ | |
| 538 /* for a possible carry digit if rounding for REMNEAR is needed] */ | |
| 539 for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) { | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
540 uInt top, mid, rem; /* work */ |
| 0 | 541 if (*uq==0) { /* no split needed */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
542 UBFROMUI(ub, 0); /* clear 9 BCD8s */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
543 UBFROMUI(ub+4, 0); /* .. */ |
| 0 | 544 *(ub+8)=0; /* .. */ |
| 545 continue; | |
| 546 } | |
| 547 /* *uq is non-zero -- split the base-billion digit into */ | |
| 548 /* hi, mid, and low three-digits */ | |
| 549 #define divsplit9 1000000 /* divisor */ | |
| 550 #define divsplit6 1000 /* divisor */ | |
| 551 /* The splitting is done by simple divides and remainders, */ | |
| 552 /* assuming the compiler will optimize these [GCC does] */ | |
| 553 top=*uq/divsplit9; | |
| 554 rem=*uq%divsplit9; | |
| 555 mid=rem/divsplit6; | |
| 556 rem=rem%divsplit6; | |
| 557 /* lay out the nine BCD digits (plus one unwanted byte) */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
558 UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4])); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
559 UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4])); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
560 UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4])); |
| 0 | 561 } /* BCD conversion loop */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
562 ub--; /* -> lsu */ |
| 0 | 563 |
| 564 /* complete the bcdnum; quodigits is correct, so the position of */ | |
| 565 /* the first non-zero is known */ | |
| 566 num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits; | |
| 567 num.lsd=ub; | |
| 568 | |
| 569 /* make exponent adjustments, etc */ | |
| 570 if (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */ | |
| 571 num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9); | |
| 572 /* if the result was exact then there may be up to 8 extra */ | |
| 573 /* trailing zeros in the overflowed quotient final unit */ | |
| 574 if (*msua==0) { | |
| 575 for (; *ub==0;) ub--; /* drop zeros */ | |
| 576 num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */ | |
| 577 num.lsd=ub; | |
| 578 } | |
| 579 } /* adjustment needed */ | |
| 580 | |
| 581 #if DIVCOUNT | |
| 582 if (divcount>maxcount) { /* new high-water nark */ | |
| 583 maxcount=divcount; | |
| 584 printf("DivNewMaxCount: %ld\n", (LI)maxcount); | |
| 585 } | |
| 586 #endif | |
| 587 | |
| 588 if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */ | |
| 589 | |
| 590 /* Is DIVIDEINT or a remainder; there is more to do -- first form */ | |
| 591 /* the integer (this is done 'after the fact', unlike as in */ | |
| 592 /* decNumber, so as not to tax DIVIDE) */ | |
| 593 | |
| 594 /* The first non-zero digit will be in the first 9 digits, known */ | |
| 595 /* from quodigits and num.msd, so there is always space for DECPMAX */ | |
| 596 /* digits */ | |
| 597 | |
| 598 length=(Int)(num.lsd-num.msd+1); | |
| 599 /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */ | |
| 600 | |
| 601 if (length+num.exponent>DECPMAX) { /* cannot fit */ | |
| 602 decFloatZero(result); | |
| 603 DFWORD(result, 0)=DECFLOAT_qNaN; | |
| 604 set->status|=DEC_Division_impossible; | |
| 605 return result; | |
| 606 } | |
| 607 | |
| 608 if (num.exponent>=0) { /* already an int, or need pad zeros */ | |
| 609 for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0; | |
| 610 num.lsd+=num.exponent; | |
| 611 } | |
| 612 else { /* too long: round or truncate needed */ | |
| 613 Int drop=-num.exponent; | |
| 614 if (!(op&REMNEAR)) { /* simple truncate */ | |
| 615 num.lsd-=drop; | |
| 616 if (num.lsd<num.msd) { /* truncated all */ | |
| 617 num.lsd=num.msd; /* make 0 */ | |
| 618 *num.lsd=0; /* .. [sign still relevant] */ | |
| 619 } | |
| 620 } | |
| 621 else { /* round to nearest even [sigh] */ | |
| 622 /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */ | |
| 623 /* (this is a special case of Quantize -- q.v. for commentary) */ | |
| 624 uByte *roundat; /* -> re-round digit */ | |
| 625 uByte reround; /* reround value */ | |
| 626 *(num.msd-1)=0; /* in case of left carry, or make 0 */ | |
| 627 if (drop<length) roundat=num.lsd-drop+1; | |
| 628 else if (drop==length) roundat=num.msd; | |
| 629 else roundat=num.msd-1; /* [-> 0] */ | |
| 630 reround=*roundat; | |
| 631 for (ub=roundat+1; ub<=num.lsd; ub++) { | |
| 632 if (*ub!=0) { | |
| 633 reround=DECSTICKYTAB[reround]; | |
| 634 break; | |
| 635 } | |
| 636 } /* check stickies */ | |
| 637 if (roundat>num.msd) num.lsd=roundat-1; | |
| 638 else { | |
| 639 num.msd--; /* use the 0 .. */ | |
| 640 num.lsd=num.msd; /* .. at the new MSD place */ | |
| 641 } | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
642 if (reround!=0) { /* discarding non-zero */ |
| 0 | 643 uInt bump=0; |
| 644 /* rounding is DEC_ROUND_HALF_EVEN always */ | |
| 645 if (reround>5) bump=1; /* >0.5 goes up */ | |
| 646 else if (reround==5) /* exactly 0.5000 .. */ | |
| 647 bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */ | |
| 648 if (bump!=0) { /* need increment */ | |
| 649 /* increment the coefficient; this might end up with 1000... */ | |
| 650 ub=num.lsd; | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
651 for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0); |
| 0 | 652 for (; *ub==9; ub--) *ub=0; /* at most 3 more */ |
| 653 *ub+=1; | |
| 654 if (ub<num.msd) num.msd--; /* carried */ | |
| 655 } /* bump needed */ | |
| 656 } /* reround!=0 */ | |
| 657 } /* remnear */ | |
| 658 } /* round or truncate needed */ | |
| 659 num.exponent=0; /* all paths */ | |
| 660 /*decShowNum(&num, "int"); */ | |
| 661 | |
| 662 if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */ | |
| 663 | |
| 664 /* Have a remainder to calculate */ | |
| 665 decFinalize("ient, &num, set); /* lay out the integer so far */ | |
| 666 DFWORD("ient, 0)^=DECFLOAT_Sign; /* negate it */ | |
| 667 sign=DFWORD(dfl, 0); /* save sign of dfl */ | |
| 668 decFloatFMA(result, "ient, dfr, dfl, set); | |
| 669 if (!DFISZERO(result)) return result; | |
| 670 /* if the result is zero the sign shall be sign of dfl */ | |
| 671 DFWORD("ient, 0)=sign; /* construct decFloat of sign */ | |
| 672 return decFloatCopySign(result, result, "ient); | |
| 673 } /* decDivide */ | |
| 674 | |
| 675 /* ------------------------------------------------------------------ */ | |
| 676 /* decFiniteMultiply -- multiply two finite decFloats */ | |
| 677 /* */ | |
| 678 /* num gets the result of multiplying dfl and dfr */ | |
| 679 /* bcdacc .. with the coefficient in this array */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
680 /* dfl is the first decFloat (lhs) */ |
| 0 | 681 /* dfr is the second decFloat (rhs) */ |
| 682 /* */ | |
| 683 /* This effects the multiplication of two decFloats, both known to be */ | |
| 684 /* finite, leaving the result in a bcdnum ready for decFinalize (for */ | |
| 685 /* use in Multiply) or in a following addition (FMA). */ | |
| 686 /* */ | |
| 687 /* bcdacc must have space for at least DECPMAX9*18+1 bytes. */ | |
| 688 /* No error is possible and no status is set. */ | |
| 689 /* ------------------------------------------------------------------ */ | |
| 690 /* This routine has two separate implementations of the core */ | |
| 691 /* multiplication; both using base-billion. One uses only 32-bit */ | |
| 692 /* variables (Ints and uInts) or smaller; the other uses uLongs (for */ | |
| 693 /* multiplication and addition only). Both implementations cover */ | |
| 694 /* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
695 /* comparisons. In any one compilation only one implementation for */ |
| 0 | 696 /* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */ |
| 697 /* version is forced. */ | |
| 698 /* */ | |
| 699 /* Historical note: an earlier version of this code also supported the */ | |
| 700 /* 256-bit format and has been preserved. That is somewhat trickier */ | |
| 701 /* during lazy carry splitting because the initial quotient estimate */ | |
| 702 /* (est) can exceed 32 bits. */ | |
| 703 | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
704 #define MULTBASE ((uInt)BILLION) /* the base used for multiply */ |
| 0 | 705 #define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */ |
| 706 #define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */ | |
| 707 #define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */ | |
| 708 | |
| 709 /* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */ | |
| 710 #if DECEMAXD>9 | |
| 711 #error Exponent may overflow when doubled for Multiply | |
| 712 #endif | |
| 713 #if MULACCLEN!=(MULACCLEN/4)*4 | |
| 714 /* This assumption is used below only for initialization */ | |
| 715 #error MULACCLEN is not a multiple of 4 | |
| 716 #endif | |
| 717 | |
| 718 static void decFiniteMultiply(bcdnum *num, uByte *bcdacc, | |
| 719 const decFloat *dfl, const decFloat *dfr) { | |
| 720 uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */ | |
| 721 uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */ | |
| 722 uInt *ui, *uj; /* work */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
723 uByte *ub; /* .. */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
724 uInt uiwork; /* for macros */ |
| 0 | 725 |
| 726 #if DECUSE64 | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
727 uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
728 uLong *pl; /* work -> lazy accumulator */ |
| 0 | 729 uInt acc[MULACCLEN]; /* coefficent in base-billion .. */ |
| 730 #else | |
| 731 uInt acc[MULACCLEN*2]; /* accumulator in base-billion .. */ | |
| 732 #endif | |
| 733 uInt *pa; /* work -> accumulator */ | |
| 734 /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */ | |
| 735 | |
| 736 /* Calculate sign and exponent */ | |
| 737 num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign; | |
| 738 num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */ | |
| 739 | |
| 740 /* Extract the coefficients and prepare the accumulator */ | |
| 741 /* the coefficients of the operands are decoded into base-billion */ | |
| 742 /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */ | |
| 743 /* appropriate size. */ | |
| 744 GETCOEFFBILL(dfl, bufl); | |
| 745 GETCOEFFBILL(dfr, bufr); | |
| 746 #if DECTRACE && 0 | |
| 747 printf("CoeffbL:"); | |
| 748 for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui); | |
| 749 printf("\n"); | |
| 750 printf("CoeffbR:"); | |
| 751 for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj); | |
| 752 printf("\n"); | |
| 753 #endif | |
| 754 | |
| 755 /* start the 64-bit/32-bit differing paths... */ | |
| 756 #if DECUSE64 | |
| 757 | |
| 758 /* zero the accumulator */ | |
| 759 #if MULACCLEN==4 | |
| 760 accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0; | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
761 #else /* use a loop */ |
| 0 | 762 /* MULACCLEN is a multiple of four, asserted above */ |
| 763 for (pl=accl; pl<accl+MULACCLEN; pl+=4) { | |
| 764 *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */ | |
| 765 } /* pl */ | |
| 766 #endif | |
| 767 | |
| 768 /* Effect the multiplication */ | |
| 769 /* The multiplcation proceeds using MFC's lazy-carry resolution */ | |
| 770 /* algorithm from decNumber. First, the multiplication is */ | |
| 771 /* effected, allowing accumulation of the partial products (which */ | |
| 772 /* are in base-billion at each column position) into 64 bits */ | |
| 773 /* without resolving back to base=billion after each addition. */ | |
| 774 /* These 64-bit numbers (which may contain up to 19 decimal digits) */ | |
| 775 /* are then split using the Clark & Cowlishaw algorithm (see below). */ | |
| 776 /* [Testing for 0 in the inner loop is not really a 'win'] */ | |
| 777 for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */ | |
| 778 if (*ui==0) continue; /* product cannot affect result */ | |
| 779 pl=accl+(ui-bufr); /* where to add the lhs */ | |
| 780 for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */ | |
| 781 /* if (*uj==0) continue; // product cannot affect result */ | |
| 782 *pl+=((uLong)*ui)*(*uj); | |
| 783 } /* uj */ | |
| 784 } /* ui */ | |
| 785 | |
| 786 /* The 64-bit carries must now be resolved; this means that a */ | |
| 787 /* quotient/remainder has to be calculated for base-billion (1E+9). */ | |
| 788 /* For this, Clark & Cowlishaw's quotient estimation approach (also */ | |
| 789 /* used in decNumber) is needed, because 64-bit divide is generally */ | |
| 790 /* extremely slow on 32-bit machines, and may be slower than this */ | |
| 791 /* approach even on 64-bit machines. This algorithm splits X */ | |
| 792 /* using: */ | |
| 793 /* */ | |
| 794 /* magic=2**(A+B)/1E+9; // 'magic number' */ | |
| 795 /* hop=X/2**A; // high order part of X (by shift) */ | |
| 796 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */ | |
| 797 /* */ | |
| 798 /* A and B are quite constrained; hop and magic must fit in 32 bits, */ | |
| 799 /* and 2**(A+B) must be as large as possible (which is 2**61 if */ | |
| 800 /* magic is to fit). Further, maxX increases with the length of */ | |
| 801 /* the operands (and hence the number of partial products */ | |
| 802 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */ | |
| 803 /* */ | |
| 804 /* It can be shown that when OPLEN is 2 then the maximum error in */ | |
| 805 /* the estimated quotient is <1, but for larger maximum x the */ | |
| 806 /* maximum error is above 1 so a correction that is >1 may be */ | |
| 807 /* needed. Values of A and B are chosen to satisfy the constraints */ | |
| 808 /* just mentioned while minimizing the maximum error (and hence the */ | |
| 809 /* maximum correction), as shown in the following table: */ | |
| 810 /* */ | |
| 811 /* Type OPLEN A B maxX maxError maxCorrection */ | |
| 812 /* --------------------------------------------------------- */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
813 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
814 /* QUAD 4 30 31 <4*10**18 1.17 2 */ |
| 0 | 815 /* */ |
| 816 /* In the OPLEN==2 case there is most choice, but the value for B */ | |
| 817 /* of 32 has a big advantage as then the calculation of the */ | |
| 818 /* estimate requires no shifting; the compiler can extract the high */ | |
| 819 /* word directly after multiplying magic*hop. */ | |
| 820 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */ | |
| 821 #if DOUBLE | |
| 822 #define MULSHIFTA 29 | |
| 823 #define MULSHIFTB 32 | |
| 824 #elif QUAD | |
| 825 #define MULSHIFTA 30 | |
| 826 #define MULSHIFTB 31 | |
| 827 #else | |
| 828 #error Unexpected type | |
| 829 #endif | |
| 830 | |
| 831 #if DECTRACE | |
| 832 printf("MulAccl:"); | |
| 833 for (pl=accl+MULACCLEN-1; pl>=accl; pl--) | |
| 834 printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff)); | |
| 835 printf("\n"); | |
| 836 #endif | |
| 837 | |
| 838 for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */ | |
| 839 uInt lo, hop; /* work */ | |
| 840 uInt est; /* cannot exceed 4E+9 */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
841 if (*pl>=MULTBASE) { |
| 0 | 842 /* *pl holds a binary number which needs to be split */ |
| 843 hop=(uInt)(*pl>>MULSHIFTA); | |
| 844 est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB); | |
| 845 /* the estimate is now in est; now calculate hi:lo-est*10**9; */ | |
| 846 /* happily the top word of the result is irrelevant because it */ | |
| 847 /* will always be zero so this needs only one multiplication */ | |
| 848 lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */ | |
| 849 /* If QUAD, the correction here could be +2 */ | |
| 850 if (lo>=MULTBASE) { | |
| 851 lo-=MULTBASE; /* correct by +1 */ | |
| 852 est++; | |
| 853 #if QUAD | |
| 854 /* may need to correct by +2 */ | |
| 855 if (lo>=MULTBASE) { | |
| 856 lo-=MULTBASE; | |
| 857 est++; | |
| 858 } | |
| 859 #endif | |
| 860 } | |
| 861 /* finally place lo as the new coefficient 'digit' and add est to */ | |
| 862 /* the next place up [this is safe because this path is never */ | |
| 863 /* taken on the final iteration as *pl will fit] */ | |
| 864 *pa=lo; | |
| 865 *(pl+1)+=est; | |
| 866 } /* *pl needed split */ | |
| 867 else { /* *pl<MULTBASE */ | |
| 868 *pa=(uInt)*pl; /* just copy across */ | |
| 869 } | |
| 870 } /* pl loop */ | |
| 871 | |
| 872 #else /* 32-bit */ | |
| 873 for (pa=acc;; pa+=4) { /* zero the accumulator */ | |
| 874 *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */ | |
| 875 if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */ | |
| 876 } /* pa */ | |
| 877 | |
| 878 /* Effect the multiplication */ | |
| 879 /* uLongs are not available (and in particular, there is no uLong */ | |
| 880 /* divide) but it is still possible to use MFC's lazy-carry */ | |
| 881 /* resolution algorithm from decNumber. First, the multiplication */ | |
| 882 /* is effected, allowing accumulation of the partial products */ | |
| 883 /* (which are in base-billion at each column position) into 64 bits */ | |
| 884 /* [with the high-order 32 bits in each position being held at */ | |
| 885 /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */ | |
| 886 /* These 64-bit numbers (which may contain up to 19 decimal digits) */ | |
| 887 /* are then split using the Clark & Cowlishaw algorithm (see */ | |
| 888 /* below). */ | |
| 889 for (ui=bufr;; ui++) { /* over each item in rhs */ | |
| 890 uInt hi, lo; /* words of exact multiply result */ | |
| 891 pa=acc+(ui-bufr); /* where to add the lhs */ | |
| 892 for (uj=bufl;; uj++, pa++) { /* over each item in lhs */ | |
| 893 LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */ | |
| 894 lo=(*ui)*(*uj); /* .. */ | |
| 895 *pa+=lo; /* accumulate low bits and .. */ | |
| 896 *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */ | |
| 897 if (uj==bufl+MULOPLEN-1) break; | |
| 898 } | |
| 899 if (ui==bufr+MULOPLEN-1) break; | |
| 900 } | |
| 901 | |
| 902 /* The 64-bit carries must now be resolved; this means that a */ | |
| 903 /* quotient/remainder has to be calculated for base-billion (1E+9). */ | |
| 904 /* For this, Clark & Cowlishaw's quotient estimation approach (also */ | |
| 905 /* used in decNumber) is needed, because 64-bit divide is generally */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
906 /* extremely slow on 32-bit machines. This algorithm splits X */ |
| 0 | 907 /* using: */ |
| 908 /* */ | |
| 909 /* magic=2**(A+B)/1E+9; // 'magic number' */ | |
| 910 /* hop=X/2**A; // high order part of X (by shift) */ | |
| 911 /* est=magic*hop/2**B // quotient estimate (may be low by 1) */ | |
| 912 /* */ | |
| 913 /* A and B are quite constrained; hop and magic must fit in 32 bits, */ | |
| 914 /* and 2**(A+B) must be as large as possible (which is 2**61 if */ | |
| 915 /* magic is to fit). Further, maxX increases with the length of */ | |
| 916 /* the operands (and hence the number of partial products */ | |
| 917 /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */ | |
| 918 /* */ | |
| 919 /* It can be shown that when OPLEN is 2 then the maximum error in */ | |
| 920 /* the estimated quotient is <1, but for larger maximum x the */ | |
| 921 /* maximum error is above 1 so a correction that is >1 may be */ | |
| 922 /* needed. Values of A and B are chosen to satisfy the constraints */ | |
| 923 /* just mentioned while minimizing the maximum error (and hence the */ | |
| 924 /* maximum correction), as shown in the following table: */ | |
| 925 /* */ | |
| 926 /* Type OPLEN A B maxX maxError maxCorrection */ | |
| 927 /* --------------------------------------------------------- */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
928 /* DOUBLE 2 29 32 <2*10**18 0.63 1 */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
929 /* QUAD 4 30 31 <4*10**18 1.17 2 */ |
| 0 | 930 /* */ |
| 931 /* In the OPLEN==2 case there is most choice, but the value for B */ | |
| 932 /* of 32 has a big advantage as then the calculation of the */ | |
| 933 /* estimate requires no shifting; the high word is simply */ | |
| 934 /* calculated from multiplying magic*hop. */ | |
| 935 #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */ | |
| 936 #if DOUBLE | |
| 937 #define MULSHIFTA 29 | |
| 938 #define MULSHIFTB 32 | |
| 939 #elif QUAD | |
| 940 #define MULSHIFTA 30 | |
| 941 #define MULSHIFTB 31 | |
| 942 #else | |
| 943 #error Unexpected type | |
| 944 #endif | |
| 945 | |
| 946 #if DECTRACE | |
| 947 printf("MulHiLo:"); | |
| 948 for (pa=acc+MULACCLEN-1; pa>=acc; pa--) | |
| 949 printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa); | |
| 950 printf("\n"); | |
| 951 #endif | |
| 952 | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
953 for (pa=acc;; pa++) { /* each low uInt */ |
| 0 | 954 uInt hi, lo; /* words of exact multiply result */ |
| 955 uInt hop, estlo; /* work */ | |
| 956 #if QUAD | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
957 uInt esthi; /* .. */ |
| 0 | 958 #endif |
| 959 | |
| 960 lo=*pa; | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
961 hi=*(pa+MULACCLEN); /* top 32 bits */ |
| 0 | 962 /* hi and lo now hold a binary number which needs to be split */ |
| 963 | |
| 964 #if DOUBLE | |
| 965 hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */ | |
| 966 LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */ | |
| 967 /* [MULSHIFTB is 32, so estlo can be used directly] */ | |
| 968 /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */ | |
| 969 /* happily the top word of the result is irrelevant because it */ | |
| 970 /* will always be zero so this needs only one multiplication */ | |
| 971 lo-=(estlo*MULTBASE); | |
| 972 /* esthi=0; // high word is ignored below */ | |
| 973 /* the correction here will be at most +1; do it */ | |
| 974 if (lo>=MULTBASE) { | |
| 975 lo-=MULTBASE; | |
| 976 estlo++; | |
| 977 } | |
| 978 #elif QUAD | |
| 979 hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */ | |
| 980 LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */ | |
| 981 estlo=hop*MULMAGIC; /* .. so low word needed */ | |
| 982 estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */ | |
| 983 /* esthi=0; // high word is ignored below */ | |
| 984 lo-=(estlo*MULTBASE); /* as above */ | |
| 985 /* the correction here could be +1 or +2 */ | |
| 986 if (lo>=MULTBASE) { | |
| 987 lo-=MULTBASE; | |
| 988 estlo++; | |
| 989 } | |
| 990 if (lo>=MULTBASE) { | |
| 991 lo-=MULTBASE; | |
| 992 estlo++; | |
| 993 } | |
| 994 #else | |
| 995 #error Unexpected type | |
| 996 #endif | |
| 997 | |
| 998 /* finally place lo as the new accumulator digit and add est to */ | |
| 999 /* the next place up; this latter add could cause a carry of 1 */ | |
| 1000 /* to the high word of the next place */ | |
| 1001 *pa=lo; | |
| 1002 *(pa+1)+=estlo; | |
| 1003 /* esthi is always 0 for DOUBLE and QUAD so this is skipped */ | |
| 1004 /* *(pa+1+MULACCLEN)+=esthi; */ | |
| 1005 if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */ | |
| 1006 if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */ | |
| 1007 } /* pa loop */ | |
| 1008 #endif | |
| 1009 | |
| 1010 /* At this point, whether using the 64-bit or the 32-bit paths, the */ | |
| 1011 /* accumulator now holds the (unrounded) result in base-billion; */ | |
| 1012 /* one base-billion 'digit' per uInt. */ | |
| 1013 #if DECTRACE | |
| 1014 printf("MultAcc:"); | |
| 1015 for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa); | |
| 1016 printf("\n"); | |
| 1017 #endif | |
| 1018 | |
| 1019 /* Now convert to BCD for rounding and cleanup, starting from the */ | |
| 1020 /* most significant end */ | |
| 1021 pa=acc+MULACCLEN-1; | |
| 1022 if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */ | |
| 1023 else { /* >=1 word of leading zeros */ | |
| 1024 num->msd=bcdacc; /* known leading zeros are gone */ | |
| 1025 pa--; /* skip first word .. */ | |
| 1026 for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */ | |
| 1027 } | |
| 1028 for (ub=bcdacc;; pa--, ub+=9) { | |
| 1029 if (*pa!=0) { /* split(s) needed */ | |
| 1030 uInt top, mid, rem; /* work */ | |
| 1031 /* *pa is non-zero -- split the base-billion acc digit into */ | |
| 1032 /* hi, mid, and low three-digits */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1033 #define mulsplit9 1000000 /* divisor */ |
| 0 | 1034 #define mulsplit6 1000 /* divisor */ |
| 1035 /* The splitting is done by simple divides and remainders, */ | |
| 1036 /* assuming the compiler will optimize these where useful */ | |
| 1037 /* [GCC does] */ | |
| 1038 top=*pa/mulsplit9; | |
| 1039 rem=*pa%mulsplit9; | |
| 1040 mid=rem/mulsplit6; | |
| 1041 rem=rem%mulsplit6; | |
| 1042 /* lay out the nine BCD digits (plus one unwanted byte) */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1043 UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4])); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1044 UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4])); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1045 UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4])); |
| 0 | 1046 } |
| 1047 else { /* *pa==0 */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1048 UBFROMUI(ub, 0); /* clear 9 BCD8s */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1049 UBFROMUI(ub+4, 0); /* .. */ |
| 0 | 1050 *(ub+8)=0; /* .. */ |
| 1051 } | |
| 1052 if (pa==acc) break; | |
| 1053 } /* BCD conversion loop */ | |
| 1054 | |
| 1055 num->lsd=ub+8; /* complete the bcdnum .. */ | |
| 1056 | |
| 1057 #if DECTRACE | |
| 1058 decShowNum(num, "postmult"); | |
| 1059 decFloatShow(dfl, "dfl"); | |
| 1060 decFloatShow(dfr, "dfr"); | |
| 1061 #endif | |
| 1062 return; | |
| 1063 } /* decFiniteMultiply */ | |
| 1064 | |
| 1065 /* ------------------------------------------------------------------ */ | |
| 1066 /* decFloatAbs -- absolute value, heeding NaNs, etc. */ | |
| 1067 /* */ | |
| 1068 /* result gets the canonicalized df with sign 0 */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1069 /* df is the decFloat to abs */ |
| 0 | 1070 /* set is the context */ |
| 1071 /* returns result */ | |
| 1072 /* */ | |
| 1073 /* This has the same effect as decFloatPlus unless df is negative, */ | |
| 1074 /* in which case it has the same effect as decFloatMinus. The */ | |
| 1075 /* effect is also the same as decFloatCopyAbs except that NaNs are */ | |
| 1076 /* handled normally (the sign of a NaN is not affected, and an sNaN */ | |
| 1077 /* will signal) and the result will be canonical. */ | |
| 1078 /* ------------------------------------------------------------------ */ | |
| 1079 decFloat * decFloatAbs(decFloat *result, const decFloat *df, | |
| 1080 decContext *set) { | |
| 1081 if (DFISNAN(df)) return decNaNs(result, df, NULL, set); | |
| 1082 decCanonical(result, df); /* copy and check */ | |
| 1083 DFBYTE(result, 0)&=~0x80; /* zero sign bit */ | |
| 1084 return result; | |
| 1085 } /* decFloatAbs */ | |
| 1086 | |
| 1087 /* ------------------------------------------------------------------ */ | |
| 1088 /* decFloatAdd -- add two decFloats */ | |
| 1089 /* */ | |
| 1090 /* result gets the result of adding dfl and dfr: */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1091 /* dfl is the first decFloat (lhs) */ |
| 0 | 1092 /* dfr is the second decFloat (rhs) */ |
| 1093 /* set is the context */ | |
| 1094 /* returns result */ | |
| 1095 /* */ | |
| 1096 /* ------------------------------------------------------------------ */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1097 #if QUAD |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1098 /* Table for testing MSDs for fastpath elimination; returns the MSD of */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1099 /* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1100 /* Infinities return -32 and NaNs return -128 so that summing the two */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1101 /* MSDs also allows rapid tests for the Specials (see code below). */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1102 const Int DECTESTMSD[64]={ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1103 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1104 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128, |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1105 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1106 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128}; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1107 #else |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1108 /* The table for testing MSDs is shared between the modules */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1109 extern const Int DECTESTMSD[64]; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1110 #endif |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1111 |
| 0 | 1112 decFloat * decFloatAdd(decFloat *result, |
| 1113 const decFloat *dfl, const decFloat *dfr, | |
| 1114 decContext *set) { | |
| 1115 bcdnum num; /* for final conversion */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1116 Int bexpl, bexpr; /* left and right biased exponents */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1117 uByte *ub, *us, *ut; /* work */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1118 uInt uiwork; /* for macros */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1119 #if QUAD |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1120 uShort uswork; /* .. */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1121 #endif |
| 0 | 1122 |
| 1123 uInt sourhil, sourhir; /* top words from source decFloats */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1124 /* [valid only through end of */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1125 /* fastpath code -- before swap] */ |
| 0 | 1126 uInt diffsign; /* non-zero if signs differ */ |
| 1127 uInt carry; /* carry: 0 or 1 before add loop */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1128 Int overlap; /* coefficient overlap (if full) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1129 Int summ; /* sum of the MSDs */ |
| 0 | 1130 /* the following buffers hold coefficients with various alignments */ |
| 1131 /* (see commentary and diagrams below) */ | |
| 1132 uByte acc[4+2+DECPMAX*3+8]; | |
| 1133 uByte buf[4+2+DECPMAX*2]; | |
| 1134 uByte *umsd, *ulsd; /* local MSD and LSD pointers */ | |
| 1135 | |
| 1136 #if DECLITEND | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1137 #define CARRYPAT 0x01000000 /* carry=1 pattern */ |
| 0 | 1138 #else |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1139 #define CARRYPAT 0x00000001 /* carry=1 pattern */ |
| 0 | 1140 #endif |
| 1141 | |
| 1142 /* Start decoding the arguments */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1143 /* The initial exponents are placed into the opposite Ints to */ |
| 0 | 1144 /* that which might be expected; there are two sets of data to */ |
| 1145 /* keep track of (each decFloat and the corresponding exponent), */ | |
| 1146 /* and this scheme means that at the swap point (after comparing */ | |
| 1147 /* exponents) only one pair of words needs to be swapped */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1148 /* whichever path is taken (thereby minimising worst-case path). */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1149 /* The calculated exponents will be nonsense when the arguments are */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1150 /* Special, but are not used in that path */ |
| 0 | 1151 sourhil=DFWORD(dfl, 0); /* LHS top word */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1152 summ=DECTESTMSD[sourhil>>26]; /* get first MSD for testing */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1153 bexpr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1154 bexpr+=GETECON(dfl); /* .. + continuation */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1155 |
| 0 | 1156 sourhir=DFWORD(dfr, 0); /* RHS top word */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1157 summ+=DECTESTMSD[sourhir>>26]; /* sum MSDs for testing */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1158 bexpl=DECCOMBEXP[sourhir>>26]; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1159 bexpl+=GETECON(dfr); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1160 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1161 /* here bexpr has biased exponent from lhs, and vice versa */ |
| 0 | 1162 |
| 1163 diffsign=(sourhil^sourhir)&DECFLOAT_Sign; | |
| 1164 | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1165 /* now determine whether to take a fast path or the full-function */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1166 /* slow path. The slow path must be taken when: */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1167 /* -- both numbers are finite, and: */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1168 /* the exponents are different, or */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1169 /* the signs are different, or */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1170 /* the sum of the MSDs is >8 (hence might overflow) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1171 /* specialness and the sum of the MSDs can be tested at once using */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1172 /* the summ value just calculated, so the test for specials is no */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1173 /* longer on the worst-case path (as of 3.60) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1174 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1175 if (summ<=8) { /* MSD+MSD is good, or there is a special */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1176 if (summ<0) { /* there is a special */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1177 /* Inf+Inf would give -64; Inf+finite is -32 or higher */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1178 if (summ<-64) return decNaNs(result, dfl, dfr, set); /* one or two NaNs */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1179 /* two infinities with different signs is invalid */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1180 if (summ==-64 && diffsign) return decInvalid(result, set); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1181 if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1182 return decInfinity(result, dfr); /* RHS must be Inf */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1183 } |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1184 /* Here when both arguments are finite; fast path is possible */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1185 /* (currently only for aligned and same-sign) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1186 if (bexpr==bexpl && !diffsign) { |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1187 uInt tac[DECLETS+1]; /* base-1000 coefficient */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1188 uInt encode; /* work */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1189 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1190 /* Get one coefficient as base-1000 and add the other */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1191 GETCOEFFTHOU(dfl, tac); /* least-significant goes to [0] */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1192 ADDCOEFFTHOU(dfr, tac); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1193 /* here the sum of the MSDs (plus any carry) will be <10 due to */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1194 /* the fastpath test earlier */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1195 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1196 /* construct the result; low word is the same for both formats */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1197 encode =BIN2DPD[tac[0]]; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1198 encode|=BIN2DPD[tac[1]]<<10; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1199 encode|=BIN2DPD[tac[2]]<<20; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1200 encode|=BIN2DPD[tac[3]]<<30; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1201 DFWORD(result, (DECBYTES/4)-1)=encode; |
| 0 | 1202 |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1203 /* collect next two declets (all that remains, for Double) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1204 encode =BIN2DPD[tac[3]]>>2; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1205 encode|=BIN2DPD[tac[4]]<<8; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1206 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1207 #if QUAD |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1208 /* complete and lay out middling words */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1209 encode|=BIN2DPD[tac[5]]<<18; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1210 encode|=BIN2DPD[tac[6]]<<28; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1211 DFWORD(result, 2)=encode; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1212 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1213 encode =BIN2DPD[tac[6]]>>4; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1214 encode|=BIN2DPD[tac[7]]<<6; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1215 encode|=BIN2DPD[tac[8]]<<16; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1216 encode|=BIN2DPD[tac[9]]<<26; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1217 DFWORD(result, 1)=encode; |
| 0 | 1218 |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1219 /* and final two declets */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1220 encode =BIN2DPD[tac[9]]>>6; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1221 encode|=BIN2DPD[tac[10]]<<4; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1222 #endif |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1223 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1224 /* add exponent continuation and sign (from either argument) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1225 encode|=sourhil & (ECONMASK | DECFLOAT_Sign); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1226 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1227 /* create lookup index = MSD + top two bits of biased exponent <<4 */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1228 tac[DECLETS]|=(bexpl>>DECECONL)<<4; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1229 encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1230 DFWORD(result, 0)=encode; /* complete */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1231 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1232 /* decFloatShow(result, ">"); */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1233 return result; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1234 } /* fast path OK */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1235 /* drop through to slow path */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1236 } /* low sum or Special(s) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1237 |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1238 /* Slow path required -- arguments are finite and might overflow, */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1239 /* or require alignment, or might have different signs */ |
| 0 | 1240 |
| 1241 /* now swap either exponents or argument pointers */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1242 if (bexpl<=bexpr) { |
| 0 | 1243 /* original left is bigger */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1244 Int bexpswap=bexpl; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1245 bexpl=bexpr; |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1246 bexpr=bexpswap; |
| 0 | 1247 /* printf("left bigger\n"); */ |
| 1248 } | |
| 1249 else { | |
| 1250 const decFloat *dfswap=dfl; | |
| 1251 dfl=dfr; | |
| 1252 dfr=dfswap; | |
| 1253 /* printf("right bigger\n"); */ | |
| 1254 } | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1255 /* [here dfl and bexpl refer to the datum with the larger exponent, */ |
| 0 | 1256 /* of if the exponents are equal then the original LHS argument] */ |
| 1257 | |
| 1258 /* if lhs is zero then result will be the rhs (now known to have */ | |
| 1259 /* the smaller exponent), which also may need to be tested for zero */ | |
| 1260 /* for the weird IEEE 754 sign rules */ | |
| 1261 if (DFISZERO(dfl)) { | |
| 1262 decCanonical(result, dfr); /* clean copy */ | |
| 1263 /* "When the sum of two operands with opposite signs is */ | |
| 1264 /* exactly zero, the sign of that sum shall be '+' in all */ | |
| 1265 /* rounding modes except round toward -Infinity, in which */ | |
| 1266 /* mode that sign shall be '-'." */ | |
| 1267 if (diffsign && DFISZERO(result)) { | |
| 1268 DFWORD(result, 0)&=~DECFLOAT_Sign; /* assume sign 0 */ | |
| 1269 if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign; | |
| 1270 } | |
| 1271 return result; | |
| 1272 } /* numfl is zero */ | |
| 1273 /* [here, LHS is non-zero; code below assumes that] */ | |
| 1274 | |
| 1275 /* Coefficients layout during the calculations to follow: */ | |
| 1276 /* */ | |
| 1277 /* Overlap case: */ | |
| 1278 /* +------------------------------------------------+ */ | |
| 1279 /* acc: |0000| coeffa | tail B | | */ | |
| 1280 /* +------------------------------------------------+ */ | |
| 1281 /* buf: |0000| pad0s | coeffb | | */ | |
| 1282 /* +------------------------------------------------+ */ | |
| 1283 /* */ | |
| 1284 /* Touching coefficients or gap: */ | |
| 1285 /* +------------------------------------------------+ */ | |
| 1286 /* acc: |0000| coeffa | gap | coeffb | */ | |
| 1287 /* +------------------------------------------------+ */ | |
| 1288 /* [buf not used or needed; gap clamped to Pmax] */ | |
| 1289 | |
| 1290 /* lay out lhs coefficient into accumulator; this starts at acc+4 */ | |
| 1291 /* for decDouble or acc+6 for decQuad so the LSD is word- */ | |
| 1292 /* aligned; the top word gap is there only in case a carry digit */ | |
| 1293 /* is prefixed after the add -- it does not need to be zeroed */ | |
| 1294 #if DOUBLE | |
| 1295 #define COFF 4 /* offset into acc */ | |
| 1296 #elif QUAD | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1297 UBFROMUS(acc+4, 0); /* prefix 00 */ |
| 0 | 1298 #define COFF 6 /* offset into acc */ |
| 1299 #endif | |
| 1300 | |
| 1301 GETCOEFF(dfl, acc+COFF); /* decode from decFloat */ | |
| 1302 ulsd=acc+COFF+DECPMAX-1; | |
| 1303 umsd=acc+4; /* [having this here avoids */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1304 |
| 0 | 1305 #if DECTRACE |
| 1306 {bcdnum tum; | |
| 1307 tum.msd=umsd; | |
| 1308 tum.lsd=ulsd; | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1309 tum.exponent=bexpl-DECBIAS; |
| 0 | 1310 tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign; |
| 1311 decShowNum(&tum, "dflx");} | |
| 1312 #endif | |
| 1313 | |
| 1314 /* if signs differ, take ten's complement of lhs (here the */ | |
| 1315 /* coefficient is subtracted from all-nines; the 1 is added during */ | |
| 1316 /* the later add cycle -- zeros to the right do not matter because */ | |
| 1317 /* the complement of zero is zero); these are fixed-length inverts */ | |
| 1318 /* where the lsd is known to be at a 4-byte boundary (so no borrow */ | |
| 1319 /* possible) */ | |
| 1320 carry=0; /* assume no carry */ | |
| 1321 if (diffsign) { | |
| 1322 carry=CARRYPAT; /* for +1 during add */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1323 UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1324 UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1325 UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1326 UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16)); |
| 0 | 1327 #if QUAD |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1328 UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1329 UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1330 UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1331 UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1332 UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36)); |
| 0 | 1333 #endif |
| 1334 } /* diffsign */ | |
| 1335 | |
| 1336 /* now process the rhs coefficient; if it cannot overlap lhs then */ | |
| 1337 /* it can be put straight into acc (with an appropriate gap, if */ | |
| 1338 /* needed) because no actual addition will be needed (except */ | |
| 1339 /* possibly to complete ten's complement) */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1340 overlap=DECPMAX-(bexpl-bexpr); |
| 0 | 1341 #if DECTRACE |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1342 printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS)); |
| 0 | 1343 printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry); |
| 1344 #endif | |
| 1345 | |
| 1346 if (overlap<=0) { /* no overlap possible */ | |
| 1347 uInt gap; /* local work */ | |
| 1348 /* since a full addition is not needed, a ten's complement */ | |
| 1349 /* calculation started above may need to be completed */ | |
| 1350 if (carry) { | |
| 1351 for (ub=ulsd; *ub==9; ub--) *ub=0; | |
| 1352 *ub+=1; | |
| 1353 carry=0; /* taken care of */ | |
| 1354 } | |
| 1355 /* up to DECPMAX-1 digits of the final result can extend down */ | |
| 1356 /* below the LSD of the lhs, so if the gap is >DECPMAX then the */ | |
| 1357 /* rhs will be simply sticky bits. In this case the gap is */ | |
| 1358 /* clamped to DECPMAX and the exponent adjusted to suit [this is */ | |
| 1359 /* safe because the lhs is non-zero]. */ | |
| 1360 gap=-overlap; | |
| 1361 if (gap>DECPMAX) { | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1362 bexpr+=gap-1; |
| 0 | 1363 gap=DECPMAX; |
| 1364 } | |
| 1365 ub=ulsd+gap+1; /* where MSD will go */ | |
| 1366 /* Fill the gap with 0s; note that there is no addition to do */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1367 ut=acc+COFF+DECPMAX; /* start of gap */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1368 for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */ |
| 0 | 1369 if (overlap<-DECPMAX) { /* gap was > DECPMAX */ |
| 1370 *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */ | |
| 1371 } | |
| 1372 else { /* need full coefficient */ | |
| 1373 GETCOEFF(dfr, ub); /* decode from decFloat */ | |
| 1374 ub+=DECPMAX-1; /* new LSD... */ | |
| 1375 } | |
| 1376 ulsd=ub; /* save new LSD */ | |
| 1377 } /* no overlap possible */ | |
| 1378 | |
| 1379 else { /* overlap>0 */ | |
| 1380 /* coefficients overlap (perhaps completely, although also */ | |
| 1381 /* perhaps only where zeros) */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1382 if (overlap==DECPMAX) { /* aligned */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1383 ub=buf+COFF; /* where msd will go */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1384 #if QUAD |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1385 UBFROMUS(buf+4, 0); /* clear quad's 00 */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1386 #endif |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1387 GETCOEFF(dfr, ub); /* decode from decFloat */ |
| 0 | 1388 } |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1389 else { /* unaligned */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1390 ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1391 /* Fill the prefix gap with 0s; 8 will cover most common */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1392 /* unalignments, so start with direct assignments (a loop is */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1393 /* then used for any remaining -- the loop (and the one in a */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1394 /* moment) is not then on the critical path because the number */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1395 /* of additions is reduced by (at least) two in this case) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1396 UBFROMUI(buf+4, 0); /* [clears decQuad 00 too] */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1397 UBFROMUI(buf+8, 0); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1398 if (ub>buf+12) { |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1399 ut=buf+12; /* start any remaining */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1400 for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1401 } |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1402 GETCOEFF(dfr, ub); /* decode from decFloat */ |
| 0 | 1403 |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1404 /* now move tail of rhs across to main acc; again use direct */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1405 /* copies for 8 digits-worth */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1406 UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1407 UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1408 if (buf+COFF+DECPMAX+8<ub+DECPMAX) { |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1409 us=buf+COFF+DECPMAX+8; /* source */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1410 ut=acc+COFF+DECPMAX+8; /* target */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1411 for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us)); |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1412 } |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1413 } /* unaligned */ |
| 0 | 1414 |
| 1415 ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */ | |
| 1416 | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1417 /* Now do the add of the non-tail; this is all nicely aligned, */ |
| 0 | 1418 /* and is over a multiple of four digits (because for Quad two */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1419 /* zero digits were added on the left); words in both acc and */ |
| 0 | 1420 /* buf (buf especially) will often be zero */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1421 /* [byte-by-byte add, here, is about 15% slower total effect than */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1422 /* the by-fours] */ |
| 0 | 1423 |
| 1424 /* Now effect the add; this is harder on a little-endian */ | |
| 1425 /* machine as the inter-digit carry cannot use the usual BCD */ | |
| 1426 /* addition trick because the bytes are loaded in the wrong order */ | |
| 1427 /* [this loop could be unrolled, but probably scarcely worth it] */ | |
| 1428 | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1429 ut=acc+COFF+DECPMAX-4; /* target LSW (acc) */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1430 us=buf+COFF+DECPMAX-4; /* source LSW (buf, to add to acc) */ |
| 0 | 1431 |
| 1432 #if !DECLITEND | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1433 for (; ut>=acc+4; ut-=4, us-=4) { /* big-endian add loop */ |
| 0 | 1434 /* bcd8 add */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1435 carry+=UBTOUI(us); /* rhs + carry */ |
| 0 | 1436 if (carry==0) continue; /* no-op */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1437 carry+=UBTOUI(ut); /* lhs */ |
| 0 | 1438 /* Big-endian BCD adjust (uses internal carry) */ |
| 1439 carry+=0x76f6f6f6; /* note top nibble not all bits */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1440 /* apply BCD adjust and save */ |
|
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1441 UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4)); |
| 0 | 1442 carry>>=31; /* true carry was at far left */ |
| 1443 } /* add loop */ | |
| 1444 #else | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1445 for (; ut>=acc+4; ut-=4, us-=4) { /* little-endian add loop */ |
| 0 | 1446 /* bcd8 add */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1447 carry+=UBTOUI(us); /* rhs + carry */ |
| 0 | 1448 if (carry==0) continue; /* no-op [common if unaligned] */ |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1449 carry+=UBTOUI(ut); /* lhs */ |
| 0 | 1450 /* Little-endian BCD adjust; inter-digit carry must be manual */ |
| 1451 /* because the lsb from the array will be in the most-significant */ | |
| 1452 /* byte of carry */ | |
| 1453 carry+=0x76767676; /* note no inter-byte carries */ | |
| 1454 carry+=(carry & 0x80000000)>>15; | |
| 1455 carry+=(carry & 0x00800000)>>15; | |
| 1456 carry+=(carry & 0x00008000)>>15; | |
| 1457 carry-=(carry & 0x60606060)>>4; /* BCD adjust back */ | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1458 UBFROMUI(ut, carry & 0x0f0f0f0f); /* clear debris and save */ |
| 0 | 1459 /* here, final carry-out bit is at 0x00000080; move it ready */ |
| 1460 /* for next word-add (i.e., to 0x01000000) */ | |
| 1461 carry=(carry & 0x00000080)<<17; | |
| 1462 } /* add loop */ | |
| 1463 #endif | |
|
55
77e2b8dfacca
update it from 4.4.3 to 4.5.0
ryoma <e075725@ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
1464 |
| 0 | 1465 #if DECTRACE |
| 1466 {bcdnum tum; | |
|
a06113de4d67
first commit
|