You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
3909 lines
176 KiB
C
3909 lines
176 KiB
C
9 years ago
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decBasic.c -- common base code for Basic decimal types */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* Copyright (c) IBM Corporation, 2000, 2010. All rights reserved. */
|
||
|
|
/* */
|
||
|
|
/* This software is made available under the terms of the */
|
||
|
|
/* ICU License -- ICU 1.8.1 and later. */
|
||
|
|
/* */
|
||
|
|
/* The description and User's Guide ("The decNumber C Library") for */
|
||
|
|
/* this software is included in the package as decNumber.pdf. This */
|
||
|
|
/* document is also available in HTML, together with specifications, */
|
||
|
|
/* testcases, and Web links, on the General Decimal Arithmetic page. */
|
||
|
|
/* */
|
||
|
|
/* Please send comments, suggestions, and corrections to the author: */
|
||
|
|
/* mfc@uk.ibm.com */
|
||
|
|
/* Mike Cowlishaw, IBM Fellow */
|
||
|
|
/* IBM UK, PO Box 31, Birmingham Road, Warwick CV34 5JL, UK */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* This module comprises code that is shared between decDouble and */
|
||
|
|
/* decQuad (but not decSingle). The main arithmetic operations are */
|
||
|
|
/* here (Add, Subtract, Multiply, FMA, and Division operators). */
|
||
|
|
/* */
|
||
|
|
/* Unlike decNumber, parameterization takes place at compile time */
|
||
|
|
/* rather than at runtime. The parameters are set in the decDouble.c */
|
||
|
|
/* (etc.) files, which then include this one to produce the compiled */
|
||
|
|
/* code. The functions here, therefore, are code shared between */
|
||
|
|
/* multiple formats. */
|
||
|
|
/* */
|
||
|
|
/* This must be included after decCommon.c. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
// Names here refer to decFloat rather than to decDouble, etc., and
|
||
|
|
// the functions are in strict alphabetical order.
|
||
|
|
|
||
|
|
// The compile-time flags SINGLE, DOUBLE, and QUAD are set up in
|
||
|
|
// decCommon.c
|
||
|
|
#if !defined(QUAD)
|
||
|
|
#error decBasic.c must be included after decCommon.c
|
||
|
|
#endif
|
||
|
|
#if SINGLE
|
||
|
|
#error Routines in decBasic.c are for decDouble and decQuad only
|
||
|
|
#endif
|
||
|
|
|
||
|
|
/* Private constants */
|
||
|
|
#define DIVIDE 0x80000000 // Divide operations [as flags]
|
||
|
|
#define REMAINDER 0x40000000 // ..
|
||
|
|
#define DIVIDEINT 0x20000000 // ..
|
||
|
|
#define REMNEAR 0x10000000 // ..
|
||
|
|
|
||
|
|
/* Private functions (local, used only by routines in this module) */
|
||
|
|
static decFloat *decDivide(decFloat *, const decFloat *,
|
||
|
|
const decFloat *, decContext *, uInt);
|
||
|
|
static decFloat *decCanonical(decFloat *, const decFloat *);
|
||
|
|
static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
|
||
|
|
const decFloat *);
|
||
|
|
static decFloat *decInfinity(decFloat *, const decFloat *);
|
||
|
|
static decFloat *decInvalid(decFloat *, decContext *);
|
||
|
|
static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
|
||
|
|
decContext *);
|
||
|
|
static Int decNumCompare(const decFloat *, const decFloat *, Flag);
|
||
|
|
static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
|
||
|
|
enum rounding, Flag);
|
||
|
|
static uInt decToInt32(const decFloat *, decContext *, enum rounding,
|
||
|
|
Flag, Flag);
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decCanonical -- copy a decFloat, making canonical */
|
||
|
|
/* */
|
||
|
|
/* result gets the canonicalized df */
|
||
|
|
/* df is the decFloat to copy and make canonical */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is exposed via decFloatCanonical for Double and Quad only. */
|
||
|
|
/* This works on specials, too; no error or exception is possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
static decFloat * decCanonical(decFloat *result, const decFloat *df) {
|
||
|
|
uInt encode, precode, dpd; // work
|
||
|
|
uInt inword, uoff, canon; // ..
|
||
|
|
Int n; // counter (down)
|
||
|
|
if (df!=result) *result=*df; // effect copy if needed
|
||
|
|
if (DFISSPECIAL(result)) {
|
||
|
|
if (DFISINF(result)) return decInfinity(result, df); // clean Infinity
|
||
|
|
// is a NaN
|
||
|
|
DFWORD(result, 0)&=~ECONNANMASK; // clear ECON except selector
|
||
|
|
if (DFISCCZERO(df)) return result; // coefficient continuation is 0
|
||
|
|
// drop through to check payload
|
||
|
|
}
|
||
|
|
// return quickly if the coefficient continuation is canonical
|
||
|
|
{ // declare block
|
||
|
|
#if DOUBLE
|
||
|
|
uInt sourhi=DFWORD(df, 0);
|
||
|
|
uInt sourlo=DFWORD(df, 1);
|
||
|
|
if (CANONDPDOFF(sourhi, 8)
|
||
|
|
&& CANONDPDTWO(sourhi, sourlo, 30)
|
||
|
|
&& CANONDPDOFF(sourlo, 20)
|
||
|
|
&& CANONDPDOFF(sourlo, 10)
|
||
|
|
&& CANONDPDOFF(sourlo, 0)) return result;
|
||
|
|
#elif QUAD
|
||
|
|
uInt sourhi=DFWORD(df, 0);
|
||
|
|
uInt sourmh=DFWORD(df, 1);
|
||
|
|
uInt sourml=DFWORD(df, 2);
|
||
|
|
uInt sourlo=DFWORD(df, 3);
|
||
|
|
if (CANONDPDOFF(sourhi, 4)
|
||
|
|
&& CANONDPDTWO(sourhi, sourmh, 26)
|
||
|
|
&& CANONDPDOFF(sourmh, 16)
|
||
|
|
&& CANONDPDOFF(sourmh, 6)
|
||
|
|
&& CANONDPDTWO(sourmh, sourml, 28)
|
||
|
|
&& CANONDPDOFF(sourml, 18)
|
||
|
|
&& CANONDPDOFF(sourml, 8)
|
||
|
|
&& CANONDPDTWO(sourml, sourlo, 30)
|
||
|
|
&& CANONDPDOFF(sourlo, 20)
|
||
|
|
&& CANONDPDOFF(sourlo, 10)
|
||
|
|
&& CANONDPDOFF(sourlo, 0)) return result;
|
||
|
|
#endif
|
||
|
|
} // block
|
||
|
|
|
||
|
|
// Loop to repair a non-canonical coefficent, as needed
|
||
|
|
inword=DECWORDS-1; // current input word
|
||
|
|
uoff=0; // bit offset of declet
|
||
|
|
encode=DFWORD(result, inword);
|
||
|
|
for (n=DECLETS-1; n>=0; n--) { // count down declets of 10 bits
|
||
|
|
dpd=encode>>uoff;
|
||
|
|
uoff+=10;
|
||
|
|
if (uoff>32) { // crossed uInt boundary
|
||
|
|
inword--;
|
||
|
|
encode=DFWORD(result, inword);
|
||
|
|
uoff-=32;
|
||
|
|
dpd|=encode<<(10-uoff); // get pending bits
|
||
|
|
}
|
||
|
|
dpd&=0x3ff; // clear uninteresting bits
|
||
|
|
if (dpd<0x16e) continue; // must be canonical
|
||
|
|
canon=BIN2DPD[DPD2BIN[dpd]]; // determine canonical declet
|
||
|
|
if (canon==dpd) continue; // have canonical declet
|
||
|
|
// need to replace declet
|
||
|
|
if (uoff>=10) { // all within current word
|
||
|
|
encode&=~(0x3ff<<(uoff-10)); // clear the 10 bits ready for replace
|
||
|
|
encode|=canon<<(uoff-10); // insert the canonical form
|
||
|
|
DFWORD(result, inword)=encode; // .. and save
|
||
|
|
continue;
|
||
|
|
}
|
||
|
|
// straddled words
|
||
|
|
precode=DFWORD(result, inword+1); // get previous
|
||
|
|
precode&=0xffffffff>>(10-uoff); // clear top bits
|
||
|
|
DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
|
||
|
|
encode&=0xffffffff<<uoff; // clear bottom bits
|
||
|
|
encode|=canon>>(10-uoff); // insert canonical
|
||
|
|
DFWORD(result, inword)=encode; // .. and save
|
||
|
|
} // n
|
||
|
|
return result;
|
||
|
|
} // decCanonical
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decDivide -- divide operations */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of dividing dfl by dfr: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* op is the operation selector */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
#define DIVCOUNT 0 // 1 to instrument subtractions counter
|
||
|
|
#define DIVBASE ((uInt)BILLION) // the base used for divide
|
||
|
|
#define DIVOPLEN DECPMAX9 // operand length ('digits' base 10**9)
|
||
|
|
#define DIVACCLEN (DIVOPLEN*3) // accumulator length (ditto)
|
||
|
|
static decFloat * decDivide(decFloat *result, const decFloat *dfl,
|
||
|
|
const decFloat *dfr, decContext *set, uInt op) {
|
||
|
|
decFloat quotient; // for remainders
|
||
|
|
bcdnum num; // for final conversion
|
||
|
|
uInt acc[DIVACCLEN]; // coefficent in base-billion ..
|
||
|
|
uInt div[DIVOPLEN]; // divisor in base-billion ..
|
||
|
|
uInt quo[DIVOPLEN+1]; // quotient in base-billion ..
|
||
|
|
uByte bcdacc[(DIVOPLEN+1)*9+2]; // for quotient in BCD, +1, +1
|
||
|
|
uInt *msua, *msud, *msuq; // -> msu of acc, div, and quo
|
||
|
|
Int divunits, accunits; // lengths
|
||
|
|
Int quodigits; // digits in quotient
|
||
|
|
uInt *lsua, *lsuq; // -> current acc and quo lsus
|
||
|
|
Int length, multiplier; // work
|
||
|
|
uInt carry, sign; // ..
|
||
|
|
uInt *ua, *ud, *uq; // ..
|
||
|
|
uByte *ub; // ..
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
uInt divtop; // top unit of div adjusted for estimating
|
||
|
|
#if DIVCOUNT
|
||
|
|
static uInt maxcount=0; // worst-seen subtractions count
|
||
|
|
uInt divcount=0; // subtractions count [this divide]
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// calculate sign
|
||
|
|
num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
|
||
|
|
|
||
|
|
if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { // either is special?
|
||
|
|
// NaNs are handled as usual
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
// one or two infinities
|
||
|
|
if (DFISINF(dfl)) {
|
||
|
|
if (DFISINF(dfr)) return decInvalid(result, set); // Two infinities bad
|
||
|
|
if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); // as is rem
|
||
|
|
// Infinity/x is infinite and quiet, even if x=0
|
||
|
|
DFWORD(result, 0)=num.sign;
|
||
|
|
return decInfinity(result, result);
|
||
|
|
}
|
||
|
|
// must be x/Infinity -- remainders are lhs
|
||
|
|
if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
|
||
|
|
// divides: return zero with correct sign and exponent depending
|
||
|
|
// on op (Etiny for divide, 0 for divideInt)
|
||
|
|
decFloatZero(result);
|
||
|
|
if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; // add sign
|
||
|
|
else DFWORD(result, 0)=num.sign; // zeros the exponent, too
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
// next, handle zero operands (x/0 and 0/x)
|
||
|
|
if (DFISZERO(dfr)) { // x/0
|
||
|
|
if (DFISZERO(dfl)) { // 0/0 is undefined
|
||
|
|
decFloatZero(result);
|
||
|
|
DFWORD(result, 0)=DECFLOAT_qNaN;
|
||
|
|
set->status|=DEC_Division_undefined;
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); // bad rem
|
||
|
|
set->status|=DEC_Division_by_zero;
|
||
|
|
DFWORD(result, 0)=num.sign;
|
||
|
|
return decInfinity(result, result); // x/0 -> signed Infinity
|
||
|
|
}
|
||
|
|
num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); // ideal exponent
|
||
|
|
if (DFISZERO(dfl)) { // 0/x (x!=0)
|
||
|
|
// if divide, result is 0 with ideal exponent; divideInt has
|
||
|
|
// exponent=0, remainders give zero with lower exponent
|
||
|
|
if (op&DIVIDEINT) {
|
||
|
|
decFloatZero(result);
|
||
|
|
DFWORD(result, 0)|=num.sign; // add sign
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
if (!(op&DIVIDE)) { // a remainder
|
||
|
|
// exponent is the minimum of the operands
|
||
|
|
num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
|
||
|
|
// if the result is zero the sign shall be sign of dfl
|
||
|
|
num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
|
||
|
|
}
|
||
|
|
bcdacc[0]=0;
|
||
|
|
num.msd=bcdacc; // -> 0
|
||
|
|
num.lsd=bcdacc; // ..
|
||
|
|
return decFinalize(result, &num, set); // [divide may clamp exponent]
|
||
|
|
} // 0/x
|
||
|
|
// [here, both operands are known to be finite and non-zero]
|
||
|
|
|
||
|
|
// extract the operand coefficents into 'units' which are
|
||
|
|
// base-billion; the lhs is high-aligned in acc and the msu of both
|
||
|
|
// acc and div is at the right-hand end of array (offset length-1);
|
||
|
|
// the quotient can need one more unit than the operands as digits
|
||
|
|
// in it are not necessarily aligned neatly; further, the quotient
|
||
|
|
// may not start accumulating until after the end of the initial
|
||
|
|
// operand in acc if that is small (e.g., 1) so the accumulator
|
||
|
|
// must have at least that number of units extra (at the ls end)
|
||
|
|
GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
|
||
|
|
GETCOEFFBILL(dfr, div);
|
||
|
|
// zero the low uInts of acc
|
||
|
|
acc[0]=0;
|
||
|
|
acc[1]=0;
|
||
|
|
acc[2]=0;
|
||
|
|
acc[3]=0;
|
||
|
|
#if DOUBLE
|
||
|
|
#if DIVOPLEN!=2
|
||
|
|
#error Unexpected Double DIVOPLEN
|
||
|
|
#endif
|
||
|
|
#elif QUAD
|
||
|
|
acc[4]=0;
|
||
|
|
acc[5]=0;
|
||
|
|
acc[6]=0;
|
||
|
|
acc[7]=0;
|
||
|
|
#if DIVOPLEN!=4
|
||
|
|
#error Unexpected Quad DIVOPLEN
|
||
|
|
#endif
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// set msu and lsu pointers
|
||
|
|
msua=acc+DIVACCLEN-1; // [leading zeros removed below]
|
||
|
|
msuq=quo+DIVOPLEN;
|
||
|
|
//[loop for div will terminate because operands are non-zero]
|
||
|
|
for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
|
||
|
|
// the initial least-significant unit of acc is set so acc appears
|
||
|
|
// to have the same length as div.
|
||
|
|
// This moves one position towards the least possible for each
|
||
|
|
// iteration
|
||
|
|
divunits=(Int)(msud-div+1); // precalculate
|
||
|
|
lsua=msua-divunits+1; // initial working lsu of acc
|
||
|
|
lsuq=msuq; // and of quo
|
||
|
|
|
||
|
|
// set up the estimator for the multiplier; this is the msu of div,
|
||
|
|
// plus two bits from the unit below (if any) rounded up by one if
|
||
|
|
// there are any non-zero bits or units below that [the extra two
|
||
|
|
// bits makes for a much better estimate when the top unit is small]
|
||
|
|
divtop=*msud<<2;
|
||
|
|
if (divunits>1) {
|
||
|
|
uInt *um=msud-1;
|
||
|
|
uInt d=*um;
|
||
|
|
if (d>=750000000) {divtop+=3; d-=750000000;}
|
||
|
|
else if (d>=500000000) {divtop+=2; d-=500000000;}
|
||
|
|
else if (d>=250000000) {divtop++; d-=250000000;}
|
||
|
|
if (d) divtop++;
|
||
|
|
else for (um--; um>=div; um--) if (*um) {
|
||
|
|
divtop++;
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
} // >1 unit
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
{Int i;
|
||
|
|
printf("----- div=");
|
||
|
|
for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
|
||
|
|
printf("\n");}
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// now collect up to DECPMAX+1 digits in the quotient (this may
|
||
|
|
// need OPLEN+1 uInts if unaligned)
|
||
|
|
quodigits=0; // no digits yet
|
||
|
|
for (;; lsua--) { // outer loop -- each input position
|
||
|
|
#if DECCHECK
|
||
|
|
if (lsua<acc) {
|
||
|
|
printf("Acc underrun...\n");
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
#if DECTRACE
|
||
|
|
printf("Outer: quodigits=%ld acc=", (LI)quodigits);
|
||
|
|
for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
|
||
|
|
printf("\n");
|
||
|
|
#endif
|
||
|
|
*lsuq=0; // default unit result is 0
|
||
|
|
for (;;) { // inner loop -- calculate quotient unit
|
||
|
|
// strip leading zero units from acc (either there initially or
|
||
|
|
// from subtraction below); this may strip all if exactly 0
|
||
|
|
for (; *msua==0 && msua>=lsua;) msua--;
|
||
|
|
accunits=(Int)(msua-lsua+1); // [maybe 0]
|
||
|
|
// subtraction is only necessary and possible if there are as
|
||
|
|
// least as many units remaining in acc for this iteration as
|
||
|
|
// there are in div
|
||
|
|
if (accunits<divunits) {
|
||
|
|
if (accunits==0) msua++; // restore
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
|
||
|
|
// If acc is longer than div then subtraction is definitely
|
||
|
|
// possible (as msu of both is non-zero), but if they are the
|
||
|
|
// same length a comparison is needed.
|
||
|
|
// If a subtraction is needed then a good estimate of the
|
||
|
|
// multiplier for the subtraction is also needed in order to
|
||
|
|
// minimise the iterations of this inner loop because the
|
||
|
|
// subtractions needed dominate division performance.
|
||
|
|
if (accunits==divunits) {
|
||
|
|
// compare the high divunits of acc and div:
|
||
|
|
// acc<div: this quotient unit is unchanged; subtraction
|
||
|
|
// will be possible on the next iteration
|
||
|
|
// acc==div: quotient gains 1, set acc=0
|
||
|
|
// acc>div: subtraction necessary at this position
|
||
|
|
for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
|
||
|
|
// [now at first mismatch or lsu]
|
||
|
|
if (*ud>*ua) break; // next time...
|
||
|
|
if (*ud==*ua) { // all compared equal
|
||
|
|
*lsuq+=1; // increment result
|
||
|
|
msua=lsua; // collapse acc units
|
||
|
|
*msua=0; // .. to a zero
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
|
||
|
|
// subtraction necessary; estimate multiplier [see above]
|
||
|
|
// if both *msud and *msua are small it is cost-effective to
|
||
|
|
// bring in part of the following units (if any) to get a
|
||
|
|
// better estimate (assume some other non-zero in div)
|
||
|
|
#define DIVLO 1000000U
|
||
|
|
#define DIVHI (DIVBASE/DIVLO)
|
||
|
|
#if DECUSE64
|
||
|
|
if (divunits>1) {
|
||
|
|
// there cannot be a *(msud-2) for DECDOUBLE so next is
|
||
|
|
// an exact calculation unless DECQUAD (which needs to
|
||
|
|
// assume bits out there if divunits>2)
|
||
|
|
uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
|
||
|
|
uLong div=(uLong)*msud * DIVBASE + *(msud-1);
|
||
|
|
#if QUAD
|
||
|
|
if (divunits>2) div++;
|
||
|
|
#endif
|
||
|
|
mul/=div;
|
||
|
|
multiplier=(Int)mul;
|
||
|
|
}
|
||
|
|
else multiplier=*msua/(*msud);
|
||
|
|
#else
|
||
|
|
if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
|
||
|
|
multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
|
||
|
|
/(*msud*DIVHI + *(msud-1)/DIVLO +1);
|
||
|
|
}
|
||
|
|
else multiplier=(*msua<<2)/divtop;
|
||
|
|
#endif
|
||
|
|
}
|
||
|
|
else { // accunits>divunits
|
||
|
|
// msud is one unit 'lower' than msua, so estimate differently
|
||
|
|
#if DECUSE64
|
||
|
|
uLong mul;
|
||
|
|
// as before, bring in extra digits if possible
|
||
|
|
if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
|
||
|
|
mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
|
||
|
|
+ *(msua-2)/DIVLO;
|
||
|
|
mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
|
||
|
|
}
|
||
|
|
else if (divunits==1) {
|
||
|
|
mul=(uLong)*msua * DIVBASE + *(msua-1);
|
||
|
|
mul/=*msud; // no more to the right
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
mul=(uLong)(*msua) * (uInt)(DIVBASE<<2)
|
||
|
|
+ (*(msua-1)<<2);
|
||
|
|
mul/=divtop; // [divtop already allows for sticky bits]
|
||
|
|
}
|
||
|
|
multiplier=(Int)mul;
|
||
|
|
#else
|
||
|
|
multiplier=*msua * ((DIVBASE<<2)/divtop);
|
||
|
|
#endif
|
||
|
|
}
|
||
|
|
if (multiplier==0) multiplier=1; // marginal case
|
||
|
|
*lsuq+=multiplier;
|
||
|
|
|
||
|
|
#if DIVCOUNT
|
||
|
|
// printf("Multiplier: %ld\n", (LI)multiplier);
|
||
|
|
divcount++;
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// Carry out the subtraction acc-(div*multiplier); for each
|
||
|
|
// unit in div, do the multiply, split to units (see
|
||
|
|
// decFloatMultiply for the algorithm), and subtract from acc
|
||
|
|
#define DIVMAGIC 2305843009U // 2**61/10**9
|
||
|
|
#define DIVSHIFTA 29
|
||
|
|
#define DIVSHIFTB 32
|
||
|
|
carry=0;
|
||
|
|
for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
|
||
|
|
uInt lo, hop;
|
||
|
|
#if DECUSE64
|
||
|
|
uLong sub=(uLong)multiplier*(*ud)+carry;
|
||
|
|
if (sub<DIVBASE) {
|
||
|
|
carry=0;
|
||
|
|
lo=(uInt)sub;
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
hop=(uInt)(sub>>DIVSHIFTA);
|
||
|
|
carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
|
||
|
|
// the estimate is now in hi; now calculate sub-hi*10**9
|
||
|
|
// to get the remainder (which will be <DIVBASE))
|
||
|
|
lo=(uInt)sub;
|
||
|
|
lo-=carry*DIVBASE; // low word of result
|
||
|
|
if (lo>=DIVBASE) {
|
||
|
|
lo-=DIVBASE; // correct by +1
|
||
|
|
carry++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
#else // 32-bit
|
||
|
|
uInt hi;
|
||
|
|
// calculate multiplier*(*ud) into hi and lo
|
||
|
|
LONGMUL32HI(hi, *ud, multiplier); // get the high word
|
||
|
|
lo=multiplier*(*ud); // .. and the low
|
||
|
|
lo+=carry; // add the old hi
|
||
|
|
carry=hi+(lo<carry); // .. with any carry
|
||
|
|
if (carry || lo>=DIVBASE) { // split is needed
|
||
|
|
hop=(carry<<3)+(lo>>DIVSHIFTA); // hi:lo/2**29
|
||
|
|
LONGMUL32HI(carry, hop, DIVMAGIC); // only need the high word
|
||
|
|
// [DIVSHIFTB is 32, so carry can be used directly]
|
||
|
|
// the estimate is now in carry; now calculate hi:lo-est*10**9;
|
||
|
|
// happily the top word of the result is irrelevant because it
|
||
|
|
// will always be zero so this needs only one multiplication
|
||
|
|
lo-=(carry*DIVBASE);
|
||
|
|
// the correction here will be at most +1; do it
|
||
|
|
if (lo>=DIVBASE) {
|
||
|
|
lo-=DIVBASE;
|
||
|
|
carry++;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
if (lo>*ua) { // borrow needed
|
||
|
|
*ua+=DIVBASE;
|
||
|
|
carry++;
|
||
|
|
}
|
||
|
|
*ua-=lo;
|
||
|
|
} // ud loop
|
||
|
|
if (carry) *ua-=carry; // accdigits>divdigits [cannot borrow]
|
||
|
|
} // inner loop
|
||
|
|
|
||
|
|
// the outer loop terminates when there is either an exact result
|
||
|
|
// or enough digits; first update the quotient digit count and
|
||
|
|
// pointer (if any significant digits)
|
||
|
|
#if DECTRACE
|
||
|
|
if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
|
||
|
|
#endif
|
||
|
|
if (quodigits) {
|
||
|
|
quodigits+=9; // had leading unit earlier
|
||
|
|
lsuq--;
|
||
|
|
if (quodigits>DECPMAX+1) break; // have enough
|
||
|
|
}
|
||
|
|
else if (*lsuq) { // first quotient digits
|
||
|
|
const uInt *pow;
|
||
|
|
for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
|
||
|
|
lsuq--;
|
||
|
|
// [cannot have >DECPMAX+1 on first unit]
|
||
|
|
}
|
||
|
|
|
||
|
|
if (*msua!=0) continue; // not an exact result
|
||
|
|
// acc is zero iff used all of original units and zero down to lsua
|
||
|
|
// (must also continue to original lsu for correct quotient length)
|
||
|
|
if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
|
||
|
|
for (; msua>lsua && *msua==0;) msua--;
|
||
|
|
if (*msua==0 && msua==lsua) break;
|
||
|
|
} // outer loop
|
||
|
|
|
||
|
|
// all of the original operand in acc has been covered at this point
|
||
|
|
// quotient now has at least DECPMAX+2 digits
|
||
|
|
// *msua is now non-0 if inexact and sticky bits
|
||
|
|
// lsuq is one below the last uint of the quotient
|
||
|
|
lsuq++; // set -> true lsu of quo
|
||
|
|
if (*msua) *lsuq|=1; // apply sticky bit
|
||
|
|
|
||
|
|
// quo now holds the (unrounded) quotient in base-billion; one
|
||
|
|
// base-billion 'digit' per uInt.
|
||
|
|
#if DECTRACE
|
||
|
|
printf("DivQuo:");
|
||
|
|
for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
|
||
|
|
printf("\n");
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// Now convert to BCD for rounding and cleanup, starting from the
|
||
|
|
// most significant end [offset by one into bcdacc to leave room
|
||
|
|
// for a possible carry digit if rounding for REMNEAR is needed]
|
||
|
|
for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
|
||
|
|
uInt top, mid, rem; // work
|
||
|
|
if (*uq==0) { // no split needed
|
||
|
|
UBFROMUI(ub, 0); // clear 9 BCD8s
|
||
|
|
UBFROMUI(ub+4, 0); // ..
|
||
|
|
*(ub+8)=0; // ..
|
||
|
|
continue;
|
||
|
|
}
|
||
|
|
// *uq is non-zero -- split the base-billion digit into
|
||
|
|
// hi, mid, and low three-digits
|
||
|
|
#define divsplit9 1000000 // divisor
|
||
|
|
#define divsplit6 1000 // divisor
|
||
|
|
// The splitting is done by simple divides and remainders,
|
||
|
|
// assuming the compiler will optimize these [GCC does]
|
||
|
|
top=*uq/divsplit9;
|
||
|
|
rem=*uq%divsplit9;
|
||
|
|
mid=rem/divsplit6;
|
||
|
|
rem=rem%divsplit6;
|
||
|
|
// lay out the nine BCD digits (plus one unwanted byte)
|
||
|
|
UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
|
||
|
|
UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
|
||
|
|
UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
|
||
|
|
} // BCD conversion loop
|
||
|
|
ub--; // -> lsu
|
||
|
|
|
||
|
|
// complete the bcdnum; quodigits is correct, so the position of
|
||
|
|
// the first non-zero is known
|
||
|
|
num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
|
||
|
|
num.lsd=ub;
|
||
|
|
|
||
|
|
// make exponent adjustments, etc
|
||
|
|
if (lsua<acc+DIVACCLEN-DIVOPLEN) { // used extra digits
|
||
|
|
num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
|
||
|
|
// if the result was exact then there may be up to 8 extra
|
||
|
|
// trailing zeros in the overflowed quotient final unit
|
||
|
|
if (*msua==0) {
|
||
|
|
for (; *ub==0;) ub--; // drop zeros
|
||
|
|
num.exponent+=(Int)(num.lsd-ub); // and adjust exponent
|
||
|
|
num.lsd=ub;
|
||
|
|
}
|
||
|
|
} // adjustment needed
|
||
|
|
|
||
|
|
#if DIVCOUNT
|
||
|
|
if (divcount>maxcount) { // new high-water nark
|
||
|
|
maxcount=divcount;
|
||
|
|
printf("DivNewMaxCount: %ld\n", (LI)maxcount);
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
|
||
|
|
if (op&DIVIDE) return decFinalize(result, &num, set); // all done
|
||
|
|
|
||
|
|
// Is DIVIDEINT or a remainder; there is more to do -- first form
|
||
|
|
// the integer (this is done 'after the fact', unlike as in
|
||
|
|
// decNumber, so as not to tax DIVIDE)
|
||
|
|
|
||
|
|
// The first non-zero digit will be in the first 9 digits, known
|
||
|
|
// from quodigits and num.msd, so there is always space for DECPMAX
|
||
|
|
// digits
|
||
|
|
|
||
|
|
length=(Int)(num.lsd-num.msd+1);
|
||
|
|
//printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent);
|
||
|
|
|
||
|
|
if (length+num.exponent>DECPMAX) { // cannot fit
|
||
|
|
decFloatZero(result);
|
||
|
|
DFWORD(result, 0)=DECFLOAT_qNaN;
|
||
|
|
set->status|=DEC_Division_impossible;
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
|
||
|
|
if (num.exponent>=0) { // already an int, or need pad zeros
|
||
|
|
for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
|
||
|
|
num.lsd+=num.exponent;
|
||
|
|
}
|
||
|
|
else { // too long: round or truncate needed
|
||
|
|
Int drop=-num.exponent;
|
||
|
|
if (!(op&REMNEAR)) { // simple truncate
|
||
|
|
num.lsd-=drop;
|
||
|
|
if (num.lsd<num.msd) { // truncated all
|
||
|
|
num.lsd=num.msd; // make 0
|
||
|
|
*num.lsd=0; // .. [sign still relevant]
|
||
|
|
}
|
||
|
|
}
|
||
|
|
else { // round to nearest even [sigh]
|
||
|
|
// round-to-nearest, in-place; msd is at or to right of bcdacc+1
|
||
|
|
// (this is a special case of Quantize -- q.v. for commentary)
|
||
|
|
uByte *roundat; // -> re-round digit
|
||
|
|
uByte reround; // reround value
|
||
|
|
*(num.msd-1)=0; // in case of left carry, or make 0
|
||
|
|
if (drop<length) roundat=num.lsd-drop+1;
|
||
|
|
else if (drop==length) roundat=num.msd;
|
||
|
|
else roundat=num.msd-1; // [-> 0]
|
||
|
|
reround=*roundat;
|
||
|
|
for (ub=roundat+1; ub<=num.lsd; ub++) {
|
||
|
|
if (*ub!=0) {
|
||
|
|
reround=DECSTICKYTAB[reround];
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
} // check stickies
|
||
|
|
if (roundat>num.msd) num.lsd=roundat-1;
|
||
|
|
else {
|
||
|
|
num.msd--; // use the 0 ..
|
||
|
|
num.lsd=num.msd; // .. at the new MSD place
|
||
|
|
}
|
||
|
|
if (reround!=0) { // discarding non-zero
|
||
|
|
uInt bump=0;
|
||
|
|
// rounding is DEC_ROUND_HALF_EVEN always
|
||
|
|
if (reround>5) bump=1; // >0.5 goes up
|
||
|
|
else if (reround==5) // exactly 0.5000 ..
|
||
|
|
bump=*(num.lsd) & 0x01; // .. up iff [new] lsd is odd
|
||
|
|
if (bump!=0) { // need increment
|
||
|
|
// increment the coefficient; this might end up with 1000...
|
||
|
|
ub=num.lsd;
|
||
|
|
for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
|
||
|
|
for (; *ub==9; ub--) *ub=0; // at most 3 more
|
||
|
|
*ub+=1;
|
||
|
|
if (ub<num.msd) num.msd--; // carried
|
||
|
|
} // bump needed
|
||
|
|
} // reround!=0
|
||
|
|
} // remnear
|
||
|
|
} // round or truncate needed
|
||
|
|
num.exponent=0; // all paths
|
||
|
|
//decShowNum(&num, "int");
|
||
|
|
|
||
|
|
if (op&DIVIDEINT) return decFinalize(result, &num, set); // all done
|
||
|
|
|
||
|
|
// Have a remainder to calculate
|
||
|
|
decFinalize("ient, &num, set); // lay out the integer so far
|
||
|
|
DFWORD("ient, 0)^=DECFLOAT_Sign; // negate it
|
||
|
|
sign=DFWORD(dfl, 0); // save sign of dfl
|
||
|
|
decFloatFMA(result, "ient, dfr, dfl, set);
|
||
|
|
if (!DFISZERO(result)) return result;
|
||
|
|
// if the result is zero the sign shall be sign of dfl
|
||
|
|
DFWORD("ient, 0)=sign; // construct decFloat of sign
|
||
|
|
return decFloatCopySign(result, result, "ient);
|
||
|
|
} // decDivide
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFiniteMultiply -- multiply two finite decFloats */
|
||
|
|
/* */
|
||
|
|
/* num gets the result of multiplying dfl and dfr */
|
||
|
|
/* bcdacc .. with the coefficient in this array */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* */
|
||
|
|
/* This effects the multiplication of two decFloats, both known to be */
|
||
|
|
/* finite, leaving the result in a bcdnum ready for decFinalize (for */
|
||
|
|
/* use in Multiply) or in a following addition (FMA). */
|
||
|
|
/* */
|
||
|
|
/* bcdacc must have space for at least DECPMAX9*18+1 bytes. */
|
||
|
|
/* No error is possible and no status is set. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
// This routine has two separate implementations of the core
|
||
|
|
// multiplication; both using base-billion. One uses only 32-bit
|
||
|
|
// variables (Ints and uInts) or smaller; the other uses uLongs (for
|
||
|
|
// multiplication and addition only). Both implementations cover
|
||
|
|
// both arithmetic sizes (DOUBLE and QUAD) in order to allow timing
|
||
|
|
// comparisons. In any one compilation only one implementation for
|
||
|
|
// each size can be used, and if DECUSE64 is 0 then use of the 32-bit
|
||
|
|
// version is forced.
|
||
|
|
//
|
||
|
|
// Historical note: an earlier version of this code also supported the
|
||
|
|
// 256-bit format and has been preserved. That is somewhat trickier
|
||
|
|
// during lazy carry splitting because the initial quotient estimate
|
||
|
|
// (est) can exceed 32 bits.
|
||
|
|
|
||
|
|
#define MULTBASE ((uInt)BILLION) // the base used for multiply
|
||
|
|
#define MULOPLEN DECPMAX9 // operand length ('digits' base 10**9)
|
||
|
|
#define MULACCLEN (MULOPLEN*2) // accumulator length (ditto)
|
||
|
|
#define LEADZEROS (MULACCLEN*9 - DECPMAX*2) // leading zeros always
|
||
|
|
|
||
|
|
// Assertions: exponent not too large and MULACCLEN is a multiple of 4
|
||
|
|
#if DECEMAXD>9
|
||
|
|
#error Exponent may overflow when doubled for Multiply
|
||
|
|
#endif
|
||
|
|
#if MULACCLEN!=(MULACCLEN/4)*4
|
||
|
|
// This assumption is used below only for initialization
|
||
|
|
#error MULACCLEN is not a multiple of 4
|
||
|
|
#endif
|
||
|
|
|
||
|
|
static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
|
||
|
|
const decFloat *dfl, const decFloat *dfr) {
|
||
|
|
uInt bufl[MULOPLEN]; // left coefficient (base-billion)
|
||
|
|
uInt bufr[MULOPLEN]; // right coefficient (base-billion)
|
||
|
|
uInt *ui, *uj; // work
|
||
|
|
uByte *ub; // ..
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
|
||
|
|
#if DECUSE64
|
||
|
|
uLong accl[MULACCLEN]; // lazy accumulator (base-billion+)
|
||
|
|
uLong *pl; // work -> lazy accumulator
|
||
|
|
uInt acc[MULACCLEN]; // coefficent in base-billion ..
|
||
|
|
#else
|
||
|
|
uInt acc[MULACCLEN*2]; // accumulator in base-billion ..
|
||
|
|
#endif
|
||
|
|
uInt *pa; // work -> accumulator
|
||
|
|
//printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN);
|
||
|
|
|
||
|
|
/* Calculate sign and exponent */
|
||
|
|
num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
|
||
|
|
num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); // [see assertion above]
|
||
|
|
|
||
|
|
/* Extract the coefficients and prepare the accumulator */
|
||
|
|
// the coefficients of the operands are decoded into base-billion
|
||
|
|
// numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the
|
||
|
|
// appropriate size.
|
||
|
|
GETCOEFFBILL(dfl, bufl);
|
||
|
|
GETCOEFFBILL(dfr, bufr);
|
||
|
|
#if DECTRACE && 0
|
||
|
|
printf("CoeffbL:");
|
||
|
|
for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
|
||
|
|
printf("\n");
|
||
|
|
printf("CoeffbR:");
|
||
|
|
for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
|
||
|
|
printf("\n");
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// start the 64-bit/32-bit differing paths...
|
||
|
|
#if DECUSE64
|
||
|
|
|
||
|
|
// zero the accumulator
|
||
|
|
#if MULACCLEN==4
|
||
|
|
accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
|
||
|
|
#else // use a loop
|
||
|
|
// MULACCLEN is a multiple of four, asserted above
|
||
|
|
for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
|
||
|
|
*pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;// [reduce overhead]
|
||
|
|
} // pl
|
||
|
|
#endif
|
||
|
|
|
||
|
|
/* Effect the multiplication */
|
||
|
|
// The multiplcation proceeds using MFC's lazy-carry resolution
|
||
|
|
// algorithm from decNumber. First, the multiplication is
|
||
|
|
// effected, allowing accumulation of the partial products (which
|
||
|
|
// are in base-billion at each column position) into 64 bits
|
||
|
|
// without resolving back to base=billion after each addition.
|
||
|
|
// These 64-bit numbers (which may contain up to 19 decimal digits)
|
||
|
|
// are then split using the Clark & Cowlishaw algorithm (see below).
|
||
|
|
// [Testing for 0 in the inner loop is not really a 'win']
|
||
|
|
for (ui=bufr; ui<bufr+MULOPLEN; ui++) { // over each item in rhs
|
||
|
|
if (*ui==0) continue; // product cannot affect result
|
||
|
|
pl=accl+(ui-bufr); // where to add the lhs
|
||
|
|
for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { // over each item in lhs
|
||
|
|
// if (*uj==0) continue; // product cannot affect result
|
||
|
|
*pl+=((uLong)*ui)*(*uj);
|
||
|
|
} // uj
|
||
|
|
} // ui
|
||
|
|
|
||
|
|
// The 64-bit carries must now be resolved; this means that a
|
||
|
|
// quotient/remainder has to be calculated for base-billion (1E+9).
|
||
|
|
// For this, Clark & Cowlishaw's quotient estimation approach (also
|
||
|
|
// used in decNumber) is needed, because 64-bit divide is generally
|
||
|
|
// extremely slow on 32-bit machines, and may be slower than this
|
||
|
|
// approach even on 64-bit machines. This algorithm splits X
|
||
|
|
// using:
|
||
|
|
//
|
||
|
|
// magic=2**(A+B)/1E+9; // 'magic number'
|
||
|
|
// hop=X/2**A; // high order part of X (by shift)
|
||
|
|
// est=magic*hop/2**B // quotient estimate (may be low by 1)
|
||
|
|
//
|
||
|
|
// A and B are quite constrained; hop and magic must fit in 32 bits,
|
||
|
|
// and 2**(A+B) must be as large as possible (which is 2**61 if
|
||
|
|
// magic is to fit). Further, maxX increases with the length of
|
||
|
|
// the operands (and hence the number of partial products
|
||
|
|
// accumulated); maxX is OPLEN*(10**18), which is up to 19 digits.
|
||
|
|
//
|
||
|
|
// It can be shown that when OPLEN is 2 then the maximum error in
|
||
|
|
// the estimated quotient is <1, but for larger maximum x the
|
||
|
|
// maximum error is above 1 so a correction that is >1 may be
|
||
|
|
// needed. Values of A and B are chosen to satisfy the constraints
|
||
|
|
// just mentioned while minimizing the maximum error (and hence the
|
||
|
|
// maximum correction), as shown in the following table:
|
||
|
|
//
|
||
|
|
// Type OPLEN A B maxX maxError maxCorrection
|
||
|
|
// ---------------------------------------------------------
|
||
|
|
// DOUBLE 2 29 32 <2*10**18 0.63 1
|
||
|
|
// QUAD 4 30 31 <4*10**18 1.17 2
|
||
|
|
//
|
||
|
|
// In the OPLEN==2 case there is most choice, but the value for B
|
||
|
|
// of 32 has a big advantage as then the calculation of the
|
||
|
|
// estimate requires no shifting; the compiler can extract the high
|
||
|
|
// word directly after multiplying magic*hop.
|
||
|
|
#define MULMAGIC 2305843009U // 2**61/10**9 [both cases]
|
||
|
|
#if DOUBLE
|
||
|
|
#define MULSHIFTA 29
|
||
|
|
#define MULSHIFTB 32
|
||
|
|
#elif QUAD
|
||
|
|
#define MULSHIFTA 30
|
||
|
|
#define MULSHIFTB 31
|
||
|
|
#else
|
||
|
|
#error Unexpected type
|
||
|
|
#endif
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
printf("MulAccl:");
|
||
|
|
for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
|
||
|
|
printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
|
||
|
|
printf("\n");
|
||
|
|
#endif
|
||
|
|
|
||
|
|
for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { // each column position
|
||
|
|
uInt lo, hop; // work
|
||
|
|
uInt est; // cannot exceed 4E+9
|
||
|
|
if (*pl>=MULTBASE) {
|
||
|
|
// *pl holds a binary number which needs to be split
|
||
|
|
hop=(uInt)(*pl>>MULSHIFTA);
|
||
|
|
est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
|
||
|
|
// the estimate is now in est; now calculate hi:lo-est*10**9;
|
||
|
|
// happily the top word of the result is irrelevant because it
|
||
|
|
// will always be zero so this needs only one multiplication
|
||
|
|
lo=(uInt)(*pl-((uLong)est*MULTBASE)); // low word of result
|
||
|
|
// If QUAD, the correction here could be +2
|
||
|
|
if (lo>=MULTBASE) {
|
||
|
|
lo-=MULTBASE; // correct by +1
|
||
|
|
est++;
|
||
|
|
#if QUAD
|
||
|
|
// may need to correct by +2
|
||
|
|
if (lo>=MULTBASE) {
|
||
|
|
lo-=MULTBASE;
|
||
|
|
est++;
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
}
|
||
|
|
// finally place lo as the new coefficient 'digit' and add est to
|
||
|
|
// the next place up [this is safe because this path is never
|
||
|
|
// taken on the final iteration as *pl will fit]
|
||
|
|
*pa=lo;
|
||
|
|
*(pl+1)+=est;
|
||
|
|
} // *pl needed split
|
||
|
|
else { // *pl<MULTBASE
|
||
|
|
*pa=(uInt)*pl; // just copy across
|
||
|
|
}
|
||
|
|
} // pl loop
|
||
|
|
|
||
|
|
#else // 32-bit
|
||
|
|
for (pa=acc;; pa+=4) { // zero the accumulator
|
||
|
|
*pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; // [reduce overhead]
|
||
|
|
if (pa==acc+MULACCLEN*2-4) break; // multiple of 4 asserted
|
||
|
|
} // pa
|
||
|
|
|
||
|
|
/* Effect the multiplication */
|
||
|
|
// uLongs are not available (and in particular, there is no uLong
|
||
|
|
// divide) but it is still possible to use MFC's lazy-carry
|
||
|
|
// resolution algorithm from decNumber. First, the multiplication
|
||
|
|
// is effected, allowing accumulation of the partial products
|
||
|
|
// (which are in base-billion at each column position) into 64 bits
|
||
|
|
// [with the high-order 32 bits in each position being held at
|
||
|
|
// offset +ACCLEN from the low-order 32 bits in the accumulator].
|
||
|
|
// These 64-bit numbers (which may contain up to 19 decimal digits)
|
||
|
|
// are then split using the Clark & Cowlishaw algorithm (see
|
||
|
|
// below).
|
||
|
|
for (ui=bufr;; ui++) { // over each item in rhs
|
||
|
|
uInt hi, lo; // words of exact multiply result
|
||
|
|
pa=acc+(ui-bufr); // where to add the lhs
|
||
|
|
for (uj=bufl;; uj++, pa++) { // over each item in lhs
|
||
|
|
LONGMUL32HI(hi, *ui, *uj); // calculate product of digits
|
||
|
|
lo=(*ui)*(*uj); // ..
|
||
|
|
*pa+=lo; // accumulate low bits and ..
|
||
|
|
*(pa+MULACCLEN)+=hi+(*pa<lo); // .. high bits with any carry
|
||
|
|
if (uj==bufl+MULOPLEN-1) break;
|
||
|
|
}
|
||
|
|
if (ui==bufr+MULOPLEN-1) break;
|
||
|
|
}
|
||
|
|
|
||
|
|
// The 64-bit carries must now be resolved; this means that a
|
||
|
|
// quotient/remainder has to be calculated for base-billion (1E+9).
|
||
|
|
// For this, Clark & Cowlishaw's quotient estimation approach (also
|
||
|
|
// used in decNumber) is needed, because 64-bit divide is generally
|
||
|
|
// extremely slow on 32-bit machines. This algorithm splits X
|
||
|
|
// using:
|
||
|
|
//
|
||
|
|
// magic=2**(A+B)/1E+9; // 'magic number'
|
||
|
|
// hop=X/2**A; // high order part of X (by shift)
|
||
|
|
// est=magic*hop/2**B // quotient estimate (may be low by 1)
|
||
|
|
//
|
||
|
|
// A and B are quite constrained; hop and magic must fit in 32 bits,
|
||
|
|
// and 2**(A+B) must be as large as possible (which is 2**61 if
|
||
|
|
// magic is to fit). Further, maxX increases with the length of
|
||
|
|
// the operands (and hence the number of partial products
|
||
|
|
// accumulated); maxX is OPLEN*(10**18), which is up to 19 digits.
|
||
|
|
//
|
||
|
|
// It can be shown that when OPLEN is 2 then the maximum error in
|
||
|
|
// the estimated quotient is <1, but for larger maximum x the
|
||
|
|
// maximum error is above 1 so a correction that is >1 may be
|
||
|
|
// needed. Values of A and B are chosen to satisfy the constraints
|
||
|
|
// just mentioned while minimizing the maximum error (and hence the
|
||
|
|
// maximum correction), as shown in the following table:
|
||
|
|
//
|
||
|
|
// Type OPLEN A B maxX maxError maxCorrection
|
||
|
|
// ---------------------------------------------------------
|
||
|
|
// DOUBLE 2 29 32 <2*10**18 0.63 1
|
||
|
|
// QUAD 4 30 31 <4*10**18 1.17 2
|
||
|
|
//
|
||
|
|
// In the OPLEN==2 case there is most choice, but the value for B
|
||
|
|
// of 32 has a big advantage as then the calculation of the
|
||
|
|
// estimate requires no shifting; the high word is simply
|
||
|
|
// calculated from multiplying magic*hop.
|
||
|
|
#define MULMAGIC 2305843009U // 2**61/10**9 [both cases]
|
||
|
|
#if DOUBLE
|
||
|
|
#define MULSHIFTA 29
|
||
|
|
#define MULSHIFTB 32
|
||
|
|
#elif QUAD
|
||
|
|
#define MULSHIFTA 30
|
||
|
|
#define MULSHIFTB 31
|
||
|
|
#else
|
||
|
|
#error Unexpected type
|
||
|
|
#endif
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
printf("MulHiLo:");
|
||
|
|
for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
|
||
|
|
printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
|
||
|
|
printf("\n");
|
||
|
|
#endif
|
||
|
|
|
||
|
|
for (pa=acc;; pa++) { // each low uInt
|
||
|
|
uInt hi, lo; // words of exact multiply result
|
||
|
|
uInt hop, estlo; // work
|
||
|
|
#if QUAD
|
||
|
|
uInt esthi; // ..
|
||
|
|
#endif
|
||
|
|
|
||
|
|
lo=*pa;
|
||
|
|
hi=*(pa+MULACCLEN); // top 32 bits
|
||
|
|
// hi and lo now hold a binary number which needs to be split
|
||
|
|
|
||
|
|
#if DOUBLE
|
||
|
|
hop=(hi<<3)+(lo>>MULSHIFTA); // hi:lo/2**29
|
||
|
|
LONGMUL32HI(estlo, hop, MULMAGIC);// only need the high word
|
||
|
|
// [MULSHIFTB is 32, so estlo can be used directly]
|
||
|
|
// the estimate is now in estlo; now calculate hi:lo-est*10**9;
|
||
|
|
// happily the top word of the result is irrelevant because it
|
||
|
|
// will always be zero so this needs only one multiplication
|
||
|
|
lo-=(estlo*MULTBASE);
|
||
|
|
// esthi=0; // high word is ignored below
|
||
|
|
// the correction here will be at most +1; do it
|
||
|
|
if (lo>=MULTBASE) {
|
||
|
|
lo-=MULTBASE;
|
||
|
|
estlo++;
|
||
|
|
}
|
||
|
|
#elif QUAD
|
||
|
|
hop=(hi<<2)+(lo>>MULSHIFTA); // hi:lo/2**30
|
||
|
|
LONGMUL32HI(esthi, hop, MULMAGIC);// shift will be 31 ..
|
||
|
|
estlo=hop*MULMAGIC; // .. so low word needed
|
||
|
|
estlo=(esthi<<1)+(estlo>>MULSHIFTB); // [just the top bit]
|
||
|
|
// esthi=0; // high word is ignored below
|
||
|
|
lo-=(estlo*MULTBASE); // as above
|
||
|
|
// the correction here could be +1 or +2
|
||
|
|
if (lo>=MULTBASE) {
|
||
|
|
lo-=MULTBASE;
|
||
|
|
estlo++;
|
||
|
|
}
|
||
|
|
if (lo>=MULTBASE) {
|
||
|
|
lo-=MULTBASE;
|
||
|
|
estlo++;
|
||
|
|
}
|
||
|
|
#else
|
||
|
|
#error Unexpected type
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// finally place lo as the new accumulator digit and add est to
|
||
|
|
// the next place up; this latter add could cause a carry of 1
|
||
|
|
// to the high word of the next place
|
||
|
|
*pa=lo;
|
||
|
|
*(pa+1)+=estlo;
|
||
|
|
// esthi is always 0 for DOUBLE and QUAD so this is skipped
|
||
|
|
// *(pa+1+MULACCLEN)+=esthi;
|
||
|
|
if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; // carry
|
||
|
|
if (pa==acc+MULACCLEN-2) break; // [MULACCLEN-1 will never need split]
|
||
|
|
} // pa loop
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// At this point, whether using the 64-bit or the 32-bit paths, the
|
||
|
|
// accumulator now holds the (unrounded) result in base-billion;
|
||
|
|
// one base-billion 'digit' per uInt.
|
||
|
|
#if DECTRACE
|
||
|
|
printf("MultAcc:");
|
||
|
|
for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
|
||
|
|
printf("\n");
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// Now convert to BCD for rounding and cleanup, starting from the
|
||
|
|
// most significant end
|
||
|
|
pa=acc+MULACCLEN-1;
|
||
|
|
if (*pa!=0) num->msd=bcdacc+LEADZEROS;// drop known lead zeros
|
||
|
|
else { // >=1 word of leading zeros
|
||
|
|
num->msd=bcdacc; // known leading zeros are gone
|
||
|
|
pa--; // skip first word ..
|
||
|
|
for (; *pa==0; pa--) if (pa==acc) break; // .. and any more leading 0s
|
||
|
|
}
|
||
|
|
for (ub=bcdacc;; pa--, ub+=9) {
|
||
|
|
if (*pa!=0) { // split(s) needed
|
||
|
|
uInt top, mid, rem; // work
|
||
|
|
// *pa is non-zero -- split the base-billion acc digit into
|
||
|
|
// hi, mid, and low three-digits
|
||
|
|
#define mulsplit9 1000000 // divisor
|
||
|
|
#define mulsplit6 1000 // divisor
|
||
|
|
// The splitting is done by simple divides and remainders,
|
||
|
|
// assuming the compiler will optimize these where useful
|
||
|
|
// [GCC does]
|
||
|
|
top=*pa/mulsplit9;
|
||
|
|
rem=*pa%mulsplit9;
|
||
|
|
mid=rem/mulsplit6;
|
||
|
|
rem=rem%mulsplit6;
|
||
|
|
// lay out the nine BCD digits (plus one unwanted byte)
|
||
|
|
UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
|
||
|
|
UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
|
||
|
|
UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
|
||
|
|
}
|
||
|
|
else { // *pa==0
|
||
|
|
UBFROMUI(ub, 0); // clear 9 BCD8s
|
||
|
|
UBFROMUI(ub+4, 0); // ..
|
||
|
|
*(ub+8)=0; // ..
|
||
|
|
}
|
||
|
|
if (pa==acc) break;
|
||
|
|
} // BCD conversion loop
|
||
|
|
|
||
|
|
num->lsd=ub+8; // complete the bcdnum ..
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
decShowNum(num, "postmult");
|
||
|
|
decFloatShow(dfl, "dfl");
|
||
|
|
decFloatShow(dfr, "dfr");
|
||
|
|
#endif
|
||
|
|
return;
|
||
|
|
} // decFiniteMultiply
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatAbs -- absolute value, heeding NaNs, etc. */
|
||
|
|
/* */
|
||
|
|
/* result gets the canonicalized df with sign 0 */
|
||
|
|
/* df is the decFloat to abs */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This has the same effect as decFloatPlus unless df is negative, */
|
||
|
|
/* in which case it has the same effect as decFloatMinus. The */
|
||
|
|
/* effect is also the same as decFloatCopyAbs except that NaNs are */
|
||
|
|
/* handled normally (the sign of a NaN is not affected, and an sNaN */
|
||
|
|
/* will signal) and the result will be canonical. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatAbs(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set) {
|
||
|
|
if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
|
||
|
|
decCanonical(result, df); // copy and check
|
||
|
|
DFBYTE(result, 0)&=~0x80; // zero sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatAbs
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatAdd -- add two decFloats */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of adding dfl and dfr: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
#if QUAD
|
||
|
|
// Table for testing MSDs for fastpath elimination; returns the MSD of
|
||
|
|
// a decDouble or decQuad (top 6 bits tested) ignoring the sign.
|
||
|
|
// Infinities return -32 and NaNs return -128 so that summing the two
|
||
|
|
// MSDs also allows rapid tests for the Specials (see code below).
|
||
|
|
const Int DECTESTMSD[64]={
|
||
|
|
0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
|
||
|
|
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128,
|
||
|
|
0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
|
||
|
|
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128};
|
||
|
|
#else
|
||
|
|
// The table for testing MSDs is shared between the modules
|
||
|
|
extern const Int DECTESTMSD[64];
|
||
|
|
#endif
|
||
|
|
|
||
|
|
decFloat * decFloatAdd(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
bcdnum num; // for final conversion
|
||
|
|
Int bexpl, bexpr; // left and right biased exponents
|
||
|
|
uByte *ub, *us, *ut; // work
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
#if QUAD
|
||
|
|
uShort uswork; // ..
|
||
|
|
#endif
|
||
|
|
|
||
|
|
uInt sourhil, sourhir; // top words from source decFloats
|
||
|
|
// [valid only through end of
|
||
|
|
// fastpath code -- before swap]
|
||
|
|
uInt diffsign; // non-zero if signs differ
|
||
|
|
uInt carry; // carry: 0 or 1 before add loop
|
||
|
|
Int overlap; // coefficient overlap (if full)
|
||
|
|
Int summ; // sum of the MSDs
|
||
|
|
// the following buffers hold coefficients with various alignments
|
||
|
|
// (see commentary and diagrams below)
|
||
|
|
uByte acc[4+2+DECPMAX*3+8];
|
||
|
|
uByte buf[4+2+DECPMAX*2];
|
||
|
|
uByte *umsd, *ulsd; // local MSD and LSD pointers
|
||
|
|
|
||
|
|
#if DECLITEND
|
||
|
|
#define CARRYPAT 0x01000000 // carry=1 pattern
|
||
|
|
#else
|
||
|
|
#define CARRYPAT 0x00000001 // carry=1 pattern
|
||
|
|
#endif
|
||
|
|
|
||
|
|
/* Start decoding the arguments */
|
||
|
|
// The initial exponents are placed into the opposite Ints to
|
||
|
|
// that which might be expected; there are two sets of data to
|
||
|
|
// keep track of (each decFloat and the corresponding exponent),
|
||
|
|
// and this scheme means that at the swap point (after comparing
|
||
|
|
// exponents) only one pair of words needs to be swapped
|
||
|
|
// whichever path is taken (thereby minimising worst-case path).
|
||
|
|
// The calculated exponents will be nonsense when the arguments are
|
||
|
|
// Special, but are not used in that path
|
||
|
|
sourhil=DFWORD(dfl, 0); // LHS top word
|
||
|
|
summ=DECTESTMSD[sourhil>>26]; // get first MSD for testing
|
||
|
|
bexpr=DECCOMBEXP[sourhil>>26]; // get exponent high bits (in place)
|
||
|
|
bexpr+=GETECON(dfl); // .. + continuation
|
||
|
|
|
||
|
|
sourhir=DFWORD(dfr, 0); // RHS top word
|
||
|
|
summ+=DECTESTMSD[sourhir>>26]; // sum MSDs for testing
|
||
|
|
bexpl=DECCOMBEXP[sourhir>>26];
|
||
|
|
bexpl+=GETECON(dfr);
|
||
|
|
|
||
|
|
// here bexpr has biased exponent from lhs, and vice versa
|
||
|
|
|
||
|
|
diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
|
||
|
|
|
||
|
|
// now determine whether to take a fast path or the full-function
|
||
|
|
// slow path. The slow path must be taken when:
|
||
|
|
// -- both numbers are finite, and:
|
||
|
|
// the exponents are different, or
|
||
|
|
// the signs are different, or
|
||
|
|
// the sum of the MSDs is >8 (hence might overflow)
|
||
|
|
// specialness and the sum of the MSDs can be tested at once using
|
||
|
|
// the summ value just calculated, so the test for specials is no
|
||
|
|
// longer on the worst-case path (as of 3.60)
|
||
|
|
|
||
|
|
if (summ<=8) { // MSD+MSD is good, or there is a special
|
||
|
|
if (summ<0) { // there is a special
|
||
|
|
// Inf+Inf would give -64; Inf+finite is -32 or higher
|
||
|
|
if (summ<-64) return decNaNs(result, dfl, dfr, set); // one or two NaNs
|
||
|
|
// two infinities with different signs is invalid
|
||
|
|
if (summ==-64 && diffsign) return decInvalid(result, set);
|
||
|
|
if (DFISINF(dfl)) return decInfinity(result, dfl); // LHS is infinite
|
||
|
|
return decInfinity(result, dfr); // RHS must be Inf
|
||
|
|
}
|
||
|
|
// Here when both arguments are finite; fast path is possible
|
||
|
|
// (currently only for aligned and same-sign)
|
||
|
|
if (bexpr==bexpl && !diffsign) {
|
||
|
|
uInt tac[DECLETS+1]; // base-1000 coefficient
|
||
|
|
uInt encode; // work
|
||
|
|
|
||
|
|
// Get one coefficient as base-1000 and add the other
|
||
|
|
GETCOEFFTHOU(dfl, tac); // least-significant goes to [0]
|
||
|
|
ADDCOEFFTHOU(dfr, tac);
|
||
|
|
// here the sum of the MSDs (plus any carry) will be <10 due to
|
||
|
|
// the fastpath test earlier
|
||
|
|
|
||
|
|
// construct the result; low word is the same for both formats
|
||
|
|
encode =BIN2DPD[tac[0]];
|
||
|
|
encode|=BIN2DPD[tac[1]]<<10;
|
||
|
|
encode|=BIN2DPD[tac[2]]<<20;
|
||
|
|
encode|=BIN2DPD[tac[3]]<<30;
|
||
|
|
DFWORD(result, (DECBYTES/4)-1)=encode;
|
||
|
|
|
||
|
|
// collect next two declets (all that remains, for Double)
|
||
|
|
encode =BIN2DPD[tac[3]]>>2;
|
||
|
|
encode|=BIN2DPD[tac[4]]<<8;
|
||
|
|
|
||
|
|
#if QUAD
|
||
|
|
// complete and lay out middling words
|
||
|
|
encode|=BIN2DPD[tac[5]]<<18;
|
||
|
|
encode|=BIN2DPD[tac[6]]<<28;
|
||
|
|
DFWORD(result, 2)=encode;
|
||
|
|
|
||
|
|
encode =BIN2DPD[tac[6]]>>4;
|
||
|
|
encode|=BIN2DPD[tac[7]]<<6;
|
||
|
|
encode|=BIN2DPD[tac[8]]<<16;
|
||
|
|
encode|=BIN2DPD[tac[9]]<<26;
|
||
|
|
DFWORD(result, 1)=encode;
|
||
|
|
|
||
|
|
// and final two declets
|
||
|
|
encode =BIN2DPD[tac[9]]>>6;
|
||
|
|
encode|=BIN2DPD[tac[10]]<<4;
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// add exponent continuation and sign (from either argument)
|
||
|
|
encode|=sourhil & (ECONMASK | DECFLOAT_Sign);
|
||
|
|
|
||
|
|
// create lookup index = MSD + top two bits of biased exponent <<4
|
||
|
|
tac[DECLETS]|=(bexpl>>DECECONL)<<4;
|
||
|
|
encode|=DECCOMBFROM[tac[DECLETS]]; // add constructed combination field
|
||
|
|
DFWORD(result, 0)=encode; // complete
|
||
|
|
|
||
|
|
// decFloatShow(result, ">");
|
||
|
|
return result;
|
||
|
|
} // fast path OK
|
||
|
|
// drop through to slow path
|
||
|
|
} // low sum or Special(s)
|
||
|
|
|
||
|
|
/* Slow path required -- arguments are finite and might overflow, */
|
||
|
|
/* or require alignment, or might have different signs */
|
||
|
|
|
||
|
|
// now swap either exponents or argument pointers
|
||
|
|
if (bexpl<=bexpr) {
|
||
|
|
// original left is bigger
|
||
|
|
Int bexpswap=bexpl;
|
||
|
|
bexpl=bexpr;
|
||
|
|
bexpr=bexpswap;
|
||
|
|
// printf("left bigger\n");
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
const decFloat *dfswap=dfl;
|
||
|
|
dfl=dfr;
|
||
|
|
dfr=dfswap;
|
||
|
|
// printf("right bigger\n");
|
||
|
|
}
|
||
|
|
// [here dfl and bexpl refer to the datum with the larger exponent,
|
||
|
|
// of if the exponents are equal then the original LHS argument]
|
||
|
|
|
||
|
|
// if lhs is zero then result will be the rhs (now known to have
|
||
|
|
// the smaller exponent), which also may need to be tested for zero
|
||
|
|
// for the weird IEEE 754 sign rules
|
||
|
|
if (DFISZERO(dfl)) {
|
||
|
|
decCanonical(result, dfr); // clean copy
|
||
|
|
// "When the sum of two operands with opposite signs is
|
||
|
|
// exactly zero, the sign of that sum shall be '+' in all
|
||
|
|
// rounding modes except round toward -Infinity, in which
|
||
|
|
// mode that sign shall be '-'."
|
||
|
|
if (diffsign && DFISZERO(result)) {
|
||
|
|
DFWORD(result, 0)&=~DECFLOAT_Sign; // assume sign 0
|
||
|
|
if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
|
||
|
|
}
|
||
|
|
return result;
|
||
|
|
} // numfl is zero
|
||
|
|
// [here, LHS is non-zero; code below assumes that]
|
||
|
|
|
||
|
|
// Coefficients layout during the calculations to follow:
|
||
|
|
//
|
||
|
|
// Overlap case:
|
||
|
|
// +------------------------------------------------+
|
||
|
|
// acc: |0000| coeffa | tail B | |
|
||
|
|
// +------------------------------------------------+
|
||
|
|
// buf: |0000| pad0s | coeffb | |
|
||
|
|
// +------------------------------------------------+
|
||
|
|
//
|
||
|
|
// Touching coefficients or gap:
|
||
|
|
// +------------------------------------------------+
|
||
|
|
// acc: |0000| coeffa | gap | coeffb |
|
||
|
|
// +------------------------------------------------+
|
||
|
|
// [buf not used or needed; gap clamped to Pmax]
|
||
|
|
|
||
|
|
// lay out lhs coefficient into accumulator; this starts at acc+4
|
||
|
|
// for decDouble or acc+6 for decQuad so the LSD is word-
|
||
|
|
// aligned; the top word gap is there only in case a carry digit
|
||
|
|
// is prefixed after the add -- it does not need to be zeroed
|
||
|
|
#if DOUBLE
|
||
|
|
#define COFF 4 // offset into acc
|
||
|
|
#elif QUAD
|
||
|
|
UBFROMUS(acc+4, 0); // prefix 00
|
||
|
|
#define COFF 6 // offset into acc
|
||
|
|
#endif
|
||
|
|
|
||
|
|
GETCOEFF(dfl, acc+COFF); // decode from decFloat
|
||
|
|
ulsd=acc+COFF+DECPMAX-1;
|
||
|
|
umsd=acc+4; // [having this here avoids
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
{bcdnum tum;
|
||
|
|
tum.msd=umsd;
|
||
|
|
tum.lsd=ulsd;
|
||
|
|
tum.exponent=bexpl-DECBIAS;
|
||
|
|
tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
|
||
|
|
decShowNum(&tum, "dflx");}
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// if signs differ, take ten's complement of lhs (here the
|
||
|
|
// coefficient is subtracted from all-nines; the 1 is added during
|
||
|
|
// the later add cycle -- zeros to the right do not matter because
|
||
|
|
// the complement of zero is zero); these are fixed-length inverts
|
||
|
|
// where the lsd is known to be at a 4-byte boundary (so no borrow
|
||
|
|
// possible)
|
||
|
|
carry=0; // assume no carry
|
||
|
|
if (diffsign) {
|
||
|
|
carry=CARRYPAT; // for +1 during add
|
||
|
|
UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4));
|
||
|
|
UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8));
|
||
|
|
UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12));
|
||
|
|
UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16));
|
||
|
|
#if QUAD
|
||
|
|
UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20));
|
||
|
|
UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24));
|
||
|
|
UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28));
|
||
|
|
UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32));
|
||
|
|
UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36));
|
||
|
|
#endif
|
||
|
|
} // diffsign
|
||
|
|
|
||
|
|
// now process the rhs coefficient; if it cannot overlap lhs then
|
||
|
|
// it can be put straight into acc (with an appropriate gap, if
|
||
|
|
// needed) because no actual addition will be needed (except
|
||
|
|
// possibly to complete ten's complement)
|
||
|
|
overlap=DECPMAX-(bexpl-bexpr);
|
||
|
|
#if DECTRACE
|
||
|
|
printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS));
|
||
|
|
printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
|
||
|
|
#endif
|
||
|
|
|
||
|
|
if (overlap<=0) { // no overlap possible
|
||
|
|
uInt gap; // local work
|
||
|
|
// since a full addition is not needed, a ten's complement
|
||
|
|
// calculation started above may need to be completed
|
||
|
|
if (carry) {
|
||
|
|
for (ub=ulsd; *ub==9; ub--) *ub=0;
|
||
|
|
*ub+=1;
|
||
|
|
carry=0; // taken care of
|
||
|
|
}
|
||
|
|
// up to DECPMAX-1 digits of the final result can extend down
|
||
|
|
// below the LSD of the lhs, so if the gap is >DECPMAX then the
|
||
|
|
// rhs will be simply sticky bits. In this case the gap is
|
||
|
|
// clamped to DECPMAX and the exponent adjusted to suit [this is
|
||
|
|
// safe because the lhs is non-zero].
|
||
|
|
gap=-overlap;
|
||
|
|
if (gap>DECPMAX) {
|
||
|
|
bexpr+=gap-1;
|
||
|
|
gap=DECPMAX;
|
||
|
|
}
|
||
|
|
ub=ulsd+gap+1; // where MSD will go
|
||
|
|
// Fill the gap with 0s; note that there is no addition to do
|
||
|
|
ut=acc+COFF+DECPMAX; // start of gap
|
||
|
|
for (; ut<ub; ut+=4) UBFROMUI(ut, 0); // mind the gap
|
||
|
|
if (overlap<-DECPMAX) { // gap was > DECPMAX
|
||
|
|
*ub=(uByte)(!DFISZERO(dfr)); // make sticky digit
|
||
|
|
}
|
||
|
|
else { // need full coefficient
|
||
|
|
GETCOEFF(dfr, ub); // decode from decFloat
|
||
|
|
ub+=DECPMAX-1; // new LSD...
|
||
|
|
}
|
||
|
|
ulsd=ub; // save new LSD
|
||
|
|
} // no overlap possible
|
||
|
|
|
||
|
|
else { // overlap>0
|
||
|
|
// coefficients overlap (perhaps completely, although also
|
||
|
|
// perhaps only where zeros)
|
||
|
|
if (overlap==DECPMAX) { // aligned
|
||
|
|
ub=buf+COFF; // where msd will go
|
||
|
|
#if QUAD
|
||
|
|
UBFROMUS(buf+4, 0); // clear quad's 00
|
||
|
|
#endif
|
||
|
|
GETCOEFF(dfr, ub); // decode from decFloat
|
||
|
|
}
|
||
|
|
else { // unaligned
|
||
|
|
ub=buf+COFF+DECPMAX-overlap; // where MSD will go
|
||
|
|
// Fill the prefix gap with 0s; 8 will cover most common
|
||
|
|
// unalignments, so start with direct assignments (a loop is
|
||
|
|
// then used for any remaining -- the loop (and the one in a
|
||
|
|
// moment) is not then on the critical path because the number
|
||
|
|
// of additions is reduced by (at least) two in this case)
|
||
|
|
UBFROMUI(buf+4, 0); // [clears decQuad 00 too]
|
||
|
|
UBFROMUI(buf+8, 0);
|
||
|
|
if (ub>buf+12) {
|
||
|
|
ut=buf+12; // start any remaining
|
||
|
|
for (; ut<ub; ut+=4) UBFROMUI(ut, 0); // fill them
|
||
|
|
}
|
||
|
|
GETCOEFF(dfr, ub); // decode from decFloat
|
||
|
|
|
||
|
|
// now move tail of rhs across to main acc; again use direct
|
||
|
|
// copies for 8 digits-worth
|
||
|
|
UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX));
|
||
|
|
UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4));
|
||
|
|
if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
|
||
|
|
us=buf+COFF+DECPMAX+8; // source
|
||
|
|
ut=acc+COFF+DECPMAX+8; // target
|
||
|
|
for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us));
|
||
|
|
}
|
||
|
|
} // unaligned
|
||
|
|
|
||
|
|
ulsd=acc+(ub-buf+DECPMAX-1); // update LSD pointer
|
||
|
|
|
||
|
|
// Now do the add of the non-tail; this is all nicely aligned,
|
||
|
|
// and is over a multiple of four digits (because for Quad two
|
||
|
|
// zero digits were added on the left); words in both acc and
|
||
|
|
// buf (buf especially) will often be zero
|
||
|
|
// [byte-by-byte add, here, is about 15% slower total effect than
|
||
|
|
// the by-fours]
|
||
|
|
|
||
|
|
// Now effect the add; this is harder on a little-endian
|
||
|
|
// machine as the inter-digit carry cannot use the usual BCD
|
||
|
|
// addition trick because the bytes are loaded in the wrong order
|
||
|
|
// [this loop could be unrolled, but probably scarcely worth it]
|
||
|
|
|
||
|
|
ut=acc+COFF+DECPMAX-4; // target LSW (acc)
|
||
|
|
us=buf+COFF+DECPMAX-4; // source LSW (buf, to add to acc)
|
||
|
|
|
||
|
|
#if !DECLITEND
|
||
|
|
for (; ut>=acc+4; ut-=4, us-=4) { // big-endian add loop
|
||
|
|
// bcd8 add
|
||
|
|
carry+=UBTOUI(us); // rhs + carry
|
||
|
|
if (carry==0) continue; // no-op
|
||
|
|
carry+=UBTOUI(ut); // lhs
|
||
|
|
// Big-endian BCD adjust (uses internal carry)
|
||
|
|
carry+=0x76f6f6f6; // note top nibble not all bits
|
||
|
|
// apply BCD adjust and save
|
||
|
|
UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4));
|
||
|
|
carry>>=31; // true carry was at far left
|
||
|
|
} // add loop
|
||
|
|
#else
|
||
|
|
for (; ut>=acc+4; ut-=4, us-=4) { // little-endian add loop
|
||
|
|
// bcd8 add
|
||
|
|
carry+=UBTOUI(us); // rhs + carry
|
||
|
|
if (carry==0) continue; // no-op [common if unaligned]
|
||
|
|
carry+=UBTOUI(ut); // lhs
|
||
|
|
// Little-endian BCD adjust; inter-digit carry must be manual
|
||
|
|
// because the lsb from the array will be in the most-significant
|
||
|
|
// byte of carry
|
||
|
|
carry+=0x76767676; // note no inter-byte carries
|
||
|
|
carry+=(carry & 0x80000000)>>15;
|
||
|
|
carry+=(carry & 0x00800000)>>15;
|
||
|
|
carry+=(carry & 0x00008000)>>15;
|
||
|
|
carry-=(carry & 0x60606060)>>4; // BCD adjust back
|
||
|
|
UBFROMUI(ut, carry & 0x0f0f0f0f); // clear debris and save
|
||
|
|
// here, final carry-out bit is at 0x00000080; move it ready
|
||
|
|
// for next word-add (i.e., to 0x01000000)
|
||
|
|
carry=(carry & 0x00000080)<<17;
|
||
|
|
} // add loop
|
||
|
|
#endif
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
{bcdnum tum;
|
||
|
|
printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
|
||
|
|
tum.msd=umsd; // acc+4;
|
||
|
|
tum.lsd=ulsd;
|
||
|
|
tum.exponent=0;
|
||
|
|
tum.sign=0;
|
||
|
|
decShowNum(&tum, "dfadd");}
|
||
|
|
#endif
|
||
|
|
} // overlap possible
|
||
|
|
|
||
|
|
// ordering here is a little strange in order to have slowest path
|
||
|
|
// first in GCC asm listing
|
||
|
|
if (diffsign) { // subtraction
|
||
|
|
if (!carry) { // no carry out means RHS<LHS
|
||
|
|
// borrowed -- take ten's complement
|
||
|
|
// sign is lhs sign
|
||
|
|
num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
|
||
|
|
|
||
|
|
// invert the coefficient first by fours, then add one; space
|
||
|
|
// at the end of the buffer ensures the by-fours is always
|
||
|
|
// safe, but lsd+1 must be cleared to prevent a borrow
|
||
|
|
// if big-endian
|
||
|
|
#if !DECLITEND
|
||
|
|
*(ulsd+1)=0;
|
||
|
|
#endif
|
||
|
|
// there are always at least four coefficient words
|
||
|
|
UBFROMUI(umsd, 0x09090909-UBTOUI(umsd));
|
||
|
|
UBFROMUI(umsd+4, 0x09090909-UBTOUI(umsd+4));
|
||
|
|
UBFROMUI(umsd+8, 0x09090909-UBTOUI(umsd+8));
|
||
|
|
UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12));
|
||
|
|
#if DOUBLE
|
||
|
|
#define BNEXT 16
|
||
|
|
#elif QUAD
|
||
|
|
UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16));
|
||
|
|
UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20));
|
||
|
|
UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24));
|
||
|
|
UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28));
|
||
|
|
UBFROMUI(umsd+32, 0x09090909-UBTOUI(umsd+32));
|
||
|
|
#define BNEXT 36
|
||
|
|
#endif
|
||
|
|
if (ulsd>=umsd+BNEXT) { // unaligned
|
||
|
|
// eight will handle most unaligments for Double; 16 for Quad
|
||
|
|
UBFROMUI(umsd+BNEXT, 0x09090909-UBTOUI(umsd+BNEXT));
|
||
|
|
UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4));
|
||
|
|
#if DOUBLE
|
||
|
|
#define BNEXTY (BNEXT+8)
|
||
|
|
#elif QUAD
|
||
|
|
UBFROMUI(umsd+BNEXT+8, 0x09090909-UBTOUI(umsd+BNEXT+8));
|
||
|
|
UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12));
|
||
|
|
#define BNEXTY (BNEXT+16)
|
||
|
|
#endif
|
||
|
|
if (ulsd>=umsd+BNEXTY) { // very unaligned
|
||
|
|
ut=umsd+BNEXTY; // -> continue
|
||
|
|
for (;;ut+=4) {
|
||
|
|
UBFROMUI(ut, 0x09090909-UBTOUI(ut)); // invert four digits
|
||
|
|
if (ut>=ulsd-3) break; // all done
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
// complete the ten's complement by adding 1
|
||
|
|
for (ub=ulsd; *ub==9; ub--) *ub=0;
|
||
|
|
*ub+=1;
|
||
|
|
} // borrowed
|
||
|
|
|
||
|
|
else { // carry out means RHS>=LHS
|
||
|
|
num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
|
||
|
|
// all done except for the special IEEE 754 exact-zero-result
|
||
|
|
// rule (see above); while testing for zero, strip leading
|
||
|
|
// zeros (which will save decFinalize doing it) (this is in
|
||
|
|
// diffsign path, so carry impossible and true umsd is
|
||
|
|
// acc+COFF)
|
||
|
|
|
||
|
|
// Check the initial coefficient area using the fast macro;
|
||
|
|
// this will often be all that needs to be done (as on the
|
||
|
|
// worst-case path when the subtraction was aligned and
|
||
|
|
// full-length)
|
||
|
|
if (ISCOEFFZERO(acc+COFF)) {
|
||
|
|
umsd=acc+COFF+DECPMAX-1; // so far, so zero
|
||
|
|
if (ulsd>umsd) { // more to check
|
||
|
|
umsd++; // to align after checked area
|
||
|
|
for (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4;
|
||
|
|
for (; *umsd==0 && umsd<ulsd;) umsd++;
|
||
|
|
}
|
||
|
|
if (*umsd==0) { // must be true zero (and diffsign)
|
||
|
|
num.sign=0; // assume +
|
||
|
|
if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
// [else was not zero, might still have leading zeros]
|
||
|
|
} // subtraction gave positive result
|
||
|
|
} // diffsign
|
||
|
|
|
||
|
|
else { // same-sign addition
|
||
|
|
num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
|
||
|
|
#if DOUBLE
|
||
|
|
if (carry) { // only possible with decDouble
|
||
|
|
*(acc+3)=1; // [Quad has leading 00]
|
||
|
|
umsd=acc+3;
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
} // same sign
|
||
|
|
|
||
|
|
num.msd=umsd; // set MSD ..
|
||
|
|
num.lsd=ulsd; // .. and LSD
|
||
|
|
num.exponent=bexpr-DECBIAS; // set exponent to smaller, unbiassed
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
decFloatShow(dfl, "dfl");
|
||
|
|
decFloatShow(dfr, "dfr");
|
||
|
|
decShowNum(&num, "postadd");
|
||
|
|
#endif
|
||
|
|
return decFinalize(result, &num, set); // round, check, and lay out
|
||
|
|
} // decFloatAdd
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatAnd -- logical digitwise AND of two decFloats */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of ANDing dfl and dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result, which will be canonical with sign=0 */
|
||
|
|
/* */
|
||
|
|
/* The operands must be positive, finite with exponent q=0, and */
|
||
|
|
/* comprise just zeros and ones; if not, Invalid operation results. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatAnd(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
|
||
|
|
|| !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
|
||
|
|
// the operands are positive finite integers (q=0) with just 0s and 1s
|
||
|
|
#if DOUBLE
|
||
|
|
DFWORD(result, 0)=ZEROWORD
|
||
|
|
|((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
|
||
|
|
DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
|
||
|
|
#elif QUAD
|
||
|
|
DFWORD(result, 0)=ZEROWORD
|
||
|
|
|((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
|
||
|
|
DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
|
||
|
|
DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
|
||
|
|
DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
|
||
|
|
#endif
|
||
|
|
return result;
|
||
|
|
} // decFloatAnd
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCanonical -- copy a decFloat, making canonical */
|
||
|
|
/* */
|
||
|
|
/* result gets the canonicalized df */
|
||
|
|
/* df is the decFloat to copy and make canonical */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This works on specials, too; no error or exception is possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
|
||
|
|
return decCanonical(result, df);
|
||
|
|
} // decFloatCanonical
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatClass -- return the class of a decFloat */
|
||
|
|
/* */
|
||
|
|
/* df is the decFloat to test */
|
||
|
|
/* returns the decClass that df falls into */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
enum decClass decFloatClass(const decFloat *df) {
|
||
|
|
Int exp; // exponent
|
||
|
|
if (DFISSPECIAL(df)) {
|
||
|
|
if (DFISQNAN(df)) return DEC_CLASS_QNAN;
|
||
|
|
if (DFISSNAN(df)) return DEC_CLASS_SNAN;
|
||
|
|
// must be an infinity
|
||
|
|
if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
|
||
|
|
return DEC_CLASS_POS_INF;
|
||
|
|
}
|
||
|
|
if (DFISZERO(df)) { // quite common
|
||
|
|
if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
|
||
|
|
return DEC_CLASS_POS_ZERO;
|
||
|
|
}
|
||
|
|
// is finite and non-zero; similar code to decFloatIsNormal, here
|
||
|
|
// [this could be speeded up slightly by in-lining decFloatDigits]
|
||
|
|
exp=GETEXPUN(df) // get unbiased exponent ..
|
||
|
|
+decFloatDigits(df)-1; // .. and make adjusted exponent
|
||
|
|
if (exp>=DECEMIN) { // is normal
|
||
|
|
if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
|
||
|
|
return DEC_CLASS_POS_NORMAL;
|
||
|
|
}
|
||
|
|
// is subnormal
|
||
|
|
if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
|
||
|
|
return DEC_CLASS_POS_SUBNORMAL;
|
||
|
|
} // decFloatClass
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatClassString -- return the class of a decFloat as a string */
|
||
|
|
/* */
|
||
|
|
/* df is the decFloat to test */
|
||
|
|
/* returns a constant string describing the class df falls into */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
const char *decFloatClassString(const decFloat *df) {
|
||
|
|
enum decClass eclass=decFloatClass(df);
|
||
|
|
if (eclass==DEC_CLASS_POS_NORMAL) return DEC_ClassString_PN;
|
||
|
|
if (eclass==DEC_CLASS_NEG_NORMAL) return DEC_ClassString_NN;
|
||
|
|
if (eclass==DEC_CLASS_POS_ZERO) return DEC_ClassString_PZ;
|
||
|
|
if (eclass==DEC_CLASS_NEG_ZERO) return DEC_ClassString_NZ;
|
||
|
|
if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
|
||
|
|
if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
|
||
|
|
if (eclass==DEC_CLASS_POS_INF) return DEC_ClassString_PI;
|
||
|
|
if (eclass==DEC_CLASS_NEG_INF) return DEC_ClassString_NI;
|
||
|
|
if (eclass==DEC_CLASS_QNAN) return DEC_ClassString_QN;
|
||
|
|
if (eclass==DEC_CLASS_SNAN) return DEC_ClassString_SN;
|
||
|
|
return DEC_ClassString_UN; // Unknown
|
||
|
|
} // decFloatClassString
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of comparing dfl and dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result, which may be -1, 0, 1, or NaN (Unordered) */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCompare(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int comp; // work
|
||
|
|
// NaNs are handled as usual
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
// numeric comparison needed
|
||
|
|
comp=decNumCompare(dfl, dfr, 0);
|
||
|
|
decFloatZero(result);
|
||
|
|
if (comp==0) return result;
|
||
|
|
DFBYTE(result, DECBYTES-1)=0x01; // LSD=1
|
||
|
|
if (comp<0) DFBYTE(result, 0)|=0x80; // set sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatCompare
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCompareSignal -- compare two decFloats; all NaNs signal */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of comparing dfl and dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result, which may be -1, 0, 1, or NaN (Unordered) */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCompareSignal(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int comp; // work
|
||
|
|
// NaNs are handled as usual, except that all NaNs signal
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) {
|
||
|
|
set->status|=DEC_Invalid_operation;
|
||
|
|
return decNaNs(result, dfl, dfr, set);
|
||
|
|
}
|
||
|
|
// numeric comparison needed
|
||
|
|
comp=decNumCompare(dfl, dfr, 0);
|
||
|
|
decFloatZero(result);
|
||
|
|
if (comp==0) return result;
|
||
|
|
DFBYTE(result, DECBYTES-1)=0x01; // LSD=1
|
||
|
|
if (comp<0) DFBYTE(result, 0)|=0x80; // set sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatCompareSignal
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCompareTotal -- compare two decFloats with total ordering */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of comparing dfl and dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* returns result, which may be -1, 0, or 1 */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCompareTotal(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr) {
|
||
|
|
Int comp; // work
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
#if QUAD
|
||
|
|
uShort uswork; // ..
|
||
|
|
#endif
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) {
|
||
|
|
Int nanl, nanr; // work
|
||
|
|
// morph NaNs to +/- 1 or 2, leave numbers as 0
|
||
|
|
nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; // quiet > signalling
|
||
|
|
if (DFISSIGNED(dfl)) nanl=-nanl;
|
||
|
|
nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
|
||
|
|
if (DFISSIGNED(dfr)) nanr=-nanr;
|
||
|
|
if (nanl>nanr) comp=+1;
|
||
|
|
else if (nanl<nanr) comp=-1;
|
||
|
|
else { // NaNs are the same type and sign .. must compare payload
|
||
|
|
// buffers need +2 for QUAD
|
||
|
|
uByte bufl[DECPMAX+4]; // for LHS coefficient + foot
|
||
|
|
uByte bufr[DECPMAX+4]; // for RHS coefficient + foot
|
||
|
|
uByte *ub, *uc; // work
|
||
|
|
Int sigl; // signum of LHS
|
||
|
|
sigl=(DFISSIGNED(dfl) ? -1 : +1);
|
||
|
|
|
||
|
|
// decode the coefficients
|
||
|
|
// (shift both right two if Quad to make a multiple of four)
|
||
|
|
#if QUAD
|
||
|
|
UBFROMUS(bufl, 0);
|
||
|
|
UBFROMUS(bufr, 0);
|
||
|
|
#endif
|
||
|
|
GETCOEFF(dfl, bufl+QUAD*2); // decode from decFloat
|
||
|
|
GETCOEFF(dfr, bufr+QUAD*2); // ..
|
||
|
|
// all multiples of four, here
|
||
|
|
comp=0; // assume equal
|
||
|
|
for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
|
||
|
|
uInt ui=UBTOUI(ub);
|
||
|
|
if (ui==UBTOUI(uc)) continue; // so far so same
|
||
|
|
// about to find a winner; go by bytes in case little-endian
|
||
|
|
for (;; ub++, uc++) {
|
||
|
|
if (*ub==*uc) continue;
|
||
|
|
if (*ub>*uc) comp=sigl; // difference found
|
||
|
|
else comp=-sigl; // ..
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
} // same NaN type and sign
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
// numeric comparison needed
|
||
|
|
comp=decNumCompare(dfl, dfr, 1); // total ordering
|
||
|
|
}
|
||
|
|
decFloatZero(result);
|
||
|
|
if (comp==0) return result;
|
||
|
|
DFBYTE(result, DECBYTES-1)=0x01; // LSD=1
|
||
|
|
if (comp<0) DFBYTE(result, 0)|=0x80; // set sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatCompareTotal
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCompareTotalMag -- compare magnitudes with total ordering */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of comparing abs(dfl) and abs(dfr) */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* returns result, which may be -1, 0, or 1 */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCompareTotalMag(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr) {
|
||
|
|
decFloat a, b; // for copy if needed
|
||
|
|
// copy and redirect signed operand(s)
|
||
|
|
if (DFISSIGNED(dfl)) {
|
||
|
|
decFloatCopyAbs(&a, dfl);
|
||
|
|
dfl=&a;
|
||
|
|
}
|
||
|
|
if (DFISSIGNED(dfr)) {
|
||
|
|
decFloatCopyAbs(&b, dfr);
|
||
|
|
dfr=&b;
|
||
|
|
}
|
||
|
|
return decFloatCompareTotal(result, dfl, dfr);
|
||
|
|
} // decFloatCompareTotalMag
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCopy -- copy a decFloat as-is */
|
||
|
|
/* */
|
||
|
|
/* result gets the copy of dfl */
|
||
|
|
/* dfl is the decFloat to copy */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is a bitwise operation; no errors or exceptions are possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
|
||
|
|
if (dfl!=result) *result=*dfl; // copy needed
|
||
|
|
return result;
|
||
|
|
} // decFloatCopy
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */
|
||
|
|
/* */
|
||
|
|
/* result gets the copy of dfl with sign bit 0 */
|
||
|
|
/* dfl is the decFloat to copy */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is a bitwise operation; no errors or exceptions are possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
|
||
|
|
if (dfl!=result) *result=*dfl; // copy needed
|
||
|
|
DFBYTE(result, 0)&=~0x80; // zero sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatCopyAbs
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
|
||
|
|
/* */
|
||
|
|
/* result gets the copy of dfl with sign bit inverted */
|
||
|
|
/* dfl is the decFloat to copy */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is a bitwise operation; no errors or exceptions are possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
|
||
|
|
if (dfl!=result) *result=*dfl; // copy needed
|
||
|
|
DFBYTE(result, 0)^=0x80; // invert sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatCopyNegate
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatCopySign -- copy a decFloat with the sign of another */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of copying dfl with the sign of dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is a bitwise operation; no errors or exceptions are possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatCopySign(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr) {
|
||
|
|
uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); // save sign bit
|
||
|
|
if (dfl!=result) *result=*dfl; // copy needed
|
||
|
|
DFBYTE(result, 0)&=~0x80; // clear sign ..
|
||
|
|
DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); // .. and set saved
|
||
|
|
return result;
|
||
|
|
} // decFloatCopySign
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatDigits -- return the number of digits in a decFloat */
|
||
|
|
/* */
|
||
|
|
/* df is the decFloat to investigate */
|
||
|
|
/* returns the number of significant digits in the decFloat; a */
|
||
|
|
/* zero coefficient returns 1 as does an infinity (a NaN returns */
|
||
|
|
/* the number of digits in the payload) */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
// private macro to extract a declet according to provided formula
|
||
|
|
// (form), and if it is non-zero then return the calculated digits
|
||
|
|
// depending on the declet number (n), where n=0 for the most
|
||
|
|
// significant declet; uses uInt dpd for work
|
||
|
|
#define dpdlenchk(n, form) dpd=(form)&0x3ff; \
|
||
|
|
if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3])
|
||
|
|
// next one is used when it is known that the declet must be
|
||
|
|
// non-zero, or is the final zero declet
|
||
|
|
#define dpdlendun(n, form) dpd=(form)&0x3ff; \
|
||
|
|
if (dpd==0) return 1; \
|
||
|
|
return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3])
|
||
|
|
|
||
|
|
uInt decFloatDigits(const decFloat *df) {
|
||
|
|
uInt dpd; // work
|
||
|
|
uInt sourhi=DFWORD(df, 0); // top word from source decFloat
|
||
|
|
#if QUAD
|
||
|
|
uInt sourmh, sourml;
|
||
|
|
#endif
|
||
|
|
uInt sourlo;
|
||
|
|
|
||
|
|
if (DFISINF(df)) return 1;
|
||
|
|
// A NaN effectively has an MSD of 0; otherwise if non-zero MSD
|
||
|
|
// then the coefficient is full-length
|
||
|
|
if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
|
||
|
|
|
||
|
|
#if DOUBLE
|
||
|
|
if (sourhi&0x0003ffff) { // ends in first
|
||
|
|
dpdlenchk(0, sourhi>>8);
|
||
|
|
sourlo=DFWORD(df, 1);
|
||
|
|
dpdlendun(1, (sourhi<<2) | (sourlo>>30));
|
||
|
|
} // [cannot drop through]
|
||
|
|
sourlo=DFWORD(df, 1); // sourhi not involved now
|
||
|
|
if (sourlo&0xfff00000) { // in one of first two
|
||
|
|
dpdlenchk(1, sourlo>>30); // very rare
|
||
|
|
dpdlendun(2, sourlo>>20);
|
||
|
|
} // [cannot drop through]
|
||
|
|
dpdlenchk(3, sourlo>>10);
|
||
|
|
dpdlendun(4, sourlo);
|
||
|
|
// [cannot drop through]
|
||
|
|
|
||
|
|
#elif QUAD
|
||
|
|
if (sourhi&0x00003fff) { // ends in first
|
||
|
|
dpdlenchk(0, sourhi>>4);
|
||
|
|
sourmh=DFWORD(df, 1);
|
||
|
|
dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
|
||
|
|
} // [cannot drop through]
|
||
|
|
sourmh=DFWORD(df, 1);
|
||
|
|
if (sourmh) {
|
||
|
|
dpdlenchk(1, sourmh>>26);
|
||
|
|
dpdlenchk(2, sourmh>>16);
|
||
|
|
dpdlenchk(3, sourmh>>6);
|
||
|
|
sourml=DFWORD(df, 2);
|
||
|
|
dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
|
||
|
|
} // [cannot drop through]
|
||
|
|
sourml=DFWORD(df, 2);
|
||
|
|
if (sourml) {
|
||
|
|
dpdlenchk(4, sourml>>28);
|
||
|
|
dpdlenchk(5, sourml>>18);
|
||
|
|
dpdlenchk(6, sourml>>8);
|
||
|
|
sourlo=DFWORD(df, 3);
|
||
|
|
dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
|
||
|
|
} // [cannot drop through]
|
||
|
|
sourlo=DFWORD(df, 3);
|
||
|
|
if (sourlo&0xfff00000) { // in one of first two
|
||
|
|
dpdlenchk(7, sourlo>>30); // very rare
|
||
|
|
dpdlendun(8, sourlo>>20);
|
||
|
|
} // [cannot drop through]
|
||
|
|
dpdlenchk(9, sourlo>>10);
|
||
|
|
dpdlendun(10, sourlo);
|
||
|
|
// [cannot drop through]
|
||
|
|
#endif
|
||
|
|
} // decFloatDigits
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatDivide -- divide a decFloat by another */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of dividing dfl by dfr: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
// This is just a wrapper.
|
||
|
|
decFloat * decFloatDivide(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
return decDivide(result, dfl, dfr, set, DIVIDE);
|
||
|
|
} // decFloatDivide
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatDivideInteger -- integer divide a decFloat by another */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of dividing dfl by dfr: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatDivideInteger(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
return decDivide(result, dfl, dfr, set, DIVIDEINT);
|
||
|
|
} // decFloatDivideInteger
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatFMA -- multiply and add three decFloats, fused */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of (dfl*dfr)+dff with a single rounding */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* dff is the final decFloat (fhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
|
||
|
|
const decFloat *dfr, const decFloat *dff,
|
||
|
|
decContext *set) {
|
||
|
|
|
||
|
|
// The accumulator has the bytes needed for FiniteMultiply, plus
|
||
|
|
// one byte to the left in case of carry, plus DECPMAX+2 to the
|
||
|
|
// right for the final addition (up to full fhs + round & sticky)
|
||
|
|
#define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2))
|
||
|
|
uByte acc[FMALEN]; // for multiplied coefficient in BCD
|
||
|
|
// .. and for final result
|
||
|
|
bcdnum mul; // for multiplication result
|
||
|
|
bcdnum fin; // for final operand, expanded
|
||
|
|
uByte coe[ROUNDUP4(DECPMAX)]; // dff coefficient in BCD
|
||
|
|
bcdnum *hi, *lo; // bcdnum with higher/lower exponent
|
||
|
|
uInt diffsign; // non-zero if signs differ
|
||
|
|
uInt hipad; // pad digit for hi if needed
|
||
|
|
Int padding; // excess exponent
|
||
|
|
uInt carry; // +1 for ten's complement and during add
|
||
|
|
uByte *ub, *uh, *ul; // work
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
|
||
|
|
// handle all the special values [any special operand leads to a
|
||
|
|
// special result]
|
||
|
|
if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
|
||
|
|
decFloat proxy; // multiplication result proxy
|
||
|
|
// NaNs are handled as usual, giving priority to sNaNs
|
||
|
|
if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
|
||
|
|
// One or more of the three is infinite
|
||
|
|
// infinity times zero is bad
|
||
|
|
decFloatZero(&proxy);
|
||
|
|
if (DFISINF(dfl)) {
|
||
|
|
if (DFISZERO(dfr)) return decInvalid(result, set);
|
||
|
|
decInfinity(&proxy, &proxy);
|
||
|
|
}
|
||
|
|
else if (DFISINF(dfr)) {
|
||
|
|
if (DFISZERO(dfl)) return decInvalid(result, set);
|
||
|
|
decInfinity(&proxy, &proxy);
|
||
|
|
}
|
||
|
|
// compute sign of multiplication and place in proxy
|
||
|
|
DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
|
||
|
|
if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
|
||
|
|
// dff is Infinite
|
||
|
|
if (!DFISINF(&proxy)) return decInfinity(result, dff);
|
||
|
|
// both sides of addition are infinite; different sign is bad
|
||
|
|
if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
|
||
|
|
return decInvalid(result, set);
|
||
|
|
return decFloatCopy(result, &proxy);
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Here when all operands are finite */
|
||
|
|
|
||
|
|
// First multiply dfl*dfr
|
||
|
|
decFiniteMultiply(&mul, acc+1, dfl, dfr);
|
||
|
|
// The multiply is complete, exact and unbounded, and described in
|
||
|
|
// mul with the coefficient held in acc[1...]
|
||
|
|
|
||
|
|
// now add in dff; the algorithm is essentially the same as
|
||
|
|
// decFloatAdd, but the code is different because the code there
|
||
|
|
// is highly optimized for adding two numbers of the same size
|
||
|
|
fin.exponent=GETEXPUN(dff); // get dff exponent and sign
|
||
|
|
fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
|
||
|
|
diffsign=mul.sign^fin.sign; // note if signs differ
|
||
|
|
fin.msd=coe;
|
||
|
|
fin.lsd=coe+DECPMAX-1;
|
||
|
|
GETCOEFF(dff, coe); // extract the coefficient
|
||
|
|
|
||
|
|
// now set hi and lo so that hi points to whichever of mul and fin
|
||
|
|
// has the higher exponent and lo points to the other [don't care,
|
||
|
|
// if the same]. One coefficient will be in acc, the other in coe.
|
||
|
|
if (mul.exponent>=fin.exponent) {
|
||
|
|
hi=&mul;
|
||
|
|
lo=&fin;
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
hi=&fin;
|
||
|
|
lo=&mul;
|
||
|
|
}
|
||
|
|
|
||
|
|
// remove leading zeros on both operands; this will save time later
|
||
|
|
// and make testing for zero trivial (tests are safe because acc
|
||
|
|
// and coe are rounded up to uInts)
|
||
|
|
for (; UBTOUI(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
|
||
|
|
for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
|
||
|
|
for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
|
||
|
|
for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
|
||
|
|
|
||
|
|
// if hi is zero then result will be lo (which has the smaller
|
||
|
|
// exponent), which also may need to be tested for zero for the
|
||
|
|
// weird IEEE 754 sign rules
|
||
|
|
if (*hi->msd==0) { // hi is zero
|
||
|
|
// "When the sum of two operands with opposite signs is
|
||
|
|
// exactly zero, the sign of that sum shall be '+' in all
|
||
|
|
// rounding modes except round toward -Infinity, in which
|
||
|
|
// mode that sign shall be '-'."
|
||
|
|
if (diffsign) {
|
||
|
|
if (*lo->msd==0) { // lo is zero
|
||
|
|
lo->sign=0;
|
||
|
|
if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
|
||
|
|
} // diffsign && lo=0
|
||
|
|
} // diffsign
|
||
|
|
return decFinalize(result, lo, set); // may need clamping
|
||
|
|
} // numfl is zero
|
||
|
|
// [here, both are minimal length and hi is non-zero]
|
||
|
|
// (if lo is zero then padding with zeros may be needed, below)
|
||
|
|
|
||
|
|
// if signs differ, take the ten's complement of hi (zeros to the
|
||
|
|
// right do not matter because the complement of zero is zero); the
|
||
|
|
// +1 is done later, as part of the addition, inserted at the
|
||
|
|
// correct digit
|
||
|
|
hipad=0;
|
||
|
|
carry=0;
|
||
|
|
if (diffsign) {
|
||
|
|
hipad=9;
|
||
|
|
carry=1;
|
||
|
|
// exactly the correct number of digits must be inverted
|
||
|
|
for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(uh));
|
||
|
|
for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
|
||
|
|
}
|
||
|
|
|
||
|
|
// ready to add; note that hi has no leading zeros so gap
|
||
|
|
// calculation does not have to be as pessimistic as in decFloatAdd
|
||
|
|
// (this is much more like the arbitrary-precision algorithm in
|
||
|
|
// Rexx and decNumber)
|
||
|
|
|
||
|
|
// padding is the number of zeros that would need to be added to hi
|
||
|
|
// for its lsd to be aligned with the lsd of lo
|
||
|
|
padding=hi->exponent-lo->exponent;
|
||
|
|
// printf("FMA pad %ld\n", (LI)padding);
|
||
|
|
|
||
|
|
// the result of the addition will be built into the accumulator,
|
||
|
|
// starting from the far right; this could be either hi or lo, and
|
||
|
|
// will be aligned
|
||
|
|
ub=acc+FMALEN-1; // where lsd of result will go
|
||
|
|
ul=lo->lsd; // lsd of rhs
|
||
|
|
|
||
|
|
if (padding!=0) { // unaligned
|
||
|
|
// if the msd of lo is more than DECPMAX+2 digits to the right of
|
||
|
|
// the original msd of hi then it can be reduced to a single
|
||
|
|
// digit at the right place, as it stays clear of hi digits
|
||
|
|
// [it must be DECPMAX+2 because during a subtraction the msd
|
||
|
|
// could become 0 after a borrow from 1.000 to 0.9999...]
|
||
|
|
|
||
|
|
Int hilen=(Int)(hi->lsd-hi->msd+1); // length of hi
|
||
|
|
Int lolen=(Int)(lo->lsd-lo->msd+1); // and of lo
|
||
|
|
|
||
|
|
if (hilen+padding-lolen > DECPMAX+2) { // can reduce lo to single
|
||
|
|
// make sure it is virtually at least DECPMAX from hi->msd, at
|
||
|
|
// least to right of hi->lsd (in case of destructive subtract),
|
||
|
|
// and separated by at least two digits from either of those
|
||
|
|
// (the tricky DOUBLE case is when hi is a 1 that will become a
|
||
|
|
// 0.9999... by subtraction:
|
||
|
|
// hi: 1 E+16
|
||
|
|
// lo: .................1000000000000000 E-16
|
||
|
|
// which for the addition pads to:
|
||
|
|
// hi: 1000000000000000000 E-16
|
||
|
|
// lo: .................1000000000000000 E-16
|
||
|
|
Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
|
||
|
|
|
||
|
|
// printf("FMA reduce: %ld\n", (LI)reduce);
|
||
|
|
lo->lsd=lo->msd; // to single digit [maybe 0]
|
||
|
|
lo->exponent=newexp; // new lowest exponent
|
||
|
|
padding=hi->exponent-lo->exponent; // recalculate
|
||
|
|
ul=lo->lsd; // .. and repoint
|
||
|
|
}
|
||
|
|
|
||
|
|
// padding is still > 0, but will fit in acc (less leading carry slot)
|
||
|
|
#if DECCHECK
|
||
|
|
if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
|
||
|
|
if (hilen+padding+1>FMALEN)
|
||
|
|
printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (LI)padding);
|
||
|
|
// printf("FMA padding: %ld\n", (LI)padding);
|
||
|
|
#endif
|
||
|
|
|
||
|
|
// padding digits can now be set in the result; one or more of
|
||
|
|
// these will come from lo; others will be zeros in the gap
|
||
|
|
for (; ul-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) {
|
||
|
|
UBFROMUI(ub-3, UBTOUI(ul-3)); // [cannot overlap]
|
||
|
|
}
|
||
|
|
for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
|
||
|
|
for (;padding>0; padding--, ub--) *ub=0; // mind the gap
|
||
|
|
}
|
||
|
|
|
||
|
|
// addition now complete to the right of the rightmost digit of hi
|
||
|
|
uh=hi->lsd;
|
||
|
|
|
||
|
|
// dow do the add from hi->lsd to the left
|
||
|
|
// [bytewise, because either operand can run out at any time]
|
||
|
|
// carry was set up depending on ten's complement above
|
||
|
|
// first assume both operands have some digits
|
||
|
|
for (;; ub--) {
|
||
|
|
if (uh<hi->msd || ul<lo->msd) break;
|
||
|
|
*ub=(uByte)(carry+(*uh--)+(*ul--));
|
||
|
|
carry=0;
|
||
|
|
if (*ub<10) continue;
|
||
|
|
*ub-=10;
|
||
|
|
carry=1;
|
||
|
|
} // both loop
|
||
|
|
|
||
|
|
if (ul<lo->msd) { // to left of lo
|
||
|
|
for (;; ub--) {
|
||
|
|
if (uh<hi->msd) break;
|
||
|
|
*ub=(uByte)(carry+(*uh--)); // [+0]
|
||
|
|
carry=0;
|
||
|
|
if (*ub<10) continue;
|
||
|
|
*ub-=10;
|
||
|
|
carry=1;
|
||
|
|
} // hi loop
|
||
|
|
}
|
||
|
|
else { // to left of hi
|
||
|
|
for (;; ub--) {
|
||
|
|
if (ul<lo->msd) break;
|
||
|
|
*ub=(uByte)(carry+hipad+(*ul--));
|
||
|
|
carry=0;
|
||
|
|
if (*ub<10) continue;
|
||
|
|
*ub-=10;
|
||
|
|
carry=1;
|
||
|
|
} // lo loop
|
||
|
|
}
|
||
|
|
|
||
|
|
// addition complete -- now handle carry, borrow, etc.
|
||
|
|
// use lo to set up the num (its exponent is already correct, and
|
||
|
|
// sign usually is)
|
||
|
|
lo->msd=ub+1;
|
||
|
|
lo->lsd=acc+FMALEN-1;
|
||
|
|
// decShowNum(lo, "lo");
|
||
|
|
if (!diffsign) { // same-sign addition
|
||
|
|
if (carry) { // carry out
|
||
|
|
*ub=1; // place the 1 ..
|
||
|
|
lo->msd--; // .. and update
|
||
|
|
}
|
||
|
|
} // same sign
|
||
|
|
else { // signs differed (subtraction)
|
||
|
|
if (!carry) { // no carry out means hi<lo
|
||
|
|
// borrowed -- take ten's complement of the right digits
|
||
|
|
lo->sign=hi->sign; // sign is lhs sign
|
||
|
|
for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(ul));
|
||
|
|
for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); // [leaves ul at lsd+1]
|
||
|
|
// complete the ten's complement by adding 1 [cannot overrun]
|
||
|
|
for (ul--; *ul==9; ul--) *ul=0;
|
||
|
|
*ul+=1;
|
||
|
|
} // borrowed
|
||
|
|
else { // carry out means hi>=lo
|
||
|
|
// sign to use is lo->sign
|
||
|
|
// all done except for the special IEEE 754 exact-zero-result
|
||
|
|
// rule (see above); while testing for zero, strip leading
|
||
|
|
// zeros (which will save decFinalize doing it)
|
||
|
|
for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
|
||
|
|
for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
|
||
|
|
if (*lo->msd==0) { // must be true zero (and diffsign)
|
||
|
|
lo->sign=0; // assume +
|
||
|
|
if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
|
||
|
|
}
|
||
|
|
// [else was not zero, might still have leading zeros]
|
||
|
|
} // subtraction gave positive result
|
||
|
|
} // diffsign
|
||
|
|
|
||
|
|
#if DECCHECK
|
||
|
|
// assert no left underrun
|
||
|
|
if (lo->msd<acc) {
|
||
|
|
printf("FMA underrun by %ld \n", (LI)(acc-lo->msd));
|
||
|
|
}
|
||
|
|
#endif
|
||
|
|
|
||
|
|
return decFinalize(result, lo, set); // round, check, and lay out
|
||
|
|
} // decFloatFMA
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatFromInt -- initialise a decFloat from an Int */
|
||
|
|
/* */
|
||
|
|
/* result gets the converted Int */
|
||
|
|
/* n is the Int to convert */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* The result is Exact; no errors or exceptions are possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatFromInt32(decFloat *result, Int n) {
|
||
|
|
uInt u=(uInt)n; // copy as bits
|
||
|
|
uInt encode; // work
|
||
|
|
DFWORD(result, 0)=ZEROWORD; // always
|
||
|
|
#if QUAD
|
||
|
|
DFWORD(result, 1)=0;
|
||
|
|
DFWORD(result, 2)=0;
|
||
|
|
#endif
|
||
|
|
if (n<0) { // handle -n with care
|
||
|
|
// [This can be done without the test, but is then slightly slower]
|
||
|
|
u=(~u)+1;
|
||
|
|
DFWORD(result, 0)|=DECFLOAT_Sign;
|
||
|
|
}
|
||
|
|
// Since the maximum value of u now is 2**31, only the low word of
|
||
|
|
// result is affected
|
||
|
|
encode=BIN2DPD[u%1000];
|
||
|
|
u/=1000;
|
||
|
|
encode|=BIN2DPD[u%1000]<<10;
|
||
|
|
u/=1000;
|
||
|
|
encode|=BIN2DPD[u%1000]<<20;
|
||
|
|
u/=1000; // now 0, 1, or 2
|
||
|
|
encode|=u<<30;
|
||
|
|
DFWORD(result, DECWORDS-1)=encode;
|
||
|
|
return result;
|
||
|
|
} // decFloatFromInt32
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatFromUInt -- initialise a decFloat from a uInt */
|
||
|
|
/* */
|
||
|
|
/* result gets the converted uInt */
|
||
|
|
/* n is the uInt to convert */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* The result is Exact; no errors or exceptions are possible. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
|
||
|
|
uInt encode; // work
|
||
|
|
DFWORD(result, 0)=ZEROWORD; // always
|
||
|
|
#if QUAD
|
||
|
|
DFWORD(result, 1)=0;
|
||
|
|
DFWORD(result, 2)=0;
|
||
|
|
#endif
|
||
|
|
encode=BIN2DPD[u%1000];
|
||
|
|
u/=1000;
|
||
|
|
encode|=BIN2DPD[u%1000]<<10;
|
||
|
|
u/=1000;
|
||
|
|
encode|=BIN2DPD[u%1000]<<20;
|
||
|
|
u/=1000; // now 0 -> 4
|
||
|
|
encode|=u<<30;
|
||
|
|
DFWORD(result, DECWORDS-1)=encode;
|
||
|
|
DFWORD(result, DECWORDS-2)|=u>>2; // rarely non-zero
|
||
|
|
return result;
|
||
|
|
} // decFloatFromUInt32
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatInvert -- logical digitwise INVERT of a decFloat */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of INVERTing df */
|
||
|
|
/* df is the decFloat to invert */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result, which will be canonical with sign=0 */
|
||
|
|
/* */
|
||
|
|
/* The operand must be positive, finite with exponent q=0, and */
|
||
|
|
/* comprise just zeros and ones; if not, Invalid operation results. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatInvert(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set) {
|
||
|
|
uInt sourhi=DFWORD(df, 0); // top word of dfs
|
||
|
|
|
||
|
|
if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
|
||
|
|
// the operand is a finite integer (q=0)
|
||
|
|
#if DOUBLE
|
||
|
|
DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
|
||
|
|
DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491;
|
||
|
|
#elif QUAD
|
||
|
|
DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
|
||
|
|
DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449;
|
||
|
|
DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124;
|
||
|
|
DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491;
|
||
|
|
#endif
|
||
|
|
return result;
|
||
|
|
} // decFloatInvert
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatIs -- decFloat tests (IsSigned, etc.) */
|
||
|
|
/* */
|
||
|
|
/* df is the decFloat to test */
|
||
|
|
/* returns 0 or 1 in a uInt */
|
||
|
|
/* */
|
||
|
|
/* Many of these could be macros, but having them as real functions */
|
||
|
|
/* is a little cleaner (and they can be referred to here by the */
|
||
|
|
/* generic names) */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
uInt decFloatIsCanonical(const decFloat *df) {
|
||
|
|
if (DFISSPECIAL(df)) {
|
||
|
|
if (DFISINF(df)) {
|
||
|
|
if (DFWORD(df, 0)&ECONMASK) return 0; // exponent continuation
|
||
|
|
if (!DFISCCZERO(df)) return 0; // coefficient continuation
|
||
|
|
return 1;
|
||
|
|
}
|
||
|
|
// is a NaN
|
||
|
|
if (DFWORD(df, 0)&ECONNANMASK) return 0; // exponent continuation
|
||
|
|
if (DFISCCZERO(df)) return 1; // coefficient continuation
|
||
|
|
// drop through to check payload
|
||
|
|
}
|
||
|
|
{ // declare block
|
||
|
|
#if DOUBLE
|
||
|
|
uInt sourhi=DFWORD(df, 0);
|
||
|
|
uInt sourlo=DFWORD(df, 1);
|
||
|
|
if (CANONDPDOFF(sourhi, 8)
|
||
|
|
&& CANONDPDTWO(sourhi, sourlo, 30)
|
||
|
|
&& CANONDPDOFF(sourlo, 20)
|
||
|
|
&& CANONDPDOFF(sourlo, 10)
|
||
|
|
&& CANONDPDOFF(sourlo, 0)) return 1;
|
||
|
|
#elif QUAD
|
||
|
|
uInt sourhi=DFWORD(df, 0);
|
||
|
|
uInt sourmh=DFWORD(df, 1);
|
||
|
|
uInt sourml=DFWORD(df, 2);
|
||
|
|
uInt sourlo=DFWORD(df, 3);
|
||
|
|
if (CANONDPDOFF(sourhi, 4)
|
||
|
|
&& CANONDPDTWO(sourhi, sourmh, 26)
|
||
|
|
&& CANONDPDOFF(sourmh, 16)
|
||
|
|
&& CANONDPDOFF(sourmh, 6)
|
||
|
|
&& CANONDPDTWO(sourmh, sourml, 28)
|
||
|
|
&& CANONDPDOFF(sourml, 18)
|
||
|
|
&& CANONDPDOFF(sourml, 8)
|
||
|
|
&& CANONDPDTWO(sourml, sourlo, 30)
|
||
|
|
&& CANONDPDOFF(sourlo, 20)
|
||
|
|
&& CANONDPDOFF(sourlo, 10)
|
||
|
|
&& CANONDPDOFF(sourlo, 0)) return 1;
|
||
|
|
#endif
|
||
|
|
} // block
|
||
|
|
return 0; // a declet is non-canonical
|
||
|
|
}
|
||
|
|
|
||
|
|
uInt decFloatIsFinite(const decFloat *df) {
|
||
|
|
return !DFISSPECIAL(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsInfinite(const decFloat *df) {
|
||
|
|
return DFISINF(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsInteger(const decFloat *df) {
|
||
|
|
return DFISINT(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsLogical(const decFloat *df) {
|
||
|
|
return DFISUINT01(df) & DFISCC01(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsNaN(const decFloat *df) {
|
||
|
|
return DFISNAN(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsNegative(const decFloat *df) {
|
||
|
|
return DFISSIGNED(df) && !DFISZERO(df) && !DFISNAN(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsNormal(const decFloat *df) {
|
||
|
|
Int exp; // exponent
|
||
|
|
if (DFISSPECIAL(df)) return 0;
|
||
|
|
if (DFISZERO(df)) return 0;
|
||
|
|
// is finite and non-zero
|
||
|
|
exp=GETEXPUN(df) // get unbiased exponent ..
|
||
|
|
+decFloatDigits(df)-1; // .. and make adjusted exponent
|
||
|
|
return (exp>=DECEMIN); // < DECEMIN is subnormal
|
||
|
|
}
|
||
|
|
uInt decFloatIsPositive(const decFloat *df) {
|
||
|
|
return !DFISSIGNED(df) && !DFISZERO(df) && !DFISNAN(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsSignaling(const decFloat *df) {
|
||
|
|
return DFISSNAN(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsSignalling(const decFloat *df) {
|
||
|
|
return DFISSNAN(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsSigned(const decFloat *df) {
|
||
|
|
return DFISSIGNED(df);
|
||
|
|
}
|
||
|
|
uInt decFloatIsSubnormal(const decFloat *df) {
|
||
|
|
if (DFISSPECIAL(df)) return 0;
|
||
|
|
// is finite
|
||
|
|
if (decFloatIsNormal(df)) return 0;
|
||
|
|
// it is <Nmin, but could be zero
|
||
|
|
if (DFISZERO(df)) return 0;
|
||
|
|
return 1; // is subnormal
|
||
|
|
}
|
||
|
|
uInt decFloatIsZero(const decFloat *df) {
|
||
|
|
return DFISZERO(df);
|
||
|
|
} // decFloatIs...
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatLogB -- return adjusted exponent, by 754 rules */
|
||
|
|
/* */
|
||
|
|
/* result gets the adjusted exponent as an integer, or a NaN etc. */
|
||
|
|
/* df is the decFloat to be examined */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* Notable cases: */
|
||
|
|
/* A<0 -> Use |A| */
|
||
|
|
/* A=0 -> -Infinity (Division by zero) */
|
||
|
|
/* A=Infinite -> +Infinity (Exact) */
|
||
|
|
/* A=1 exactly -> 0 (Exact) */
|
||
|
|
/* NaNs are propagated as usual */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatLogB(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set) {
|
||
|
|
Int ae; // adjusted exponent
|
||
|
|
if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
|
||
|
|
if (DFISINF(df)) {
|
||
|
|
DFWORD(result, 0)=0; // need +ve
|
||
|
|
return decInfinity(result, result); // canonical +Infinity
|
||
|
|
}
|
||
|
|
if (DFISZERO(df)) {
|
||
|
|
set->status|=DEC_Division_by_zero; // as per 754
|
||
|
|
DFWORD(result, 0)=DECFLOAT_Sign; // make negative
|
||
|
|
return decInfinity(result, result); // canonical -Infinity
|
||
|
|
}
|
||
|
|
ae=GETEXPUN(df) // get unbiased exponent ..
|
||
|
|
+decFloatDigits(df)-1; // .. and make adjusted exponent
|
||
|
|
// ae has limited range (3 digits for DOUBLE and 4 for QUAD), so
|
||
|
|
// it is worth using a special case of decFloatFromInt32
|
||
|
|
DFWORD(result, 0)=ZEROWORD; // always
|
||
|
|
if (ae<0) {
|
||
|
|
DFWORD(result, 0)|=DECFLOAT_Sign; // -0 so far
|
||
|
|
ae=-ae;
|
||
|
|
}
|
||
|
|
#if DOUBLE
|
||
|
|
DFWORD(result, 1)=BIN2DPD[ae]; // a single declet
|
||
|
|
#elif QUAD
|
||
|
|
DFWORD(result, 1)=0;
|
||
|
|
DFWORD(result, 2)=0;
|
||
|
|
DFWORD(result, 3)=(ae/1000)<<10; // is <10, so need no DPD encode
|
||
|
|
DFWORD(result, 3)|=BIN2DPD[ae%1000];
|
||
|
|
#endif
|
||
|
|
return result;
|
||
|
|
} // decFloatLogB
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatMax -- return maxnum of two operands */
|
||
|
|
/* */
|
||
|
|
/* result gets the chosen decFloat */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* If just one operand is a quiet NaN it is ignored. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatMax(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int comp;
|
||
|
|
if (DFISNAN(dfl)) {
|
||
|
|
// sNaN or both NaNs leads to normal NaN processing
|
||
|
|
if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
return decCanonical(result, dfr); // RHS is numeric
|
||
|
|
}
|
||
|
|
if (DFISNAN(dfr)) {
|
||
|
|
// sNaN leads to normal NaN processing (both NaN handled above)
|
||
|
|
if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
return decCanonical(result, dfl); // LHS is numeric
|
||
|
|
}
|
||
|
|
// Both operands are numeric; numeric comparison needed -- use
|
||
|
|
// total order for a well-defined choice (and +0 > -0)
|
||
|
|
comp=decNumCompare(dfl, dfr, 1);
|
||
|
|
if (comp>=0) return decCanonical(result, dfl);
|
||
|
|
return decCanonical(result, dfr);
|
||
|
|
} // decFloatMax
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatMaxMag -- return maxnummag of two operands */
|
||
|
|
/* */
|
||
|
|
/* result gets the chosen decFloat */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* Returns according to the magnitude comparisons if both numeric and */
|
||
|
|
/* unequal, otherwise returns maxnum */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatMaxMag(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int comp;
|
||
|
|
decFloat absl, absr;
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
|
||
|
|
|
||
|
|
decFloatCopyAbs(&absl, dfl);
|
||
|
|
decFloatCopyAbs(&absr, dfr);
|
||
|
|
comp=decNumCompare(&absl, &absr, 0);
|
||
|
|
if (comp>0) return decCanonical(result, dfl);
|
||
|
|
if (comp<0) return decCanonical(result, dfr);
|
||
|
|
return decFloatMax(result, dfl, dfr, set);
|
||
|
|
} // decFloatMaxMag
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatMin -- return minnum of two operands */
|
||
|
|
/* */
|
||
|
|
/* result gets the chosen decFloat */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* If just one operand is a quiet NaN it is ignored. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatMin(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int comp;
|
||
|
|
if (DFISNAN(dfl)) {
|
||
|
|
// sNaN or both NaNs leads to normal NaN processing
|
||
|
|
if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
return decCanonical(result, dfr); // RHS is numeric
|
||
|
|
}
|
||
|
|
if (DFISNAN(dfr)) {
|
||
|
|
// sNaN leads to normal NaN processing (both NaN handled above)
|
||
|
|
if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
return decCanonical(result, dfl); // LHS is numeric
|
||
|
|
}
|
||
|
|
// Both operands are numeric; numeric comparison needed -- use
|
||
|
|
// total order for a well-defined choice (and +0 > -0)
|
||
|
|
comp=decNumCompare(dfl, dfr, 1);
|
||
|
|
if (comp<=0) return decCanonical(result, dfl);
|
||
|
|
return decCanonical(result, dfr);
|
||
|
|
} // decFloatMin
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatMinMag -- return minnummag of two operands */
|
||
|
|
/* */
|
||
|
|
/* result gets the chosen decFloat */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* Returns according to the magnitude comparisons if both numeric and */
|
||
|
|
/* unequal, otherwise returns minnum */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatMinMag(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int comp;
|
||
|
|
decFloat absl, absr;
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
|
||
|
|
|
||
|
|
decFloatCopyAbs(&absl, dfl);
|
||
|
|
decFloatCopyAbs(&absr, dfr);
|
||
|
|
comp=decNumCompare(&absl, &absr, 0);
|
||
|
|
if (comp<0) return decCanonical(result, dfl);
|
||
|
|
if (comp>0) return decCanonical(result, dfr);
|
||
|
|
return decFloatMin(result, dfl, dfr, set);
|
||
|
|
} // decFloatMinMag
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatMinus -- negate value, heeding NaNs, etc. */
|
||
|
|
/* */
|
||
|
|
/* result gets the canonicalized 0-df */
|
||
|
|
/* df is the decFloat to minus */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This has the same effect as 0-df where the exponent of the zero is */
|
||
|
|
/* the same as that of df (if df is finite). */
|
||
|
|
/* The effect is also the same as decFloatCopyNegate except that NaNs */
|
||
|
|
/* are handled normally (the sign of a NaN is not affected, and an */
|
||
|
|
/* sNaN will signal), the result is canonical, and zero gets sign 0. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatMinus(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set) {
|
||
|
|
if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
|
||
|
|
decCanonical(result, df); // copy and check
|
||
|
|
if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; // turn off sign bit
|
||
|
|
else DFBYTE(result, 0)^=0x80; // flip sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatMinus
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatMultiply -- multiply two decFloats */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of multiplying dfl and dfr: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatMultiply(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
bcdnum num; // for final conversion
|
||
|
|
uByte bcdacc[DECPMAX9*18+1]; // for coefficent in BCD
|
||
|
|
|
||
|
|
if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { // either is special?
|
||
|
|
// NaNs are handled as usual
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
// infinity times zero is bad
|
||
|
|
if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
|
||
|
|
if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
|
||
|
|
// both infinite; return canonical infinity with computed sign
|
||
|
|
DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); // compute sign
|
||
|
|
return decInfinity(result, result);
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Here when both operands are finite */
|
||
|
|
decFiniteMultiply(&num, bcdacc, dfl, dfr);
|
||
|
|
return decFinalize(result, &num, set); // round, check, and lay out
|
||
|
|
} // decFloatMultiply
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatNextMinus -- next towards -Infinity */
|
||
|
|
/* */
|
||
|
|
/* result gets the next lesser decFloat */
|
||
|
|
/* dfl is the decFloat to start with */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is 754 nextdown; Invalid is the only status possible (from */
|
||
|
|
/* an sNaN). */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
|
||
|
|
decContext *set) {
|
||
|
|
decFloat delta; // tiny increment
|
||
|
|
uInt savestat; // saves status
|
||
|
|
enum rounding saveround; // .. and mode
|
||
|
|
|
||
|
|
// +Infinity is the special case
|
||
|
|
if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
|
||
|
|
DFSETNMAX(result);
|
||
|
|
return result; // [no status to set]
|
||
|
|
}
|
||
|
|
// other cases are effected by sutracting a tiny delta -- this
|
||
|
|
// should be done in a wider format as the delta is unrepresentable
|
||
|
|
// here (but can be done with normal add if the sign of zero is
|
||
|
|
// treated carefully, because no Inexactitude is interesting);
|
||
|
|
// rounding to -Infinity then pushes the result to next below
|
||
|
|
decFloatZero(&delta); // set up tiny delta
|
||
|
|
DFWORD(&delta, DECWORDS-1)=1; // coefficient=1
|
||
|
|
DFWORD(&delta, 0)=DECFLOAT_Sign; // Sign=1 + biased exponent=0
|
||
|
|
// set up for the directional round
|
||
|
|
saveround=set->round; // save mode
|
||
|
|
set->round=DEC_ROUND_FLOOR; // .. round towards -Infinity
|
||
|
|
savestat=set->status; // save status
|
||
|
|
decFloatAdd(result, dfl, &delta, set);
|
||
|
|
// Add rules mess up the sign when going from +Ntiny to 0
|
||
|
|
if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; // correct
|
||
|
|
set->status&=DEC_Invalid_operation; // preserve only sNaN status
|
||
|
|
set->status|=savestat; // restore pending flags
|
||
|
|
set->round=saveround; // .. and mode
|
||
|
|
return result;
|
||
|
|
} // decFloatNextMinus
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatNextPlus -- next towards +Infinity */
|
||
|
|
/* */
|
||
|
|
/* result gets the next larger decFloat */
|
||
|
|
/* dfl is the decFloat to start with */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is 754 nextup; Invalid is the only status possible (from */
|
||
|
|
/* an sNaN). */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
|
||
|
|
decContext *set) {
|
||
|
|
uInt savestat; // saves status
|
||
|
|
enum rounding saveround; // .. and mode
|
||
|
|
decFloat delta; // tiny increment
|
||
|
|
|
||
|
|
// -Infinity is the special case
|
||
|
|
if (DFISINF(dfl) && DFISSIGNED(dfl)) {
|
||
|
|
DFSETNMAX(result);
|
||
|
|
DFWORD(result, 0)|=DECFLOAT_Sign; // make negative
|
||
|
|
return result; // [no status to set]
|
||
|
|
}
|
||
|
|
// other cases are effected by sutracting a tiny delta -- this
|
||
|
|
// should be done in a wider format as the delta is unrepresentable
|
||
|
|
// here (but can be done with normal add if the sign of zero is
|
||
|
|
// treated carefully, because no Inexactitude is interesting);
|
||
|
|
// rounding to +Infinity then pushes the result to next above
|
||
|
|
decFloatZero(&delta); // set up tiny delta
|
||
|
|
DFWORD(&delta, DECWORDS-1)=1; // coefficient=1
|
||
|
|
DFWORD(&delta, 0)=0; // Sign=0 + biased exponent=0
|
||
|
|
// set up for the directional round
|
||
|
|
saveround=set->round; // save mode
|
||
|
|
set->round=DEC_ROUND_CEILING; // .. round towards +Infinity
|
||
|
|
savestat=set->status; // save status
|
||
|
|
decFloatAdd(result, dfl, &delta, set);
|
||
|
|
// Add rules mess up the sign when going from -Ntiny to -0
|
||
|
|
if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; // correct
|
||
|
|
set->status&=DEC_Invalid_operation; // preserve only sNaN status
|
||
|
|
set->status|=savestat; // restore pending flags
|
||
|
|
set->round=saveround; // .. and mode
|
||
|
|
return result;
|
||
|
|
} // decFloatNextPlus
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatNextToward -- next towards a decFloat */
|
||
|
|
/* */
|
||
|
|
/* result gets the next decFloat */
|
||
|
|
/* dfl is the decFloat to start with */
|
||
|
|
/* dfr is the decFloat to move toward */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is 754-1985 nextafter, as modified during revision (dropped */
|
||
|
|
/* from 754-2008); status may be set unless the result is a normal */
|
||
|
|
/* number. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatNextToward(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
decFloat delta; // tiny increment or decrement
|
||
|
|
decFloat pointone; // 1e-1
|
||
|
|
uInt savestat; // saves status
|
||
|
|
enum rounding saveround; // .. and mode
|
||
|
|
uInt deltatop; // top word for delta
|
||
|
|
Int comp; // work
|
||
|
|
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
// Both are numeric, so Invalid no longer a possibility
|
||
|
|
comp=decNumCompare(dfl, dfr, 0);
|
||
|
|
if (comp==0) return decFloatCopySign(result, dfl, dfr); // equal
|
||
|
|
// unequal; do NextPlus or NextMinus but with different status rules
|
||
|
|
|
||
|
|
if (comp<0) { // lhs<rhs, do NextPlus, see above for commentary
|
||
|
|
if (DFISINF(dfl) && DFISSIGNED(dfl)) { // -Infinity special case
|
||
|
|
DFSETNMAX(result);
|
||
|
|
DFWORD(result, 0)|=DECFLOAT_Sign;
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
saveround=set->round; // save mode
|
||
|
|
set->round=DEC_ROUND_CEILING; // .. round towards +Infinity
|
||
|
|
deltatop=0; // positive delta
|
||
|
|
}
|
||
|
|
else { // lhs>rhs, do NextMinus, see above for commentary
|
||
|
|
if (DFISINF(dfl) && !DFISSIGNED(dfl)) { // +Infinity special case
|
||
|
|
DFSETNMAX(result);
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
saveround=set->round; // save mode
|
||
|
|
set->round=DEC_ROUND_FLOOR; // .. round towards -Infinity
|
||
|
|
deltatop=DECFLOAT_Sign; // negative delta
|
||
|
|
}
|
||
|
|
savestat=set->status; // save status
|
||
|
|
// Here, Inexact is needed where appropriate (and hence Underflow,
|
||
|
|
// etc.). Therefore the tiny delta which is otherwise
|
||
|
|
// unrepresentable (see NextPlus and NextMinus) is constructed
|
||
|
|
// using the multiplication of FMA.
|
||
|
|
decFloatZero(&delta); // set up tiny delta
|
||
|
|
DFWORD(&delta, DECWORDS-1)=1; // coefficient=1
|
||
|
|
DFWORD(&delta, 0)=deltatop; // Sign + biased exponent=0
|
||
|
|
decFloatFromString(&pointone, "1E-1", set); // set up multiplier
|
||
|
|
decFloatFMA(result, &delta, &pointone, dfl, set);
|
||
|
|
// [Delta is truly tiny, so no need to correct sign of zero]
|
||
|
|
// use new status unless the result is normal
|
||
|
|
if (decFloatIsNormal(result)) set->status=savestat; // else goes forward
|
||
|
|
set->round=saveround; // restore mode
|
||
|
|
return result;
|
||
|
|
} // decFloatNextToward
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatOr -- logical digitwise OR of two decFloats */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of ORing dfl and dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result, which will be canonical with sign=0 */
|
||
|
|
/* */
|
||
|
|
/* The operands must be positive, finite with exponent q=0, and */
|
||
|
|
/* comprise just zeros and ones; if not, Invalid operation results. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatOr(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
|
||
|
|
|| !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
|
||
|
|
// the operands are positive finite integers (q=0) with just 0s and 1s
|
||
|
|
#if DOUBLE
|
||
|
|
DFWORD(result, 0)=ZEROWORD
|
||
|
|
|((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
|
||
|
|
DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
|
||
|
|
#elif QUAD
|
||
|
|
DFWORD(result, 0)=ZEROWORD
|
||
|
|
|((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
|
||
|
|
DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
|
||
|
|
DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
|
||
|
|
DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
|
||
|
|
#endif
|
||
|
|
return result;
|
||
|
|
} // decFloatOr
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatPlus -- add value to 0, heeding NaNs, etc. */
|
||
|
|
/* */
|
||
|
|
/* result gets the canonicalized 0+df */
|
||
|
|
/* df is the decFloat to plus */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This has the same effect as 0+df where the exponent of the zero is */
|
||
|
|
/* the same as that of df (if df is finite). */
|
||
|
|
/* The effect is also the same as decFloatCopy except that NaNs */
|
||
|
|
/* are handled normally (the sign of a NaN is not affected, and an */
|
||
|
|
/* sNaN will signal), the result is canonical, and zero gets sign 0. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatPlus(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set) {
|
||
|
|
if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
|
||
|
|
decCanonical(result, df); // copy and check
|
||
|
|
if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; // turn off sign bit
|
||
|
|
return result;
|
||
|
|
} // decFloatPlus
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatQuantize -- quantize a decFloat */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of quantizing dfl to match dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs), which sets the exponent */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* Unless there is an error or the result is infinite, the exponent */
|
||
|
|
/* of result is guaranteed to be the same as that of dfr. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatQuantize(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int explb, exprb; // left and right biased exponents
|
||
|
|
uByte *ulsd; // local LSD pointer
|
||
|
|
uByte *ub, *uc; // work
|
||
|
|
Int drop; // ..
|
||
|
|
uInt dpd; // ..
|
||
|
|
uInt encode; // encoding accumulator
|
||
|
|
uInt sourhil, sourhir; // top words from source decFloats
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
#if QUAD
|
||
|
|
uShort uswork; // ..
|
||
|
|
#endif
|
||
|
|
// the following buffer holds the coefficient for manipulation
|
||
|
|
uByte buf[4+DECPMAX*3+2*QUAD]; // + space for zeros to left or right
|
||
|
|
#if DECTRACE
|
||
|
|
bcdnum num; // for trace displays
|
||
|
|
#endif
|
||
|
|
|
||
|
|
/* Start decoding the arguments */
|
||
|
|
sourhil=DFWORD(dfl, 0); // LHS top word
|
||
|
|
explb=DECCOMBEXP[sourhil>>26]; // get exponent high bits (in place)
|
||
|
|
sourhir=DFWORD(dfr, 0); // RHS top word
|
||
|
|
exprb=DECCOMBEXP[sourhir>>26];
|
||
|
|
|
||
|
|
if (EXPISSPECIAL(explb | exprb)) { // either is special?
|
||
|
|
// NaNs are handled as usual
|
||
|
|
if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
// one infinity but not both is bad
|
||
|
|
if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
|
||
|
|
// both infinite; return canonical infinity with sign of LHS
|
||
|
|
return decInfinity(result, dfl);
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Here when both arguments are finite */
|
||
|
|
// complete extraction of the exponents [no need to unbias]
|
||
|
|
explb+=GETECON(dfl); // + continuation
|
||
|
|
exprb+=GETECON(dfr); // ..
|
||
|
|
|
||
|
|
// calculate the number of digits to drop from the coefficient
|
||
|
|
drop=exprb-explb; // 0 if nothing to do
|
||
|
|
if (drop==0) return decCanonical(result, dfl); // return canonical
|
||
|
|
|
||
|
|
// the coefficient is needed; lay it out into buf, offset so zeros
|
||
|
|
// can be added before or after as needed -- an extra heading is
|
||
|
|
// added so can safely pad Quad DECPMAX-1 zeros to the left by
|
||
|
|
// fours
|
||
|
|
#define BUFOFF (buf+4+DECPMAX)
|
||
|
|
GETCOEFF(dfl, BUFOFF); // decode from decFloat
|
||
|
|
// [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1]
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
num.msd=BUFOFF;
|
||
|
|
num.lsd=BUFOFF+DECPMAX-1;
|
||
|
|
num.exponent=explb-DECBIAS;
|
||
|
|
num.sign=sourhil & DECFLOAT_Sign;
|
||
|
|
decShowNum(&num, "dfl");
|
||
|
|
#endif
|
||
|
|
|
||
|
|
if (drop>0) { // [most common case]
|
||
|
|
// (this code is very similar to that in decFloatFinalize, but
|
||
|
|
// has many differences so is duplicated here -- so any changes
|
||
|
|
// may need to be made there, too)
|
||
|
|
uByte *roundat; // -> re-round digit
|
||
|
|
uByte reround; // reround value
|
||
|
|
// printf("Rounding; drop=%ld\n", (LI)drop);
|
||
|
|
|
||
|
|
// there is at least one zero needed to the left, in all but one
|
||
|
|
// exceptional (all-nines) case, so place four zeros now; this is
|
||
|
|
// needed almost always and makes rounding all-nines by fours safe
|
||
|
|
UBFROMUI(BUFOFF-4, 0);
|
||
|
|
|
||
|
|
// Three cases here:
|
||
|
|
// 1. new LSD is in coefficient (almost always)
|
||
|
|
// 2. new LSD is digit to left of coefficient (so MSD is
|
||
|
|
// round-for-reround digit)
|
||
|
|
// 3. new LSD is to left of case 2 (whole coefficient is sticky)
|
||
|
|
// Note that leading zeros can safely be treated as useful digits
|
||
|
|
|
||
|
|
// [duplicate check-stickies code to save a test]
|
||
|
|
// [by-digit check for stickies as runs of zeros are rare]
|
||
|
|
if (drop<DECPMAX) { // NB lengths not addresses
|
||
|
|
roundat=BUFOFF+DECPMAX-drop;
|
||
|
|
reround=*roundat;
|
||
|
|
for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
|
||
|
|
if (*ub!=0) { // non-zero to be discarded
|
||
|
|
reround=DECSTICKYTAB[reround]; // apply sticky bit
|
||
|
|
break; // [remainder don't-care]
|
||
|
|
}
|
||
|
|
} // check stickies
|
||
|
|
ulsd=roundat-1; // set LSD
|
||
|
|
}
|
||
|
|
else { // edge case
|
||
|
|
if (drop==DECPMAX) {
|
||
|
|
roundat=BUFOFF;
|
||
|
|
reround=*roundat;
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
roundat=BUFOFF-1;
|
||
|
|
reround=0;
|
||
|
|
}
|
||
|
|
for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
|
||
|
|
if (*ub!=0) { // non-zero to be discarded
|
||
|
|
reround=DECSTICKYTAB[reround]; // apply sticky bit
|
||
|
|
break; // [remainder don't-care]
|
||
|
|
}
|
||
|
|
} // check stickies
|
||
|
|
*BUFOFF=0; // make a coefficient of 0
|
||
|
|
ulsd=BUFOFF; // .. at the MSD place
|
||
|
|
}
|
||
|
|
|
||
|
|
if (reround!=0) { // discarding non-zero
|
||
|
|
uInt bump=0;
|
||
|
|
set->status|=DEC_Inexact;
|
||
|
|
|
||
|
|
// next decide whether to increment the coefficient
|
||
|
|
if (set->round==DEC_ROUND_HALF_EVEN) { // fastpath slowest case
|
||
|
|
if (reround>5) bump=1; // >0.5 goes up
|
||
|
|
else if (reround==5) // exactly 0.5000 ..
|
||
|
|
bump=*ulsd & 0x01; // .. up iff [new] lsd is odd
|
||
|
|
} // r-h-e
|
||
|
|
else switch (set->round) {
|
||
|
|
case DEC_ROUND_DOWN: {
|
||
|
|
// no change
|
||
|
|
break;} // r-d
|
||
|
|
case DEC_ROUND_HALF_DOWN: {
|
||
|
|
if (reround>5) bump=1;
|
||
|
|
break;} // r-h-d
|
||
|
|
case DEC_ROUND_HALF_UP: {
|
||
|
|
if (reround>=5) bump=1;
|
||
|
|
break;} // r-h-u
|
||
|
|
case DEC_ROUND_UP: {
|
||
|
|
if (reround>0) bump=1;
|
||
|
|
break;} // r-u
|
||
|
|
case DEC_ROUND_CEILING: {
|
||
|
|
// same as _UP for positive numbers, and as _DOWN for negatives
|
||
|
|
if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
|
||
|
|
break;} // r-c
|
||
|
|
case DEC_ROUND_FLOOR: {
|
||
|
|
// same as _UP for negative numbers, and as _DOWN for positive
|
||
|
|
// [negative reround cannot occur on 0]
|
||
|
|
if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
|
||
|
|
break;} // r-f
|
||
|
|
case DEC_ROUND_05UP: {
|
||
|
|
if (reround>0) { // anything out there is 'sticky'
|
||
|
|
// bump iff lsd=0 or 5; this cannot carry so it could be
|
||
|
|
// effected immediately with no bump -- but the code
|
||
|
|
// is clearer if this is done the same way as the others
|
||
|
|
if (*ulsd==0 || *ulsd==5) bump=1;
|
||
|
|
}
|
||
|
|
break;} // r-r
|
||
|
|
default: { // e.g., DEC_ROUND_MAX
|
||
|
|
set->status|=DEC_Invalid_context;
|
||
|
|
#if DECCHECK
|
||
|
|
printf("Unknown rounding mode: %ld\n", (LI)set->round);
|
||
|
|
#endif
|
||
|
|
break;}
|
||
|
|
} // switch (not r-h-e)
|
||
|
|
// printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump);
|
||
|
|
|
||
|
|
if (bump!=0) { // need increment
|
||
|
|
// increment the coefficient; this could give 1000... (after
|
||
|
|
// the all nines case)
|
||
|
|
ub=ulsd;
|
||
|
|
for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
|
||
|
|
// now at most 3 digits left to non-9 (usually just the one)
|
||
|
|
for (; *ub==9; ub--) *ub=0;
|
||
|
|
*ub+=1;
|
||
|
|
// [the all-nines case will have carried one digit to the
|
||
|
|
// left of the original MSD -- just where it is needed]
|
||
|
|
} // bump needed
|
||
|
|
} // inexact rounding
|
||
|
|
|
||
|
|
// now clear zeros to the left so exactly DECPMAX digits will be
|
||
|
|
// available in the coefficent -- the first word to the left was
|
||
|
|
// cleared earlier for safe carry; now add any more needed
|
||
|
|
if (drop>4) {
|
||
|
|
UBFROMUI(BUFOFF-8, 0); // must be at least 5
|
||
|
|
for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 0);
|
||
|
|
}
|
||
|
|
} // need round (drop>0)
|
||
|
|
|
||
|
|
else { // drop<0; padding with -drop digits is needed
|
||
|
|
// This is the case where an error can occur if the padded
|
||
|
|
// coefficient will not fit; checking for this can be done in the
|
||
|
|
// same loop as padding for zeros if the no-hope and zero cases
|
||
|
|
// are checked first
|
||
|
|
if (-drop>DECPMAX-1) { // cannot fit unless 0
|
||
|
|
if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
|
||
|
|
// a zero can have any exponent; just drop through and use it
|
||
|
|
ulsd=BUFOFF+DECPMAX-1;
|
||
|
|
}
|
||
|
|
else { // padding will fit (but may still be too long)
|
||
|
|
// final-word mask depends on endianess
|
||
|
|
#if DECLITEND
|
||
|
|
static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
|
||
|
|
#else
|
||
|
|
static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
|
||
|
|
#endif
|
||
|
|
// note that here zeros to the right are added by fours, so in
|
||
|
|
// the Quad case this could write 36 zeros if the coefficient has
|
||
|
|
// fewer than three significant digits (hence the +2*QUAD for buf)
|
||
|
|
for (uc=BUFOFF+DECPMAX;; uc+=4) {
|
||
|
|
UBFROMUI(uc, 0);
|
||
|
|
if (UBTOUI(uc-DECPMAX)!=0) { // could be bad
|
||
|
|
// if all four digits should be zero, definitely bad
|
||
|
|
if (uc<=BUFOFF+DECPMAX+(-drop)-4)
|
||
|
|
return decInvalid(result, set);
|
||
|
|
// must be a 1- to 3-digit sequence; check more carefully
|
||
|
|
if ((UBTOUI(uc-DECPMAX)&dmask[(-drop)%4])!=0)
|
||
|
|
return decInvalid(result, set);
|
||
|
|
break; // no need for loop end test
|
||
|
|
}
|
||
|
|
if (uc>=BUFOFF+DECPMAX+(-drop)-4) break; // done
|
||
|
|
}
|
||
|
|
ulsd=BUFOFF+DECPMAX+(-drop)-1;
|
||
|
|
} // pad and check leading zeros
|
||
|
|
} // drop<0
|
||
|
|
|
||
|
|
#if DECTRACE
|
||
|
|
num.msd=ulsd-DECPMAX+1;
|
||
|
|
num.lsd=ulsd;
|
||
|
|
num.exponent=explb-DECBIAS;
|
||
|
|
num.sign=sourhil & DECFLOAT_Sign;
|
||
|
|
decShowNum(&num, "res");
|
||
|
|
#endif
|
||
|
|
|
||
|
|
/*------------------------------------------------------------------*/
|
||
|
|
/* At this point the result is DECPMAX digits, ending at ulsd, so */
|
||
|
|
/* fits the encoding exactly; there is no possibility of error */
|
||
|
|
/*------------------------------------------------------------------*/
|
||
|
|
encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); // make index
|
||
|
|
encode=DECCOMBFROM[encode]; // indexed by (0-2)*16+msd
|
||
|
|
// the exponent continuation can be extracted from the original RHS
|
||
|
|
encode|=sourhir & ECONMASK;
|
||
|
|
encode|=sourhil&DECFLOAT_Sign; // add the sign from LHS
|
||
|
|
|
||
|
|
// finally encode the coefficient
|
||
|
|
// private macro to encode a declet; this version can be used
|
||
|
|
// because all coefficient digits exist
|
||
|
|
#define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \
|
||
|
|
dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
|
||
|
|
|
||
|
|
#if DOUBLE
|
||
|
|
getDPD3q(dpd, 4); encode|=dpd<<8;
|
||
|
|
getDPD3q(dpd, 3); encode|=dpd>>2;
|
||
|
|
DFWORD(result, 0)=encode;
|
||
|
|
encode=dpd<<30;
|
||
|
|
getDPD3q(dpd, 2); encode|=dpd<<20;
|
||
|
|
getDPD3q(dpd, 1); encode|=dpd<<10;
|
||
|
|
getDPD3q(dpd, 0); encode|=dpd;
|
||
|
|
DFWORD(result, 1)=encode;
|
||
|
|
|
||
|
|
#elif QUAD
|
||
|
|
getDPD3q(dpd,10); encode|=dpd<<4;
|
||
|
|
getDPD3q(dpd, 9); encode|=dpd>>6;
|
||
|
|
DFWORD(result, 0)=encode;
|
||
|
|
encode=dpd<<26;
|
||
|
|
getDPD3q(dpd, 8); encode|=dpd<<16;
|
||
|
|
getDPD3q(dpd, 7); encode|=dpd<<6;
|
||
|
|
getDPD3q(dpd, 6); encode|=dpd>>4;
|
||
|
|
DFWORD(result, 1)=encode;
|
||
|
|
encode=dpd<<28;
|
||
|
|
getDPD3q(dpd, 5); encode|=dpd<<18;
|
||
|
|
getDPD3q(dpd, 4); encode|=dpd<<8;
|
||
|
|
getDPD3q(dpd, 3); encode|=dpd>>2;
|
||
|
|
DFWORD(result, 2)=encode;
|
||
|
|
encode=dpd<<30;
|
||
|
|
getDPD3q(dpd, 2); encode|=dpd<<20;
|
||
|
|
getDPD3q(dpd, 1); encode|=dpd<<10;
|
||
|
|
getDPD3q(dpd, 0); encode|=dpd;
|
||
|
|
DFWORD(result, 3)=encode;
|
||
|
|
#endif
|
||
|
|
return result;
|
||
|
|
} // decFloatQuantize
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatReduce -- reduce finite coefficient to minimum length */
|
||
|
|
/* */
|
||
|
|
/* result gets the reduced decFloat */
|
||
|
|
/* df is the source decFloat */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result, which will be canonical */
|
||
|
|
/* */
|
||
|
|
/* This removes all possible trailing zeros from the coefficient; */
|
||
|
|
/* some may remain when the number is very close to Nmax. */
|
||
|
|
/* Special values are unchanged and no status is set unless df=sNaN. */
|
||
|
|
/* Reduced zero has an exponent q=0. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatReduce(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set) {
|
||
|
|
bcdnum num; // work
|
||
|
|
uByte buf[DECPMAX], *ub; // coefficient and pointer
|
||
|
|
if (df!=result) *result=*df; // copy, if needed
|
||
|
|
if (DFISNAN(df)) return decNaNs(result, df, NULL, set); // sNaN
|
||
|
|
// zeros and infinites propagate too
|
||
|
|
if (DFISINF(df)) return decInfinity(result, df); // canonical
|
||
|
|
if (DFISZERO(df)) {
|
||
|
|
uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
|
||
|
|
decFloatZero(result);
|
||
|
|
DFWORD(result, 0)|=sign;
|
||
|
|
return result; // exponent dropped, sign OK
|
||
|
|
}
|
||
|
|
// non-zero finite
|
||
|
|
GETCOEFF(df, buf);
|
||
|
|
ub=buf+DECPMAX-1; // -> lsd
|
||
|
|
if (*ub) return result; // no trailing zeros
|
||
|
|
for (ub--; *ub==0;) ub--; // terminates because non-zero
|
||
|
|
// *ub is the first non-zero from the right
|
||
|
|
num.sign=DFWORD(df, 0)&DECFLOAT_Sign; // set up number...
|
||
|
|
num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); // adjusted exponent
|
||
|
|
num.msd=buf;
|
||
|
|
num.lsd=ub;
|
||
|
|
return decFinalize(result, &num, set);
|
||
|
|
} // decFloatReduce
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatRemainder -- integer divide and return remainder */
|
||
|
|
/* */
|
||
|
|
/* result gets the remainder of dividing dfl by dfr: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatRemainder(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
return decDivide(result, dfl, dfr, set, REMAINDER);
|
||
|
|
} // decFloatRemainder
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatRemainderNear -- integer divide to nearest and remainder */
|
||
|
|
/* */
|
||
|
|
/* result gets the remainder of dividing dfl by dfr: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This is the IEEE remainder, where the nearest integer is used. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatRemainderNear(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
return decDivide(result, dfl, dfr, set, REMNEAR);
|
||
|
|
} // decFloatRemainderNear
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatRotate -- rotate the coefficient of a decFloat left/right */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of rotating dfl */
|
||
|
|
/* dfl is the source decFloat to rotate */
|
||
|
|
/* dfr is the count of digits to rotate, an integer (with q=0) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* The digits of the coefficient of dfl are rotated to the left (if */
|
||
|
|
/* dfr is positive) or to the right (if dfr is negative) without */
|
||
|
|
/* adjusting the exponent or the sign of dfl. */
|
||
|
|
/* */
|
||
|
|
/* dfr must be in the range -DECPMAX through +DECPMAX. */
|
||
|
|
/* NaNs are propagated as usual. An infinite dfl is unaffected (but */
|
||
|
|
/* dfr must be valid). No status is set unless dfr is invalid or an */
|
||
|
|
/* operand is an sNaN. The result is canonical. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
#define PHALF (ROUNDUP(DECPMAX/2, 4)) // half length, rounded up
|
||
|
|
decFloat * decFloatRotate(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int rotate; // dfr as an Int
|
||
|
|
uByte buf[DECPMAX+PHALF]; // coefficient + half
|
||
|
|
uInt digits, savestat; // work
|
||
|
|
bcdnum num; // ..
|
||
|
|
uByte *ub; // ..
|
||
|
|
|
||
|
|
if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
if (!DFISINT(dfr)) return decInvalid(result, set);
|
||
|
|
digits=decFloatDigits(dfr); // calculate digits
|
||
|
|
if (digits>2) return decInvalid(result, set); // definitely out of range
|
||
|
|
rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; // is in bottom declet
|
||
|
|
if (rotate>DECPMAX) return decInvalid(result, set); // too big
|
||
|
|
// [from here on no error or status change is possible]
|
||
|
|
if (DFISINF(dfl)) return decInfinity(result, dfl); // canonical
|
||
|
|
// handle no-rotate cases
|
||
|
|
if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
|
||
|
|
// a real rotate is needed: 0 < rotate < DECPMAX
|
||
|
|
// reduce the rotation to no more than half to reduce copying later
|
||
|
|
// (for QUAD in fact half + 2 digits)
|
||
|
|
if (DFISSIGNED(dfr)) rotate=-rotate;
|
||
|
|
if (abs(rotate)>PHALF) {
|
||
|
|
if (rotate<0) rotate=DECPMAX+rotate;
|
||
|
|
else rotate=rotate-DECPMAX;
|
||
|
|
}
|
||
|
|
// now lay out the coefficient, leaving room to the right or the
|
||
|
|
// left depending on the direction of rotation
|
||
|
|
ub=buf;
|
||
|
|
if (rotate<0) ub+=PHALF; // rotate right, so space to left
|
||
|
|
GETCOEFF(dfl, ub);
|
||
|
|
// copy half the digits to left or right, and set num.msd
|
||
|
|
if (rotate<0) {
|
||
|
|
memcpy(buf, buf+DECPMAX, PHALF);
|
||
|
|
num.msd=buf+PHALF+rotate;
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
memcpy(buf+DECPMAX, buf, PHALF);
|
||
|
|
num.msd=buf+rotate;
|
||
|
|
}
|
||
|
|
// fill in rest of num
|
||
|
|
num.lsd=num.msd+DECPMAX-1;
|
||
|
|
num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
|
||
|
|
num.exponent=GETEXPUN(dfl);
|
||
|
|
savestat=set->status; // record
|
||
|
|
decFinalize(result, &num, set);
|
||
|
|
set->status=savestat; // restore
|
||
|
|
return result;
|
||
|
|
} // decFloatRotate
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatSameQuantum -- test decFloats for same quantum */
|
||
|
|
/* */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* returns 1 if the operands have the same quantum, 0 otherwise */
|
||
|
|
/* */
|
||
|
|
/* No error is possible and no status results. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
|
||
|
|
if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
|
||
|
|
if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
|
||
|
|
if (DFISINF(dfl) && DFISINF(dfr)) return 1;
|
||
|
|
return 0; // any other special mixture gives false
|
||
|
|
}
|
||
|
|
if (GETEXP(dfl)==GETEXP(dfr)) return 1; // biased exponents match
|
||
|
|
return 0;
|
||
|
|
} // decFloatSameQuantum
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatScaleB -- multiply by a power of 10, as per 754 */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of the operation */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs), am integer (with q=0) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* This computes result=dfl x 10**dfr where dfr is an integer in the */
|
||
|
|
/* range +/-2*(emax+pmax), typically resulting from LogB. */
|
||
|
|
/* Underflow and Overflow (with Inexact) may occur. NaNs propagate */
|
||
|
|
/* as usual. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
#define SCALEBMAX 2*(DECEMAX+DECPMAX) // D=800, Q=12356
|
||
|
|
decFloat * decFloatScaleB(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
uInt digits; // work
|
||
|
|
Int expr; // dfr as an Int
|
||
|
|
|
||
|
|
if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
if (!DFISINT(dfr)) return decInvalid(result, set);
|
||
|
|
digits=decFloatDigits(dfr); // calculate digits
|
||
|
|
|
||
|
|
#if DOUBLE
|
||
|
|
if (digits>3) return decInvalid(result, set); // definitely out of range
|
||
|
|
expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff]; // must be in bottom declet
|
||
|
|
#elif QUAD
|
||
|
|
if (digits>5) return decInvalid(result, set); // definitely out of range
|
||
|
|
expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff] // in bottom 2 declets ..
|
||
|
|
+DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000; // ..
|
||
|
|
#endif
|
||
|
|
if (expr>SCALEBMAX) return decInvalid(result, set); // oops
|
||
|
|
// [from now on no error possible]
|
||
|
|
if (DFISINF(dfl)) return decInfinity(result, dfl); // canonical
|
||
|
|
if (DFISSIGNED(dfr)) expr=-expr;
|
||
|
|
// dfl is finite and expr is valid
|
||
|
|
*result=*dfl; // copy to target
|
||
|
|
return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
|
||
|
|
} // decFloatScaleB
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatShift -- shift the coefficient of a decFloat left or right */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of shifting dfl */
|
||
|
|
/* dfl is the source decFloat to shift */
|
||
|
|
/* dfr is the count of digits to shift, an integer (with q=0) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* The digits of the coefficient of dfl are shifted to the left (if */
|
||
|
|
/* dfr is positive) or to the right (if dfr is negative) without */
|
||
|
|
/* adjusting the exponent or the sign of dfl. */
|
||
|
|
/* */
|
||
|
|
/* dfr must be in the range -DECPMAX through +DECPMAX. */
|
||
|
|
/* NaNs are propagated as usual. An infinite dfl is unaffected (but */
|
||
|
|
/* dfr must be valid). No status is set unless dfr is invalid or an */
|
||
|
|
/* operand is an sNaN. The result is canonical. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatShift(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
Int shift; // dfr as an Int
|
||
|
|
uByte buf[DECPMAX*2]; // coefficient + padding
|
||
|
|
uInt digits, savestat; // work
|
||
|
|
bcdnum num; // ..
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
|
||
|
|
if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
|
||
|
|
if (!DFISINT(dfr)) return decInvalid(result, set);
|
||
|
|
digits=decFloatDigits(dfr); // calculate digits
|
||
|
|
if (digits>2) return decInvalid(result, set); // definitely out of range
|
||
|
|
shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; // is in bottom declet
|
||
|
|
if (shift>DECPMAX) return decInvalid(result, set); // too big
|
||
|
|
// [from here on no error or status change is possible]
|
||
|
|
|
||
|
|
if (DFISINF(dfl)) return decInfinity(result, dfl); // canonical
|
||
|
|
// handle no-shift and all-shift (clear to zero) cases
|
||
|
|
if (shift==0) return decCanonical(result, dfl);
|
||
|
|
if (shift==DECPMAX) { // zero with sign
|
||
|
|
uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); // save sign bit
|
||
|
|
decFloatZero(result); // make +0
|
||
|
|
DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); // and set sign
|
||
|
|
// [cannot safely use CopySign]
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
// a real shift is needed: 0 < shift < DECPMAX
|
||
|
|
num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
|
||
|
|
num.exponent=GETEXPUN(dfl);
|
||
|
|
num.msd=buf;
|
||
|
|
GETCOEFF(dfl, buf);
|
||
|
|
if (DFISSIGNED(dfr)) { // shift right
|
||
|
|
// edge cases are taken care of, so this is easy
|
||
|
|
num.lsd=buf+DECPMAX-shift-1;
|
||
|
|
}
|
||
|
|
else { // shift left -- zero padding needed to right
|
||
|
|
UBFROMUI(buf+DECPMAX, 0); // 8 will handle most cases
|
||
|
|
UBFROMUI(buf+DECPMAX+4, 0); // ..
|
||
|
|
if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); // all other cases
|
||
|
|
num.msd+=shift;
|
||
|
|
num.lsd=num.msd+DECPMAX-1;
|
||
|
|
}
|
||
|
|
savestat=set->status; // record
|
||
|
|
decFinalize(result, &num, set);
|
||
|
|
set->status=savestat; // restore
|
||
|
|
return result;
|
||
|
|
} // decFloatShift
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatSubtract -- subtract a decFloat from another */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of subtracting dfr from dfl: */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatSubtract(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
decFloat temp;
|
||
|
|
// NaNs must propagate without sign change
|
||
|
|
if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
|
||
|
|
temp=*dfr; // make a copy
|
||
|
|
DFBYTE(&temp, 0)^=0x80; // flip sign
|
||
|
|
return decFloatAdd(result, dfl, &temp, set); // and add to the lhs
|
||
|
|
} // decFloatSubtract
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
|
||
|
|
/* */
|
||
|
|
/* df is the decFloat to round */
|
||
|
|
/* set is the context */
|
||
|
|
/* round is the rounding mode to use */
|
||
|
|
/* returns a uInt or an Int, rounded according to the name */
|
||
|
|
/* */
|
||
|
|
/* Invalid will always be signaled if df is a NaN, is Infinite, or is */
|
||
|
|
/* outside the range of the target; Inexact will not be signaled for */
|
||
|
|
/* simple rounding unless 'Exact' appears in the name. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
uInt decFloatToUInt32(const decFloat *df, decContext *set,
|
||
|
|
enum rounding round) {
|
||
|
|
return decToInt32(df, set, round, 0, 1);}
|
||
|
|
|
||
|
|
uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
|
||
|
|
enum rounding round) {
|
||
|
|
return decToInt32(df, set, round, 1, 1);}
|
||
|
|
|
||
|
|
Int decFloatToInt32(const decFloat *df, decContext *set,
|
||
|
|
enum rounding round) {
|
||
|
|
return (Int)decToInt32(df, set, round, 0, 0);}
|
||
|
|
|
||
|
|
Int decFloatToInt32Exact(const decFloat *df, decContext *set,
|
||
|
|
enum rounding round) {
|
||
|
|
return (Int)decToInt32(df, set, round, 1, 0);}
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatToIntegral -- round to integral value (two flavours) */
|
||
|
|
/* */
|
||
|
|
/* result gets the result */
|
||
|
|
/* df is the decFloat to round */
|
||
|
|
/* set is the context */
|
||
|
|
/* round is the rounding mode to use */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* No exceptions, even Inexact, are raised except for sNaN input, or */
|
||
|
|
/* if 'Exact' appears in the name. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set, enum rounding round) {
|
||
|
|
return decToIntegral(result, df, set, round, 0);}
|
||
|
|
|
||
|
|
decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set) {
|
||
|
|
return decToIntegral(result, df, set, set->round, 1);}
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decFloatXor -- logical digitwise XOR of two decFloats */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of XORing dfl and dfr */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result, which will be canonical with sign=0 */
|
||
|
|
/* */
|
||
|
|
/* The operands must be positive, finite with exponent q=0, and */
|
||
|
|
/* comprise just zeros and ones; if not, Invalid operation results. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
decFloat * decFloatXor(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
|
||
|
|
|| !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
|
||
|
|
// the operands are positive finite integers (q=0) with just 0s and 1s
|
||
|
|
#if DOUBLE
|
||
|
|
DFWORD(result, 0)=ZEROWORD
|
||
|
|
|((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
|
||
|
|
DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
|
||
|
|
#elif QUAD
|
||
|
|
DFWORD(result, 0)=ZEROWORD
|
||
|
|
|((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
|
||
|
|
DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
|
||
|
|
DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
|
||
|
|
DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
|
||
|
|
#endif
|
||
|
|
return result;
|
||
|
|
} // decFloatXor
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decInvalid -- set Invalid_operation result */
|
||
|
|
/* */
|
||
|
|
/* result gets a canonical NaN */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* status has Invalid_operation added */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
static decFloat *decInvalid(decFloat *result, decContext *set) {
|
||
|
|
decFloatZero(result);
|
||
|
|
DFWORD(result, 0)=DECFLOAT_qNaN;
|
||
|
|
set->status|=DEC_Invalid_operation;
|
||
|
|
return result;
|
||
|
|
} // decInvalid
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decInfinity -- set canonical Infinity with sign from a decFloat */
|
||
|
|
/* */
|
||
|
|
/* result gets a canonical Infinity */
|
||
|
|
/* df is source decFloat (only the sign is used) */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* df may be the same as result */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
static decFloat *decInfinity(decFloat *result, const decFloat *df) {
|
||
|
|
uInt sign=DFWORD(df, 0); // save source signword
|
||
|
|
decFloatZero(result); // clear everything
|
||
|
|
DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
|
||
|
|
return result;
|
||
|
|
} // decInfinity
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decNaNs -- handle NaN argument(s) */
|
||
|
|
/* */
|
||
|
|
/* result gets the result of handling dfl and dfr, one or both of */
|
||
|
|
/* which is a NaN */
|
||
|
|
/* dfl is the first decFloat (lhs) */
|
||
|
|
/* dfr is the second decFloat (rhs) -- may be NULL for a single- */
|
||
|
|
/* operand operation */
|
||
|
|
/* set is the context */
|
||
|
|
/* returns result */
|
||
|
|
/* */
|
||
|
|
/* Called when one or both operands is a NaN, and propagates the */
|
||
|
|
/* appropriate result to res. When an sNaN is found, it is changed */
|
||
|
|
/* to a qNaN and Invalid operation is set. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
static decFloat *decNaNs(decFloat *result,
|
||
|
|
const decFloat *dfl, const decFloat *dfr,
|
||
|
|
decContext *set) {
|
||
|
|
// handle sNaNs first
|
||
|
|
if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; // use RHS
|
||
|
|
if (DFISSNAN(dfl)) {
|
||
|
|
decCanonical(result, dfl); // propagate canonical sNaN
|
||
|
|
DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); // quiet
|
||
|
|
set->status|=DEC_Invalid_operation;
|
||
|
|
return result;
|
||
|
|
}
|
||
|
|
// one or both is a quiet NaN
|
||
|
|
if (!DFISNAN(dfl)) dfl=dfr; // RHS must be NaN, use it
|
||
|
|
return decCanonical(result, dfl); // propagate canonical qNaN
|
||
|
|
} // decNaNs
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decNumCompare -- numeric comparison of two decFloats */
|
||
|
|
/* */
|
||
|
|
/* dfl is the left-hand decFloat, which is not a NaN */
|
||
|
|
/* dfr is the right-hand decFloat, which is not a NaN */
|
||
|
|
/* tot is 1 for total order compare, 0 for simple numeric */
|
||
|
|
/* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
|
||
|
|
/* */
|
||
|
|
/* No error is possible; status and mode are unchanged. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
|
||
|
|
Int sigl, sigr; // LHS and RHS non-0 signums
|
||
|
|
Int shift; // shift needed to align operands
|
||
|
|
uByte *ub, *uc; // work
|
||
|
|
uInt uiwork; // for macros
|
||
|
|
// buffers +2 if Quad (36 digits), need double plus 4 for safe padding
|
||
|
|
uByte bufl[DECPMAX*2+QUAD*2+4]; // for LHS coefficient + padding
|
||
|
|
uByte bufr[DECPMAX*2+QUAD*2+4]; // for RHS coefficient + padding
|
||
|
|
|
||
|
|
sigl=1;
|
||
|
|
if (DFISSIGNED(dfl)) {
|
||
|
|
if (!DFISSIGNED(dfr)) { // -LHS +RHS
|
||
|
|
if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
|
||
|
|
return -1; // RHS wins
|
||
|
|
}
|
||
|
|
sigl=-1;
|
||
|
|
}
|
||
|
|
if (DFISSIGNED(dfr)) {
|
||
|
|
if (!DFISSIGNED(dfl)) { // +LHS -RHS
|
||
|
|
if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
|
||
|
|
return +1; // LHS wins
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// signs are the same; operand(s) could be zero
|
||
|
|
sigr=-sigl; // sign to return if abs(RHS) wins
|
||
|
|
|
||
|
|
if (DFISINF(dfl)) {
|
||
|
|
if (DFISINF(dfr)) return 0; // both infinite & same sign
|
||
|
|
return sigl; // inf > n
|
||
|
|
}
|
||
|
|
if (DFISINF(dfr)) return sigr; // n < inf [dfl is finite]
|
||
|
|
|
||
|
|
// here, both are same sign and finite; calculate their offset
|
||
|
|
shift=GETEXP(dfl)-GETEXP(dfr); // [0 means aligned]
|
||
|
|
// [bias can be ignored -- the absolute exponent is not relevant]
|
||
|
|
|
||
|
|
if (DFISZERO(dfl)) {
|
||
|
|
if (!DFISZERO(dfr)) return sigr; // LHS=0, RHS!=0
|
||
|
|
// both are zero, return 0 if both same exponent or numeric compare
|
||
|
|
if (shift==0 || !tot) return 0;
|
||
|
|
if (shift>0) return sigl;
|
||
|
|
return sigr; // [shift<0]
|
||
|
|
}
|
||
|
|
else { // LHS!=0
|
||
|
|
if (DFISZERO(dfr)) return sigl; // LHS!=0, RHS=0
|
||
|
|
}
|
||
|
|
// both are known to be non-zero at this point
|
||
|
|
|
||
|
|
// if the exponents are so different that the coefficients do not
|
||
|
|
// overlap (by even one digit) then a full comparison is not needed
|
||
|
|
if (abs(shift)>=DECPMAX) { // no overlap
|
||
|
|
// coefficients are known to be non-zero
|
||
|
|
if (shift>0) return sigl;
|
||
|
|
return sigr; // [shift<0]
|
||
|
|
}
|
||
|
|
|
||
|
|
// decode the coefficients
|
||
|
|
// (shift both right two if Quad to make a multiple of four)
|
||
|
|
#if QUAD
|
||
|
|
UBFROMUI(bufl, 0);
|
||
|
|
UBFROMUI(bufr, 0);
|
||
|
|
#endif
|
||
|
|
GETCOEFF(dfl, bufl+QUAD*2); // decode from decFloat
|
||
|
|
GETCOEFF(dfr, bufr+QUAD*2); // ..
|
||
|
|
if (shift==0) { // aligned; common and easy
|
||
|
|
// all multiples of four, here
|
||
|
|
for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
|
||
|
|
uInt ui=UBTOUI(ub);
|
||
|
|
if (ui==UBTOUI(uc)) continue; // so far so same
|
||
|
|
// about to find a winner; go by bytes in case little-endian
|
||
|
|
for (;; ub++, uc++) {
|
||
|
|
if (*ub>*uc) return sigl; // difference found
|
||
|
|
if (*ub<*uc) return sigr; // ..
|
||
|
|
}
|
||
|
|
}
|
||
|
|
} // aligned
|
||
|
|
else if (shift>0) { // lhs to left
|
||
|
|
ub=bufl; // RHS pointer
|
||
|
|
// pad bufl so right-aligned; most shifts will fit in 8
|
||
|
|
UBFROMUI(bufl+DECPMAX+QUAD*2, 0); // add eight zeros
|
||
|
|
UBFROMUI(bufl+DECPMAX+QUAD*2+4, 0); // ..
|
||
|
|
if (shift>8) {
|
||
|
|
// more than eight; fill the rest, and also worth doing the
|
||
|
|
// lead-in by fours
|
||
|
|
uByte *up; // work
|
||
|
|
uByte *upend=bufl+DECPMAX+QUAD*2+shift;
|
||
|
|
for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
|
||
|
|
// [pads up to 36 in all for Quad]
|
||
|
|
for (;; ub+=4) {
|
||
|
|
if (UBTOUI(ub)!=0) return sigl;
|
||
|
|
if (ub+4>bufl+shift-4) break;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
// check remaining leading digits
|
||
|
|
for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
|
||
|
|
// now start the overlapped part; bufl has been padded, so the
|
||
|
|
// comparison can go for the full length of bufr, which is a
|
||
|
|
// multiple of 4 bytes
|
||
|
|
for (uc=bufr; ; uc+=4, ub+=4) {
|
||
|
|
uInt ui=UBTOUI(ub);
|
||
|
|
if (ui!=UBTOUI(uc)) { // mismatch found
|
||
|
|
for (;; uc++, ub++) { // check from left [little-endian?]
|
||
|
|
if (*ub>*uc) return sigl; // difference found
|
||
|
|
if (*ub<*uc) return sigr; // ..
|
||
|
|
}
|
||
|
|
} // mismatch
|
||
|
|
if (uc==bufr+QUAD*2+DECPMAX-4) break; // all checked
|
||
|
|
}
|
||
|
|
} // shift>0
|
||
|
|
|
||
|
|
else { // shift<0) .. RHS is to left of LHS; mirror shift>0
|
||
|
|
uc=bufr; // RHS pointer
|
||
|
|
// pad bufr so right-aligned; most shifts will fit in 8
|
||
|
|
UBFROMUI(bufr+DECPMAX+QUAD*2, 0); // add eight zeros
|
||
|
|
UBFROMUI(bufr+DECPMAX+QUAD*2+4, 0); // ..
|
||
|
|
if (shift<-8) {
|
||
|
|
// more than eight; fill the rest, and also worth doing the
|
||
|
|
// lead-in by fours
|
||
|
|
uByte *up; // work
|
||
|
|
uByte *upend=bufr+DECPMAX+QUAD*2-shift;
|
||
|
|
for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
|
||
|
|
// [pads up to 36 in all for Quad]
|
||
|
|
for (;; uc+=4) {
|
||
|
|
if (UBTOUI(uc)!=0) return sigr;
|
||
|
|
if (uc+4>bufr-shift-4) break;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
// check remaining leading digits
|
||
|
|
for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
|
||
|
|
// now start the overlapped part; bufr has been padded, so the
|
||
|
|
// comparison can go for the full length of bufl, which is a
|
||
|
|
// multiple of 4 bytes
|
||
|
|
for (ub=bufl; ; ub+=4, uc+=4) {
|
||
|
|
uInt ui=UBTOUI(ub);
|
||
|
|
if (ui!=UBTOUI(uc)) { // mismatch found
|
||
|
|
for (;; ub++, uc++) { // check from left [little-endian?]
|
||
|
|
if (*ub>*uc) return sigl; // difference found
|
||
|
|
if (*ub<*uc) return sigr; // ..
|
||
|
|
}
|
||
|
|
} // mismatch
|
||
|
|
if (ub==bufl+QUAD*2+DECPMAX-4) break; // all checked
|
||
|
|
}
|
||
|
|
} // shift<0
|
||
|
|
|
||
|
|
// Here when compare equal
|
||
|
|
if (!tot) return 0; // numerically equal
|
||
|
|
// total ordering .. exponent matters
|
||
|
|
if (shift>0) return sigl; // total order by exponent
|
||
|
|
if (shift<0) return sigr; // ..
|
||
|
|
return 0;
|
||
|
|
} // decNumCompare
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decToInt32 -- local routine to effect ToInteger conversions */
|
||
|
|
/* */
|
||
|
|
/* df is the decFloat to convert */
|
||
|
|
/* set is the context */
|
||
|
|
/* rmode is the rounding mode to use */
|
||
|
|
/* exact is 1 if Inexact should be signalled */
|
||
|
|
/* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */
|
||
|
|
/* returns 32-bit result as a uInt */
|
||
|
|
/* */
|
||
|
|
/* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
|
||
|
|
/* these cases 0 is returned. */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
static uInt decToInt32(const decFloat *df, decContext *set,
|
||
|
|
enum rounding rmode, Flag exact, Flag unsign) {
|
||
|
|
Int exp; // exponent
|
||
|
|
uInt sourhi, sourpen, sourlo; // top word from source decFloat ..
|
||
|
|
uInt hi, lo; // .. penultimate, least, etc.
|
||
|
|
decFloat zero, result; // work
|
||
|
|
Int i; // ..
|
||
|
|
|
||
|
|
/* Start decoding the argument */
|
||
|
|
sourhi=DFWORD(df, 0); // top word
|
||
|
|
exp=DECCOMBEXP[sourhi>>26]; // get exponent high bits (in place)
|
||
|
|
if (EXPISSPECIAL(exp)) { // is special?
|
||
|
|
set->status|=DEC_Invalid_operation; // signal
|
||
|
|
return 0;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Here when the argument is finite */
|
||
|
|
if (GETEXPUN(df)==0) result=*df; // already a true integer
|
||
|
|
else { // need to round to integer
|
||
|
|
enum rounding saveround; // saver
|
||
|
|
uInt savestatus; // ..
|
||
|
|
saveround=set->round; // save rounding mode ..
|
||
|
|
savestatus=set->status; // .. and status
|
||
|
|
set->round=rmode; // set mode
|
||
|
|
decFloatZero(&zero); // make 0E+0
|
||
|
|
set->status=0; // clear
|
||
|
|
decFloatQuantize(&result, df, &zero, set); // [this may fail]
|
||
|
|
set->round=saveround; // restore rounding mode ..
|
||
|
|
if (exact) set->status|=savestatus; // include Inexact
|
||
|
|
else set->status=savestatus; // .. or just original status
|
||
|
|
}
|
||
|
|
|
||
|
|
// only the last four declets of the coefficient can contain
|
||
|
|
// non-zero; check for others (and also NaN or Infinity from the
|
||
|
|
// Quantize) first (see DFISZERO for explanation):
|
||
|
|
// decFloatShow(&result, "sofar");
|
||
|
|
#if DOUBLE
|
||
|
|
if ((DFWORD(&result, 0)&0x1c03ff00)!=0
|
||
|
|
|| (DFWORD(&result, 0)&0x60000000)==0x60000000) {
|
||
|
|
#elif QUAD
|
||
|
|
if ((DFWORD(&result, 2)&0xffffff00)!=0
|
||
|
|
|| DFWORD(&result, 1)!=0
|
||
|
|
|| (DFWORD(&result, 0)&0x1c003fff)!=0
|
||
|
|
|| (DFWORD(&result, 0)&0x60000000)==0x60000000) {
|
||
|
|
#endif
|
||
|
|
set->status|=DEC_Invalid_operation; // Invalid or out of range
|
||
|
|
return 0;
|
||
|
|
}
|
||
|
|
// get last twelve digits of the coefficent into hi & ho, base
|
||
|
|
// 10**9 (see GETCOEFFBILL):
|
||
|
|
sourlo=DFWORD(&result, DECWORDS-1);
|
||
|
|
lo=DPD2BIN0[sourlo&0x3ff]
|
||
|
|
+DPD2BINK[(sourlo>>10)&0x3ff]
|
||
|
|
+DPD2BINM[(sourlo>>20)&0x3ff];
|
||
|
|
sourpen=DFWORD(&result, DECWORDS-2);
|
||
|
|
hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
|
||
|
|
|
||
|
|
// according to request, check range carefully
|
||
|
|
if (unsign) {
|
||
|
|
if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
|
||
|
|
set->status|=DEC_Invalid_operation; // out of range
|
||
|
|
return 0;
|
||
|
|
}
|
||
|
|
return hi*BILLION+lo;
|
||
|
|
}
|
||
|
|
// signed
|
||
|
|
if (hi>2 || (hi==2 && lo>147483647)) {
|
||
|
|
// handle the usual edge case
|
||
|
|
if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
|
||
|
|
set->status|=DEC_Invalid_operation; // truly out of range
|
||
|
|
return 0;
|
||
|
|
}
|
||
|
|
i=hi*BILLION+lo;
|
||
|
|
if (DFISSIGNED(&result)) i=-i;
|
||
|
|
return (uInt)i;
|
||
|
|
} // decToInt32
|
||
|
|
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
/* decToIntegral -- local routine to effect ToIntegral value */
|
||
|
|
/* */
|
||
|
|
/* result gets the result */
|
||
|
|
/* df is the decFloat to round */
|
||
|
|
/* set is the context */
|
||
|
|
/* rmode is the rounding mode to use */
|
||
|
|
/* exact is 1 if Inexact should be signalled */
|
||
|
|
/* returns result */
|
||
|
|
/* ------------------------------------------------------------------ */
|
||
|
|
static decFloat * decToIntegral(decFloat *result, const decFloat *df,
|
||
|
|
decContext *set, enum rounding rmode,
|
||
|
|
Flag exact) {
|
||
|
|
Int exp; // exponent
|
||
|
|
uInt sourhi; // top word from source decFloat
|
||
|
|
enum rounding saveround; // saver
|
||
|
|
uInt savestatus; // ..
|
||
|
|
decFloat zero; // work
|
||
|
|
|
||
|
|
/* Start decoding the argument */
|
||
|
|
sourhi=DFWORD(df, 0); // top word
|
||
|
|
exp=DECCOMBEXP[sourhi>>26]; // get exponent high bits (in place)
|
||
|
|
|
||
|
|
if (EXPISSPECIAL(exp)) { // is special?
|
||
|
|
// NaNs are handled as usual
|
||
|
|
if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
|
||
|
|
// must be infinite; return canonical infinity with sign of df
|
||
|
|
return decInfinity(result, df);
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Here when the argument is finite */
|
||
|
|
// complete extraction of the exponent
|
||
|
|
exp+=GETECON(df)-DECBIAS; // .. + continuation and unbias
|
||
|
|
|
||
|
|
if (exp>=0) return decCanonical(result, df); // already integral
|
||
|
|
|
||
|
|
saveround=set->round; // save rounding mode ..
|
||
|
|
savestatus=set->status; // .. and status
|
||
|
|
set->round=rmode; // set mode
|
||
|
|
decFloatZero(&zero); // make 0E+0
|
||
|
|
decFloatQuantize(result, df, &zero, set); // 'integrate'; cannot fail
|
||
|
|
set->round=saveround; // restore rounding mode ..
|
||
|
|
if (!exact) set->status=savestatus; // .. and status, unless exact
|
||
|
|
return result;
|
||
|
|
} // decToIntegral
|