diff options
| author | Leonid S. Usov <leonid@practi.net> | 2018-10-27 00:04:35 +0300 |
|---|---|---|
| committer | William Langford <wlangfor@gmail.com> | 2019-10-22 14:11:04 -0400 |
| commit | b6be13d5de6dd7d8aad5fd871eb6b0b30fc7d7f6 (patch) | |
| tree | 8fc224d511d7ef46d0c55284d7b304657669ffec | |
| parent | 37b2d2129e5ff5d79c0f4ef08b031fa257b0bf28 (diff) | |
Add decNumber library
The library adds support for decimal numbers of arbitrary length.
Downloaded from ICU, under ICU 1.8.1 license
http://download.icu-project.org/files/decNumber/decNumber-icu-368.zip
36 files changed, 20503 insertions, 0 deletions
@@ -68,3 +68,41 @@ WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. + + +jq uses parts of the open source C library "decNumber", which is distribured +under the following license: + + +ICU License - ICU 1.8.1 and later + +COPYRIGHT AND PERMISSION NOTICE + +Copyright (c) 1995-2005 International Business Machines Corporation and others +All rights reserved. + +Permission is hereby granted, free of charge, to any person obtaining a +copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, and/or sell copies of the Software, and to permit persons +to whom the Software is furnished to do so, provided that the above +copyright notice(s) and this permission notice appear in all copies of +the Software and that both the above copyright notice(s) and this +permission notice appear in supporting documentation. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT +OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR +HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL +INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING +FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, +NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION +WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + +Except as contained in this notice, the name of a copyright holder +shall not be used in advertising or otherwise to promote the sale, use +or other dealings in this Software without prior written authorization +of the copyright holder. + diff --git a/docs/content/download/default.yml b/docs/content/download/default.yml index 8fb1216b..8ff45804 100644 --- a/docs/content/download/default.yml +++ b/docs/content/download/default.yml @@ -13,6 +13,10 @@ body: jq is licensed under the MIT license. For all of the gory details, read the file `COPYING` in the source distribution. + jq uses a C library for decimal number support. This is an ICU 1.8.1 + licensed code obtained from the ICU downloads archive + http://download.icu-project.org/files/decNumber/decNumber-icu-368.zip. + ### Linux diff --git a/src/decNumber/ICU-license.html b/src/decNumber/ICU-license.html new file mode 100644 index 00000000..bc99e8bf --- /dev/null +++ b/src/decNumber/ICU-license.html @@ -0,0 +1,45 @@ +<html>
+
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=us-ascii"></meta>
+<title>ICU License - ICU 1.8.1 and later</title>
+</head>
+
+<body>
+<h1>ICU License - ICU 1.8.1 and later</h1>
+<pre>
+COPYRIGHT AND PERMISSION NOTICE
+
+Copyright (c) 1995-2005 International Business Machines Corporation and others
+All rights reserved.
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, and/or sell copies of the Software, and to permit persons
+to whom the Software is furnished to do so, provided that the above
+copyright notice(s) and this permission notice appear in all copies of
+the Software and that both the above copyright notice(s) and this
+permission notice appear in supporting documentation.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT
+OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
+HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL
+INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING
+FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
+NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
+WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+
+Except as contained in this notice, the name of a copyright holder
+shall not be used in advertising or otherwise to promote the sale, use
+or other dealings in this Software without prior written authorization
+of the copyright holder.
+
+--------------------------------------------------------------------------------
+All trademarks and registered trademarks mentioned herein are the property of their respective owners.
+</pre>
+</body>
+</html>
diff --git a/src/decNumber/decBasic.c b/src/decNumber/decBasic.c new file mode 100644 index 00000000..1f680c4a --- /dev/null +++ b/src/decNumber/decBasic.c @@ -0,0 +1,3908 @@ +/* ------------------------------------------------------------------ */
+/* 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; |