OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libdecnumber/] [decBasic.c] - Blame information for rev 753

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 731 jeremybenn
/* Common base code for the decNumber C Library.
2
   Copyright (C) 2007, 2009 Free Software Foundation, Inc.
3
   Contributed by IBM Corporation.  Author Mike Cowlishaw.
4
 
5
   This file is part of GCC.
6
 
7
   GCC is free software; you can redistribute it and/or modify it under
8
   the terms of the GNU General Public License as published by the Free
9
   Software Foundation; either version 3, or (at your option) any later
10
   version.
11
 
12
   GCC is distributed in the hope that it will be useful, but WITHOUT ANY
13
   WARRANTY; without even the implied warranty of MERCHANTABILITY or
14
   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
15
   for more details.
16
 
17
Under Section 7 of GPL version 3, you are granted additional
18
permissions described in the GCC Runtime Library Exception, version
19
3.1, as published by the Free Software Foundation.
20
 
21
You should have received a copy of the GNU General Public License and
22
a copy of the GCC Runtime Library Exception along with this program;
23
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
24
<http://www.gnu.org/licenses/>.  */
25
 
26
/* ------------------------------------------------------------------ */
27
/* decBasic.c -- common base code for Basic decimal types             */
28
/* ------------------------------------------------------------------ */
29
/* This module comprises code that is shared between decDouble and    */
30
/* decQuad (but not decSingle).  The main arithmetic operations are   */
31
/* here (Add, Subtract, Multiply, FMA, and Division operators).       */
32
/*                                                                    */
33
/* Unlike decNumber, parameterization takes place at compile time     */
34
/* rather than at runtime.  The parameters are set in the decDouble.c */
35
/* (etc.) files, which then include this one to produce the compiled  */
36
/* code.  The functions here, therefore, are code shared between      */
37
/* multiple formats.                                                  */
38
/*                                                                    */
39
/* This must be included after decCommon.c.                           */
40
/* ------------------------------------------------------------------ */
41
/* Names here refer to decFloat rather than to decDouble, etc., and */
42
/* the functions are in strict alphabetical order. */
43
 
44
/* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */
45
/* decCommon.c */
46
#if !defined(QUAD)
47
  #error decBasic.c must be included after decCommon.c
48
#endif
49
#if SINGLE
50
  #error Routines in decBasic.c are for decDouble and decQuad only
51
#endif
52
 
53
/* Private constants */
54
#define DIVIDE      0x80000000     /* Divide operations [as flags] */
55
#define REMAINDER   0x40000000     /* .. */
56
#define DIVIDEINT   0x20000000     /* .. */
57
#define REMNEAR     0x10000000     /* .. */
58
 
59
/* Private functions (local, used only by routines in this module) */
60
static decFloat *decDivide(decFloat *, const decFloat *,
61
                              const decFloat *, decContext *, uInt);
62
static decFloat *decCanonical(decFloat *, const decFloat *);
63
static void      decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
64
                              const decFloat *);
65
static decFloat *decInfinity(decFloat *, const decFloat *);
66
static decFloat *decInvalid(decFloat *, decContext *);
67
static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
68
                              decContext *);
69
static Int       decNumCompare(const decFloat *, const decFloat *, Flag);
70
static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
71
                              enum rounding, Flag);
72
static uInt      decToInt32(const decFloat *, decContext *, enum rounding,
73
                              Flag, Flag);
74
 
75
/* ------------------------------------------------------------------ */
76
/* decCanonical -- copy a decFloat, making canonical                  */
77
/*                                                                    */
78
/*   result gets the canonicalized df                                 */
79
/*   df     is the decFloat to copy and make canonical                */
80
/*   returns result                                                   */
81
/*                                                                    */
82
/* This is exposed via decFloatCanonical for Double and Quad only.    */
83
/* This works on specials, too; no error or exception is possible.    */
84
/* ------------------------------------------------------------------ */
85
static decFloat * decCanonical(decFloat *result, const decFloat *df) {
86
  uInt encode, precode, dpd;       /* work */
87
  uInt inword, uoff, canon;        /* .. */
88
  Int  n;                          /* counter (down) */
89
  if (df!=result) *result=*df;     /* effect copy if needed */
90
  if (DFISSPECIAL(result)) {
91
    if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */
92
    /* is a NaN */
93
    DFWORD(result, 0)&=~ECONNANMASK;     /* clear ECON except selector */
94
    if (DFISCCZERO(df)) return result;  /* coefficient continuation is 0 */
95
    /* drop through to check payload */
96
    }
97
  /* return quickly if the coefficient continuation is canonical */
98
  { /* declare block */
99
  #if DOUBLE
100
    uInt sourhi=DFWORD(df, 0);
101
    uInt sourlo=DFWORD(df, 1);
102
    if (CANONDPDOFF(sourhi, 8)
103
     && CANONDPDTWO(sourhi, sourlo, 30)
104
     && CANONDPDOFF(sourlo, 20)
105
     && CANONDPDOFF(sourlo, 10)
106
     && CANONDPDOFF(sourlo, 0)) return result;
107
  #elif QUAD
108
    uInt sourhi=DFWORD(df, 0);
109
    uInt sourmh=DFWORD(df, 1);
110
    uInt sourml=DFWORD(df, 2);
111
    uInt sourlo=DFWORD(df, 3);
112
    if (CANONDPDOFF(sourhi, 4)
113
     && CANONDPDTWO(sourhi, sourmh, 26)
114
     && CANONDPDOFF(sourmh, 16)
115
     && CANONDPDOFF(sourmh, 6)
116
     && CANONDPDTWO(sourmh, sourml, 28)
117
     && CANONDPDOFF(sourml, 18)
118
     && CANONDPDOFF(sourml, 8)
119
     && CANONDPDTWO(sourml, sourlo, 30)
120
     && CANONDPDOFF(sourlo, 20)
121
     && CANONDPDOFF(sourlo, 10)
122
     && CANONDPDOFF(sourlo, 0)) return result;
123
  #endif
124
  } /* block */
125
 
126
  /* Loop to repair a non-canonical coefficent, as needed */
127
  inword=DECWORDS-1;               /* current input word */
128
  uoff=0;                           /* bit offset of declet */
129
  encode=DFWORD(result, inword);
130
  for (n=DECLETS-1; n>=0; n--) {   /* count down declets of 10 bits */
131
    dpd=encode>>uoff;
132
    uoff+=10;
133
    if (uoff>32) {                 /* crossed uInt boundary */
134
      inword--;
135
      encode=DFWORD(result, inword);
136
      uoff-=32;
137
      dpd|=encode<<(10-uoff);      /* get pending bits */
138
      }
139
    dpd&=0x3ff;                    /* clear uninteresting bits */
140
    if (dpd<0x16e) continue;       /* must be canonical */
141
    canon=BIN2DPD[DPD2BIN[dpd]];   /* determine canonical declet */
142
    if (canon==dpd) continue;      /* have canonical declet */
143
    /* need to replace declet */
144
    if (uoff>=10) {                /* all within current word */
145
      encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */
146
      encode|=canon<<(uoff-10);    /* insert the canonical form */
147
      DFWORD(result, inword)=encode;    /* .. and save */
148
      continue;
149
      }
150
    /* straddled words */
151
    precode=DFWORD(result, inword+1);   /* get previous */
152
    precode&=0xffffffff>>(10-uoff);     /* clear top bits */
153
    DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
154
    encode&=0xffffffff<<uoff;           /* clear bottom bits */
155
    encode|=canon>>(10-uoff);           /* insert canonical */
156
    DFWORD(result, inword)=encode;      /* .. and save */
157
    } /* n */
158
  return result;
159
  } /* decCanonical */
160
 
161
/* ------------------------------------------------------------------ */
162
/* decDivide -- divide operations                                     */
163
/*                                                                    */
164
/*   result gets the result of dividing dfl by dfr:                   */
165
/*   dfl    is the first decFloat (lhs)                               */
166
/*   dfr    is the second decFloat (rhs)                              */
167
/*   set    is the context                                            */
168
/*   op     is the operation selector                                 */
169
/*   returns result                                                   */
170
/*                                                                    */
171
/* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR.             */
172
/* ------------------------------------------------------------------ */
173
#define DIVCOUNT  0                /* 1 to instrument subtractions counter */
174
#define DIVBASE   ((uInt)BILLION)  /* the base used for divide */
175
#define DIVOPLEN  DECPMAX9         /* operand length ('digits' base 10**9) */
176
#define DIVACCLEN (DIVOPLEN*3)     /* accumulator length (ditto) */
177
static decFloat * decDivide(decFloat *result, const decFloat *dfl,
178
                            const decFloat *dfr, decContext *set, uInt op) {
179
  decFloat quotient;               /* for remainders */
180
  bcdnum num;                      /* for final conversion */
181
  uInt   acc[DIVACCLEN];           /* coefficent in base-billion .. */
182
  uInt   div[DIVOPLEN];            /* divisor in base-billion .. */
183
  uInt   quo[DIVOPLEN+1];          /* quotient in base-billion .. */
184
  uByte  bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */
185
  uInt   *msua, *msud, *msuq;      /* -> msu of acc, div, and quo */
186
  Int    divunits, accunits;       /* lengths */
187
  Int    quodigits;                /* digits in quotient */
188
  uInt   *lsua, *lsuq;             /* -> current acc and quo lsus */
189
  Int    length, multiplier;       /* work */
190
  uInt   carry, sign;              /* .. */
191
  uInt   *ua, *ud, *uq;            /* .. */
192
  uByte  *ub;                      /* .. */
193
  uInt   uiwork;                   /* for macros */
194
  uInt   divtop;                   /* top unit of div adjusted for estimating */
195
  #if DIVCOUNT
196
  static uInt maxcount=0;           /* worst-seen subtractions count */
197
  uInt   divcount=0;                /* subtractions count [this divide] */
198
  #endif
199
 
200
  /* calculate sign */
201
  num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
202
 
203
  if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
204
    /* NaNs are handled as usual */
205
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
206
    /* one or two infinities */
207
    if (DFISINF(dfl)) {
208
      if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */
209
      if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */
210
      /* Infinity/x is infinite and quiet, even if x=0 */
211
      DFWORD(result, 0)=num.sign;
212
      return decInfinity(result, result);
213
      }
214
    /* must be x/Infinity -- remainders are lhs */
215
    if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
216
    /* divides: return zero with correct sign and exponent depending */
217
    /* on op (Etiny for divide, 0 for divideInt) */
218
    decFloatZero(result);
219
    if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */
220
     else DFWORD(result, 0)=num.sign;         /* zeros the exponent, too */
221
    return result;
222
    }
223
  /* next, handle zero operands (x/0 and 0/x) */
224
  if (DFISZERO(dfr)) {                       /* x/0 */
225
    if (DFISZERO(dfl)) {                     /* 0/0 is undefined */
226
      decFloatZero(result);
227
      DFWORD(result, 0)=DECFLOAT_qNaN;
228
      set->status|=DEC_Division_undefined;
229
      return result;
230
      }
231
    if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */
232
    set->status|=DEC_Division_by_zero;
233
    DFWORD(result, 0)=num.sign;
234
    return decInfinity(result, result);      /* x/0 -> signed Infinity */
235
    }
236
  num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr);  /* ideal exponent */
237
  if (DFISZERO(dfl)) {                       /* 0/x (x!=0) */
238
    /* if divide, result is 0 with ideal exponent; divideInt has */
239
    /* exponent=0, remainders give zero with lower exponent */
240
    if (op&DIVIDEINT) {
241
      decFloatZero(result);
242
      DFWORD(result, 0)|=num.sign;            /* add sign */
243
      return result;
244
      }
245
    if (!(op&DIVIDE)) {                      /* a remainder */
246
      /* exponent is the minimum of the operands */
247
      num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
248
      /* if the result is zero the sign shall be sign of dfl */
249
      num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
250
      }
251
    bcdacc[0]=0;
252
    num.msd=bcdacc;                          /* -> 0 */
253
    num.lsd=bcdacc;                          /* .. */
254
    return decFinalize(result, &num, set);   /* [divide may clamp exponent] */
255
    } /* 0/x */
256
  /* [here, both operands are known to be finite and non-zero] */
257
 
258
  /* extract the operand coefficents into 'units' which are */
259
  /* base-billion; the lhs is high-aligned in acc and the msu of both */
260
  /* acc and div is at the right-hand end of array (offset length-1); */
261
  /* the quotient can need one more unit than the operands as digits */
262
  /* in it are not necessarily aligned neatly; further, the quotient */
263
  /* may not start accumulating until after the end of the initial */
264
  /* operand in acc if that is small (e.g., 1) so the accumulator */
265
  /* must have at least that number of units extra (at the ls end) */
266
  GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
267
  GETCOEFFBILL(dfr, div);
268
  /* zero the low uInts of acc */
269
  acc[0]=0;
270
  acc[1]=0;
271
  acc[2]=0;
272
  acc[3]=0;
273
  #if DOUBLE
274
    #if DIVOPLEN!=2
275
      #error Unexpected Double DIVOPLEN
276
    #endif
277
  #elif QUAD
278
  acc[4]=0;
279
  acc[5]=0;
280
  acc[6]=0;
281
  acc[7]=0;
282
    #if DIVOPLEN!=4
283
      #error Unexpected Quad DIVOPLEN
284
    #endif
285
  #endif
286
 
287
  /* set msu and lsu pointers */
288
  msua=acc+DIVACCLEN-1;       /* [leading zeros removed below] */
289
  msuq=quo+DIVOPLEN;
290
  /*[loop for div will terminate because operands are non-zero] */
291
  for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
292
  /* the initial least-significant unit of acc is set so acc appears */
293
  /* to have the same length as div. */
294
  /* This moves one position towards the least possible for each */
295
  /* iteration */
296
  divunits=(Int)(msud-div+1); /* precalculate */
297
  lsua=msua-divunits+1;       /* initial working lsu of acc */
298
  lsuq=msuq;                  /* and of quo */
299
 
300
  /* set up the estimator for the multiplier; this is the msu of div, */
301
  /* plus two bits from the unit below (if any) rounded up by one if */
302
  /* there are any non-zero bits or units below that [the extra two */
303
  /* bits makes for a much better estimate when the top unit is small] */
304
  divtop=*msud<<2;
305
  if (divunits>1) {
306
    uInt *um=msud-1;
307
    uInt d=*um;
308
    if (d>=750000000) {divtop+=3; d-=750000000;}
309
     else if (d>=500000000) {divtop+=2; d-=500000000;}
310
     else if (d>=250000000) {divtop++; d-=250000000;}
311
    if (d) divtop++;
312
     else for (um--; um>=div; um--) if (*um) {
313
      divtop++;
314
      break;
315
      }
316
    } /* >1 unit */
317
 
318
  #if DECTRACE
319
  {Int i;
320
  printf("----- div=");
321
  for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
322
  printf("\n");}
323
  #endif
324
 
325
  /* now collect up to DECPMAX+1 digits in the quotient (this may */
326
  /* need OPLEN+1 uInts if unaligned) */
327
  quodigits=0;                 /* no digits yet */
328
  for (;; lsua--) {           /* outer loop -- each input position */
329
    #if DECCHECK
330
    if (lsua<acc) {
331
      printf("Acc underrun...\n");
332
      break;
333
      }
334
    #endif
335
    #if DECTRACE
336
    printf("Outer: quodigits=%ld acc=", (LI)quodigits);
337
    for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
338
    printf("\n");
339
    #endif
340
    *lsuq=0;                   /* default unit result is 0 */
341
    for (;;) {                /* inner loop -- calculate quotient unit */
342
      /* strip leading zero units from acc (either there initially or */
343
      /* from subtraction below); this may strip all if exactly 0 */
344
      for (; *msua==0 && msua>=lsua;) msua--;
345
      accunits=(Int)(msua-lsua+1);                /* [maybe 0] */
346
      /* subtraction is only necessary and possible if there are as */
347
      /* least as many units remaining in acc for this iteration as */
348
      /* there are in div */
349
      if (accunits<divunits) {
350
        if (accunits==0) msua++;           /* restore */
351
        break;
352
        }
353
 
354
      /* If acc is longer than div then subtraction is definitely */
355
      /* possible (as msu of both is non-zero), but if they are the */
356
      /* same length a comparison is needed. */
357
      /* If a subtraction is needed then a good estimate of the */
358
      /* multiplier for the subtraction is also needed in order to */
359
      /* minimise the iterations of this inner loop because the */
360
      /* subtractions needed dominate division performance. */
361
      if (accunits==divunits) {
362
        /* compare the high divunits of acc and div: */
363
        /* acc<div:  this quotient unit is unchanged; subtraction */
364
        /*           will be possible on the next iteration */
365
        /* acc==div: quotient gains 1, set acc=0 */
366
        /* acc>div:  subtraction necessary at this position */
367
        for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
368
        /* [now at first mismatch or lsu] */
369
        if (*ud>*ua) break;                       /* next time... */
370
        if (*ud==*ua) {                           /* all compared equal */
371
          *lsuq+=1;                               /* increment result */
372
          msua=lsua;                              /* collapse acc units */
373
          *msua=0;                                 /* .. to a zero */
374
          break;
375
          }
376
 
377
        /* subtraction necessary; estimate multiplier [see above] */
378
        /* if both *msud and *msua are small it is cost-effective to */
379
        /* bring in part of the following units (if any) to get a */
380
        /* better estimate (assume some other non-zero in div) */
381
        #define DIVLO 1000000U
382
        #define DIVHI (DIVBASE/DIVLO)
383
        #if DECUSE64
384
          if (divunits>1) {
385
            /* there cannot be a *(msud-2) for DECDOUBLE so next is */
386
            /* an exact calculation unless DECQUAD (which needs to */
387
            /* assume bits out there if divunits>2) */
388
            uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
389
            uLong div=(uLong)*msud * DIVBASE + *(msud-1);
390
            #if QUAD
391
            if (divunits>2) div++;
392
            #endif
393
            mul/=div;
394
            multiplier=(Int)mul;
395
            }
396
           else multiplier=*msua/(*msud);
397
        #else
398
          if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
399
            multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
400
                      /(*msud*DIVHI + *(msud-1)/DIVLO +1);
401
            }
402
           else multiplier=(*msua<<2)/divtop;
403
        #endif
404
        }
405
       else {                                     /* accunits>divunits */
406
        /* msud is one unit 'lower' than msua, so estimate differently */
407
        #if DECUSE64
408
          uLong mul;
409
          /* as before, bring in extra digits if possible */
410
          if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
411
            mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
412
               + *(msua-2)/DIVLO;
413
            mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
414
            }
415
           else if (divunits==1) {
416
            mul=(uLong)*msua * DIVBASE + *(msua-1);
417
            mul/=*msud;       /* no more to the right */
418
            }
419
           else {
420
            mul=(uLong)(*msua) * (uInt)(DIVBASE<<2)
421
                + (*(msua-1)<<2);
422
            mul/=divtop;      /* [divtop already allows for sticky bits] */
423
            }
424
          multiplier=(Int)mul;
425
        #else
426
          multiplier=*msua * ((DIVBASE<<2)/divtop);
427
        #endif
428
        }
429
      if (multiplier==0) multiplier=1;             /* marginal case */
430
      *lsuq+=multiplier;
431
 
432
      #if DIVCOUNT
433
      /* printf("Multiplier: %ld\n", (LI)multiplier); */
434
      divcount++;
435
      #endif
436
 
437
      /* Carry out the subtraction  acc-(div*multiplier); for each */
438
      /* unit in div, do the multiply, split to units (see */
439
      /* decFloatMultiply for the algorithm), and subtract from acc */
440
      #define DIVMAGIC  2305843009U               /* 2**61/10**9 */
441
      #define DIVSHIFTA 29
442
      #define DIVSHIFTB 32
443
      carry=0;
444
      for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
445
        uInt lo, hop;
446
        #if DECUSE64
447
          uLong sub=(uLong)multiplier*(*ud)+carry;
448
          if (sub<DIVBASE) {
449
            carry=0;
450
            lo=(uInt)sub;
451
            }
452
           else {
453
            hop=(uInt)(sub>>DIVSHIFTA);
454
            carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
455
            /* the estimate is now in hi; now calculate sub-hi*10**9 */
456
            /* to get the remainder (which will be <DIVBASE)) */
457
            lo=(uInt)sub;
458
            lo-=carry*DIVBASE;                    /* low word of result */
459
            if (lo>=DIVBASE) {
460
              lo-=DIVBASE;                        /* correct by +1 */
461
              carry++;
462
              }
463
            }
464
        #else /* 32-bit */
465
          uInt hi;
466
          /* calculate multiplier*(*ud) into hi and lo */
467
          LONGMUL32HI(hi, *ud, multiplier);       /* get the high word */
468
          lo=multiplier*(*ud);                    /* .. and the low */
469
          lo+=carry;                              /* add the old hi */
470
          carry=hi+(lo<carry);                    /* .. with any carry */
471
          if (carry || lo>=DIVBASE) {             /* split is needed */
472
            hop=(carry<<3)+(lo>>DIVSHIFTA);       /* hi:lo/2**29 */
473
            LONGMUL32HI(carry, hop, DIVMAGIC);    /* only need the high word */
474
            /* [DIVSHIFTB is 32, so carry can be used directly] */
475
            /* the estimate is now in carry; now calculate hi:lo-est*10**9; */
476
            /* happily the top word of the result is irrelevant because it */
477
            /* will always be zero so this needs only one multiplication */
478
            lo-=(carry*DIVBASE);
479
            /* the correction here will be at most +1; do it */
480
            if (lo>=DIVBASE) {
481
              lo-=DIVBASE;
482
              carry++;
483
              }
484
            }
485
        #endif
486
        if (lo>*ua) {              /* borrow needed */
487
          *ua+=DIVBASE;
488
          carry++;
489
          }
490
        *ua-=lo;
491
        } /* ud loop */
492
      if (carry) *ua-=carry;       /* accdigits>divdigits [cannot borrow] */
493
      } /* inner loop */
494
 
495
    /* the outer loop terminates when there is either an exact result */
496
    /* or enough digits; first update the quotient digit count and */
497
    /* pointer (if any significant digits) */
498
    #if DECTRACE
499
    if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
500
    #endif
501
    if (quodigits) {
502
      quodigits+=9;                /* had leading unit earlier */
503
      lsuq--;
504
      if (quodigits>DECPMAX+1) break;   /* have enough */
505
      }
506
     else if (*lsuq) {             /* first quotient digits */
507
      const uInt *pow;
508
      for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
509
      lsuq--;
510
      /* [cannot have >DECPMAX+1 on first unit] */
511
      }
512
 
513
    if (*msua!=0) continue;         /* not an exact result */
514
    /* acc is zero iff used all of original units and zero down to lsua */
515
    /* (must also continue to original lsu for correct quotient length) */
516
    if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
517
    for (; msua>lsua && *msua==0;) msua--;
518
    if (*msua==0 && msua==lsua) break;
519
    } /* outer loop */
520
 
521
  /* all of the original operand in acc has been covered at this point */
522
  /* quotient now has at least DECPMAX+2 digits */
523
  /* *msua is now non-0 if inexact and sticky bits */
524
  /* lsuq is one below the last uint of the quotient */
525
  lsuq++;                          /* set -> true lsu of quo */
526
  if (*msua) *lsuq|=1;             /* apply sticky bit */
527
 
528
  /* quo now holds the (unrounded) quotient in base-billion; one */
529
  /* base-billion 'digit' per uInt. */
530
  #if DECTRACE
531
  printf("DivQuo:");
532
  for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
533
  printf("\n");
534
  #endif
535
 
536
  /* Now convert to BCD for rounding and cleanup, starting from the */
537
  /* most significant end [offset by one into bcdacc to leave room */
538
  /* for a possible carry digit if rounding for REMNEAR is needed] */
539
  for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
540
    uInt top, mid, rem;                 /* work */
541
    if (*uq==0) {                        /* no split needed */
542
      UBFROMUI(ub, 0);                   /* clear 9 BCD8s */
543
      UBFROMUI(ub+4, 0);         /* .. */
544
      *(ub+8)=0;                 /* .. */
545
      continue;
546
      }
547
    /* *uq is non-zero -- split the base-billion digit into */
548
    /* hi, mid, and low three-digits */
549
    #define divsplit9 1000000           /* divisor */
550
    #define divsplit6 1000              /* divisor */
551
    /* The splitting is done by simple divides and remainders, */
552
    /* assuming the compiler will optimize these [GCC does] */
553
    top=*uq/divsplit9;
554
    rem=*uq%divsplit9;
555
    mid=rem/divsplit6;
556
    rem=rem%divsplit6;
557
    /* lay out the nine BCD digits (plus one unwanted byte) */
558
    UBFROMUI(ub,   UBTOUI(&BIN2BCD8[top*4]));
559
    UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
560
    UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
561
    } /* BCD conversion loop */
562
  ub--;                                 /* -> lsu */
563
 
564
  /* complete the bcdnum; quodigits is correct, so the position of */
565
  /* the first non-zero is known */
566
  num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
567
  num.lsd=ub;
568
 
569
  /* make exponent adjustments, etc */
570
  if (lsua<acc+DIVACCLEN-DIVOPLEN) {    /* used extra digits */
571
    num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
572
    /* if the result was exact then there may be up to 8 extra */
573
    /* trailing zeros in the overflowed quotient final unit */
574
    if (*msua==0) {
575
      for (; *ub==0;) ub--;              /* drop zeros */
576
      num.exponent+=(Int)(num.lsd-ub);  /* and adjust exponent */
577
      num.lsd=ub;
578
      }
579
    } /* adjustment needed */
580
 
581
  #if DIVCOUNT
582
  if (divcount>maxcount) {              /* new high-water nark */
583
    maxcount=divcount;
584
    printf("DivNewMaxCount: %ld\n", (LI)maxcount);
585
    }
586
  #endif
587
 
588
  if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */
589
 
590
  /* Is DIVIDEINT or a remainder; there is more to do -- first form */
591
  /* the integer (this is done 'after the fact', unlike as in */
592
  /* decNumber, so as not to tax DIVIDE) */
593
 
594
  /* The first non-zero digit will be in the first 9 digits, known */
595
  /* from quodigits and num.msd, so there is always space for DECPMAX */
596
  /* digits */
597
 
598
  length=(Int)(num.lsd-num.msd+1);
599
  /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */
600
 
601
  if (length+num.exponent>DECPMAX) { /* cannot fit */
602
    decFloatZero(result);
603
    DFWORD(result, 0)=DECFLOAT_qNaN;
604
    set->status|=DEC_Division_impossible;
605
    return result;
606
    }
607
 
608
  if (num.exponent>=0) {    /* already an int, or need pad zeros */
609
    for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
610
    num.lsd+=num.exponent;
611
    }
612
   else {                          /* too long: round or truncate needed */
613
    Int drop=-num.exponent;
614
    if (!(op&REMNEAR)) {           /* simple truncate */
615
      num.lsd-=drop;
616
      if (num.lsd<num.msd) {       /* truncated all */
617
        num.lsd=num.msd;           /* make 0 */
618
        *num.lsd=0;                 /* .. [sign still relevant] */
619
        }
620
      }
621
     else {                        /* round to nearest even [sigh] */
622
      /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */
623
      /* (this is a special case of Quantize -- q.v. for commentary) */
624
      uByte *roundat;              /* -> re-round digit */
625
      uByte reround;               /* reround value */
626
      *(num.msd-1)=0;               /* in case of left carry, or make 0 */
627
      if (drop<length) roundat=num.lsd-drop+1;
628
       else if (drop==length) roundat=num.msd;
629
       else roundat=num.msd-1;     /* [-> 0] */
630
      reround=*roundat;
631
      for (ub=roundat+1; ub<=num.lsd; ub++) {
632
        if (*ub!=0) {
633
          reround=DECSTICKYTAB[reround];
634
          break;
635
          }
636
        } /* check stickies */
637
      if (roundat>num.msd) num.lsd=roundat-1;
638
       else {
639
        num.msd--;                           /* use the 0 .. */
640
        num.lsd=num.msd;                     /* .. at the new MSD place */
641
        }
642
      if (reround!=0) {               /* discarding non-zero */
643
        uInt bump=0;
644
        /* rounding is DEC_ROUND_HALF_EVEN always */
645
        if (reround>5) bump=1;               /* >0.5 goes up */
646
         else if (reround==5)                /* exactly 0.5000 .. */
647
          bump=*(num.lsd) & 0x01;            /* .. up iff [new] lsd is odd */
648
        if (bump!=0) {                        /* need increment */
649
          /* increment the coefficient; this might end up with 1000... */
650
          ub=num.lsd;
651
          for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
652
          for (; *ub==9; ub--) *ub=0;         /* at most 3 more */
653
          *ub+=1;
654
          if (ub<num.msd) num.msd--;         /* carried */
655
          } /* bump needed */
656
        } /* reround!=0 */
657
      } /* remnear */
658
    } /* round or truncate needed */
659
  num.exponent=0;                             /* all paths */
660
  /*decShowNum(&num, "int"); */
661
 
662
  if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */
663
 
664
  /* Have a remainder to calculate */
665
  decFinalize(&quotient, &num, set);         /* lay out the integer so far */
666
  DFWORD(&quotient, 0)^=DECFLOAT_Sign;        /* negate it */
667
  sign=DFWORD(dfl, 0);                        /* save sign of dfl */
668
  decFloatFMA(result, &quotient, dfr, dfl, set);
669
  if (!DFISZERO(result)) return result;
670
  /* if the result is zero the sign shall be sign of dfl */
671
  DFWORD(&quotient, 0)=sign;                  /* construct decFloat of sign */
672
  return decFloatCopySign(result, result, &quotient);
673
  } /* decDivide */
674
 
675
/* ------------------------------------------------------------------ */
676
/* decFiniteMultiply -- multiply two finite decFloats                 */
677
/*                                                                    */
678
/*   num    gets the result of multiplying dfl and dfr                */
679
/*   bcdacc .. with the coefficient in this array                     */
680
/*   dfl    is the first decFloat (lhs)                               */
681
/*   dfr    is the second decFloat (rhs)                              */
682
/*                                                                    */
683
/* This effects the multiplication of two decFloats, both known to be */
684
/* finite, leaving the result in a bcdnum ready for decFinalize (for  */
685
/* use in Multiply) or in a following addition (FMA).                 */
686
/*                                                                    */
687
/* bcdacc must have space for at least DECPMAX9*18+1 bytes.           */
688
/* No error is possible and no status is set.                         */
689
/* ------------------------------------------------------------------ */
690
/* This routine has two separate implementations of the core */
691
/* multiplication; both using base-billion.  One uses only 32-bit */
692
/* variables (Ints and uInts) or smaller; the other uses uLongs (for */
693
/* multiplication and addition only).  Both implementations cover */
694
/* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
695
/* comparisons.  In any one compilation only one implementation for */
696
/* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
697
/* version is forced. */
698
/* */
699
/* Historical note: an earlier version of this code also supported the */
700
/* 256-bit format and has been preserved.  That is somewhat trickier */
701
/* during lazy carry splitting because the initial quotient estimate */
702
/* (est) can exceed 32 bits. */
703
 
704
#define MULTBASE  ((uInt)BILLION)  /* the base used for multiply */
705
#define MULOPLEN  DECPMAX9         /* operand length ('digits' base 10**9) */
706
#define MULACCLEN (MULOPLEN*2)              /* accumulator length (ditto) */
707
#define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */
708
 
709
/* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */
710
#if DECEMAXD>9
711
  #error Exponent may overflow when doubled for Multiply
712
#endif
713
#if MULACCLEN!=(MULACCLEN/4)*4
714
  /* This assumption is used below only for initialization */
715
  #error MULACCLEN is not a multiple of 4
716
#endif
717
 
718
static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
719
                              const decFloat *dfl, const decFloat *dfr) {
720
  uInt   bufl[MULOPLEN];           /* left  coefficient (base-billion) */
721
  uInt   bufr[MULOPLEN];           /* right coefficient (base-billion) */
722
  uInt   *ui, *uj;                 /* work */
723
  uByte  *ub;                      /* .. */
724
  uInt   uiwork;                   /* for macros */
725
 
726
  #if DECUSE64
727
  uLong  accl[MULACCLEN];          /* lazy accumulator (base-billion+) */
728
  uLong  *pl;                      /* work -> lazy accumulator */
729
  uInt   acc[MULACCLEN];           /* coefficent in base-billion .. */
730
  #else
731
  uInt   acc[MULACCLEN*2];         /* accumulator in base-billion .. */
732
  #endif
733
  uInt   *pa;                      /* work -> accumulator */
734
  /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */
735
 
736
  /* Calculate sign and exponent */
737
  num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
738
  num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */
739
 
740
  /* Extract the coefficients and prepare the accumulator */
741
  /* the coefficients of the operands are decoded into base-billion */
742
  /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */
743
  /* appropriate size. */
744
  GETCOEFFBILL(dfl, bufl);
745
  GETCOEFFBILL(dfr, bufr);
746
  #if DECTRACE && 0
747
    printf("CoeffbL:");
748
    for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
749
    printf("\n");
750
    printf("CoeffbR:");
751
    for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
752
    printf("\n");
753
  #endif
754
 
755
  /* start the 64-bit/32-bit differing paths... */
756
#if DECUSE64
757
 
758
  /* zero the accumulator */
759
  #if MULACCLEN==4
760
    accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
761
  #else                                      /* use a loop */
762
    /* MULACCLEN is a multiple of four, asserted above */
763
    for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
764
      *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */
765
      } /* pl */
766
  #endif
767
 
768
  /* Effect the multiplication */
769
  /* The multiplcation proceeds using MFC's lazy-carry resolution */
770
  /* algorithm from decNumber.  First, the multiplication is */
771
  /* effected, allowing accumulation of the partial products (which */
772
  /* are in base-billion at each column position) into 64 bits */
773
  /* without resolving back to base=billion after each addition. */
774
  /* These 64-bit numbers (which may contain up to 19 decimal digits) */
775
  /* are then split using the Clark & Cowlishaw algorithm (see below). */
776
  /* [Testing for 0 in the inner loop is not really a 'win'] */
777
  for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */
778
    if (*ui==0) continue;                  /* product cannot affect result */
779
    pl=accl+(ui-bufr);                    /* where to add the lhs */
780
    for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */
781
      /* if (*uj==0) continue;            // product cannot affect result */
782
      *pl+=((uLong)*ui)*(*uj);
783
      } /* uj */
784
    } /* ui */
785
 
786
  /* The 64-bit carries must now be resolved; this means that a */
787
  /* quotient/remainder has to be calculated for base-billion (1E+9). */
788
  /* For this, Clark & Cowlishaw's quotient estimation approach (also */
789
  /* used in decNumber) is needed, because 64-bit divide is generally */
790
  /* extremely slow on 32-bit machines, and may be slower than this */
791
  /* approach even on 64-bit machines.  This algorithm splits X */
792
  /* using: */
793
  /* */
794
  /*   magic=2**(A+B)/1E+9;   // 'magic number' */
795
  /*   hop=X/2**A;            // high order part of X (by shift) */
796
  /*   est=magic*hop/2**B     // quotient estimate (may be low by 1) */
797
  /* */
798
  /* A and B are quite constrained; hop and magic must fit in 32 bits, */
799
  /* and 2**(A+B) must be as large as possible (which is 2**61 if */
800
  /* magic is to fit).  Further, maxX increases with the length of */
801
  /* the operands (and hence the number of partial products */
802
  /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
803
  /* */
804
  /* It can be shown that when OPLEN is 2 then the maximum error in */
805
  /* the estimated quotient is <1, but for larger maximum x the */
806
  /* maximum error is above 1 so a correction that is >1 may be */
807
  /* needed.  Values of A and B are chosen to satisfy the constraints */
808
  /* just mentioned while minimizing the maximum error (and hence the */
809
  /* maximum correction), as shown in the following table: */
810
  /* */
811
  /*   Type    OPLEN   A   B     maxX    maxError  maxCorrection */
812
  /*   --------------------------------------------------------- */
813
  /*   DOUBLE    2    29  32  <2*10**18    0.63       1 */
814
  /*   QUAD      4    30  31  <4*10**18    1.17       2 */
815
  /* */
816
  /* In the OPLEN==2 case there is most choice, but the value for B */
817
  /* of 32 has a big advantage as then the calculation of the */
818
  /* estimate requires no shifting; the compiler can extract the high */
819
  /* word directly after multiplying magic*hop. */
820
  #define MULMAGIC 2305843009U          /* 2**61/10**9  [both cases] */
821
  #if DOUBLE
822
    #define MULSHIFTA 29
823
    #define MULSHIFTB 32
824
  #elif QUAD
825
    #define MULSHIFTA 30
826
    #define MULSHIFTB 31
827
  #else
828
    #error Unexpected type
829
  #endif
830
 
831
  #if DECTRACE
832
  printf("MulAccl:");
833
  for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
834
    printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
835
  printf("\n");
836
  #endif
837
 
838
  for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */
839
    uInt lo, hop;                       /* work */
840
    uInt est;                           /* cannot exceed 4E+9 */
841
    if (*pl>=MULTBASE) {
842
      /* *pl holds a binary number which needs to be split */
843
      hop=(uInt)(*pl>>MULSHIFTA);
844
      est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
845
      /* the estimate is now in est; now calculate hi:lo-est*10**9; */
846
      /* happily the top word of the result is irrelevant because it */
847
      /* will always be zero so this needs only one multiplication */
848
      lo=(uInt)(*pl-((uLong)est*MULTBASE));  /* low word of result */
849
      /* If QUAD, the correction here could be +2 */
850
      if (lo>=MULTBASE) {
851
        lo-=MULTBASE;                   /* correct by +1 */
852
        est++;
853
        #if QUAD
854
        /* may need to correct by +2 */
855
        if (lo>=MULTBASE) {
856
          lo-=MULTBASE;
857
          est++;
858
          }
859
        #endif
860
        }
861
      /* finally place lo as the new coefficient 'digit' and add est to */
862
      /* the next place up [this is safe because this path is never */
863
      /* taken on the final iteration as *pl will fit] */
864
      *pa=lo;
865
      *(pl+1)+=est;
866
      } /* *pl needed split */
867
     else {                             /* *pl<MULTBASE */
868
      *pa=(uInt)*pl;                    /* just copy across */
869
      }
870
    } /* pl loop */
871
 
872
#else  /* 32-bit */
873
  for (pa=acc;; pa+=4) {                     /* zero the accumulator */
874
    *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0;  /* [reduce overhead] */
875
    if (pa==acc+MULACCLEN*2-4) break;        /* multiple of 4 asserted */
876
    } /* pa */
877
 
878
  /* Effect the multiplication */
879
  /* uLongs are not available (and in particular, there is no uLong */
880
  /* divide) but it is still possible to use MFC's lazy-carry */
881
  /* resolution algorithm from decNumber.  First, the multiplication */
882
  /* is effected, allowing accumulation of the partial products */
883
  /* (which are in base-billion at each column position) into 64 bits */
884
  /* [with the high-order 32 bits in each position being held at */
885
  /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */
886
  /* These 64-bit numbers (which may contain up to 19 decimal digits) */
887
  /* are then split using the Clark & Cowlishaw algorithm (see */
888
  /* below). */
889
  for (ui=bufr;; ui++) {                /* over each item in rhs */
890
    uInt hi, lo;                        /* words of exact multiply result */
891
    pa=acc+(ui-bufr);                   /* where to add the lhs */
892
    for (uj=bufl;; uj++, pa++) {        /* over each item in lhs */
893
      LONGMUL32HI(hi, *ui, *uj);        /* calculate product of digits */
894
      lo=(*ui)*(*uj);                   /* .. */
895
      *pa+=lo;                          /* accumulate low bits and .. */
896
      *(pa+MULACCLEN)+=hi+(*pa<lo);     /* .. high bits with any carry */
897
      if (uj==bufl+MULOPLEN-1) break;
898
      }
899
    if (ui==bufr+MULOPLEN-1) break;
900
    }
901
 
902
  /* The 64-bit carries must now be resolved; this means that a */
903
  /* quotient/remainder has to be calculated for base-billion (1E+9). */
904
  /* For this, Clark & Cowlishaw's quotient estimation approach (also */
905
  /* used in decNumber) is needed, because 64-bit divide is generally */
906
  /* extremely slow on 32-bit machines.  This algorithm splits X */
907
  /* using: */
908
  /* */
909
  /*   magic=2**(A+B)/1E+9;   // 'magic number' */
910
  /*   hop=X/2**A;            // high order part of X (by shift) */
911
  /*   est=magic*hop/2**B     // quotient estimate (may be low by 1) */
912
  /* */
913
  /* A and B are quite constrained; hop and magic must fit in 32 bits, */
914
  /* and 2**(A+B) must be as large as possible (which is 2**61 if */
915
  /* magic is to fit).  Further, maxX increases with the length of */
916
  /* the operands (and hence the number of partial products */
917
  /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
918
  /* */
919
  /* It can be shown that when OPLEN is 2 then the maximum error in */
920
  /* the estimated quotient is <1, but for larger maximum x the */
921
  /* maximum error is above 1 so a correction that is >1 may be */
922
  /* needed.  Values of A and B are chosen to satisfy the constraints */
923
  /* just mentioned while minimizing the maximum error (and hence the */
924
  /* maximum correction), as shown in the following table: */
925
  /* */
926
  /*   Type    OPLEN   A   B     maxX    maxError  maxCorrection */
927
  /*   --------------------------------------------------------- */
928
  /*   DOUBLE    2    29  32  <2*10**18    0.63       1 */
929
  /*   QUAD      4    30  31  <4*10**18    1.17       2 */
930
  /* */
931
  /* In the OPLEN==2 case there is most choice, but the value for B */
932
  /* of 32 has a big advantage as then the calculation of the */
933
  /* estimate requires no shifting; the high word is simply */
934
  /* calculated from multiplying magic*hop. */
935
  #define MULMAGIC 2305843009U          /* 2**61/10**9  [both cases] */
936
  #if DOUBLE
937
    #define MULSHIFTA 29
938
    #define MULSHIFTB 32
939
  #elif QUAD
940
    #define MULSHIFTA 30
941
    #define MULSHIFTB 31
942
  #else
943
    #error Unexpected type
944
  #endif
945
 
946
  #if DECTRACE
947
  printf("MulHiLo:");
948
  for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
949
    printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
950
  printf("\n");
951
  #endif
952
 
953
  for (pa=acc;; pa++) {                 /* each low uInt */
954
    uInt hi, lo;                        /* words of exact multiply result */
955
    uInt hop, estlo;                    /* work */
956
    #if QUAD
957
    uInt esthi;                         /* .. */
958
    #endif
959
 
960
    lo=*pa;
961
    hi=*(pa+MULACCLEN);                 /* top 32 bits */
962
    /* hi and lo now hold a binary number which needs to be split */
963
 
964
    #if DOUBLE
965
      hop=(hi<<3)+(lo>>MULSHIFTA);      /* hi:lo/2**29 */
966
      LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */
967
      /* [MULSHIFTB is 32, so estlo can be used directly] */
968
      /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */
969
      /* happily the top word of the result is irrelevant because it */
970
      /* will always be zero so this needs only one multiplication */
971
      lo-=(estlo*MULTBASE);
972
      /* esthi=0;                       // high word is ignored below */
973
      /* the correction here will be at most +1; do it */
974
      if (lo>=MULTBASE) {
975
        lo-=MULTBASE;
976
        estlo++;
977
        }
978
    #elif QUAD
979
      hop=(hi<<2)+(lo>>MULSHIFTA);      /* hi:lo/2**30 */
980
      LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */
981
      estlo=hop*MULMAGIC;               /* .. so low word needed */
982
      estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */
983
      /* esthi=0;                       // high word is ignored below */
984
      lo-=(estlo*MULTBASE);             /* as above */
985
      /* the correction here could be +1 or +2 */
986
      if (lo>=MULTBASE) {
987
        lo-=MULTBASE;
988
        estlo++;
989
        }
990
      if (lo>=MULTBASE) {
991
        lo-=MULTBASE;
992
        estlo++;
993
        }
994
    #else
995
      #error Unexpected type
996
    #endif
997
 
998
    /* finally place lo as the new accumulator digit and add est to */
999
    /* the next place up; this latter add could cause a carry of 1 */
1000
    /* to the high word of the next place */
1001
    *pa=lo;
1002
    *(pa+1)+=estlo;
1003
    /* esthi is always 0 for DOUBLE and QUAD so this is skipped */
1004
    /* *(pa+1+MULACCLEN)+=esthi; */
1005
    if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */
1006
    if (pa==acc+MULACCLEN-2) break;          /* [MULACCLEN-1 will never need split] */
1007
    } /* pa loop */
1008
#endif
1009
 
1010
  /* At this point, whether using the 64-bit or the 32-bit paths, the */
1011
  /* accumulator now holds the (unrounded) result in base-billion; */
1012
  /* one base-billion 'digit' per uInt. */
1013
  #if DECTRACE
1014
  printf("MultAcc:");
1015
  for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
1016
  printf("\n");
1017
  #endif
1018
 
1019
  /* Now convert to BCD for rounding and cleanup, starting from the */
1020
  /* most significant end */
1021
  pa=acc+MULACCLEN-1;
1022
  if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */
1023
   else {                               /* >=1 word of leading zeros */
1024
    num->msd=bcdacc;                    /* known leading zeros are gone */
1025
    pa--;                               /* skip first word .. */
1026
    for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */
1027
    }
1028
  for (ub=bcdacc;; pa--, ub+=9) {
1029
    if (*pa!=0) {                        /* split(s) needed */
1030
      uInt top, mid, rem;               /* work */
1031
      /* *pa is non-zero -- split the base-billion acc digit into */
1032
      /* hi, mid, and low three-digits */
1033
      #define mulsplit9 1000000         /* divisor */
1034
      #define mulsplit6 1000            /* divisor */
1035
      /* The splitting is done by simple divides and remainders, */
1036
      /* assuming the compiler will optimize these where useful */
1037
      /* [GCC does] */
1038
      top=*pa/mulsplit9;
1039
      rem=*pa%mulsplit9;
1040
      mid=rem/mulsplit6;
1041
      rem=rem%mulsplit6;
1042
      /* lay out the nine BCD digits (plus one unwanted byte) */
1043
      UBFROMUI(ub,   UBTOUI(&BIN2BCD8[top*4]));
1044
      UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
1045
      UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
1046
      }
1047
     else {                             /* *pa==0 */
1048
      UBFROMUI(ub, 0);                   /* clear 9 BCD8s */
1049
      UBFROMUI(ub+4, 0);         /* .. */
1050
      *(ub+8)=0;                 /* .. */
1051
      }
1052
    if (pa==acc) break;
1053
    } /* BCD conversion loop */
1054
 
1055
  num->lsd=ub+8;                        /* complete the bcdnum .. */
1056
 
1057
  #if DECTRACE
1058
  decShowNum(num, "postmult");
1059
  decFloatShow(dfl, "dfl");
1060
  decFloatShow(dfr, "dfr");
1061
  #endif
1062
  return;
1063
  } /* decFiniteMultiply */
1064
 
1065
/* ------------------------------------------------------------------ */
1066
/* decFloatAbs -- absolute value, heeding NaNs, etc.                  */
1067
/*                                                                    */
1068
/*   result gets the canonicalized df with sign 0                     */
1069
/*   df     is the decFloat to abs                                    */
1070
/*   set    is the context                                            */
1071
/*   returns result                                                   */
1072
/*                                                                    */
1073
/* This has the same effect as decFloatPlus unless df is negative,    */
1074
/* in which case it has the same effect as decFloatMinus.  The        */
1075
/* effect is also the same as decFloatCopyAbs except that NaNs are    */
1076
/* handled normally (the sign of a NaN is not affected, and an sNaN   */
1077
/* will signal) and the result will be canonical.                     */
1078
/* ------------------------------------------------------------------ */
1079
decFloat * decFloatAbs(decFloat *result, const decFloat *df,
1080
                       decContext *set) {
1081
  if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
1082
  decCanonical(result, df);             /* copy and check */
1083
  DFBYTE(result, 0)&=~0x80;              /* zero sign bit */
1084
  return result;
1085
  } /* decFloatAbs */
1086
 
1087
/* ------------------------------------------------------------------ */
1088
/* decFloatAdd -- add two decFloats                                   */
1089
/*                                                                    */
1090
/*   result gets the result of adding dfl and dfr:                    */
1091
/*   dfl    is the first decFloat (lhs)                               */
1092
/*   dfr    is the second decFloat (rhs)                              */
1093
/*   set    is the context                                            */
1094
/*   returns result                                                   */
1095
/*                                                                    */
1096
/* ------------------------------------------------------------------ */
1097
#if QUAD
1098
/* Table for testing MSDs for fastpath elimination; returns the MSD of */
1099
/* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */
1100
/* Infinities return -32 and NaNs return -128 so that summing the two */
1101
/* MSDs also allows rapid tests for the Specials (see code below). */
1102
const Int DECTESTMSD[64]={
1103
  0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5,   6,    7,
1104
  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128,
1105
  0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5,   6,    7,
1106
  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128};
1107
#else
1108
/* The table for testing MSDs is shared between the modules */
1109
extern const Int DECTESTMSD[64];
1110
#endif
1111
 
1112
decFloat * decFloatAdd(decFloat *result,
1113
                       const decFloat *dfl, const decFloat *dfr,
1114
                       decContext *set) {
1115
  bcdnum num;                      /* for final conversion */
1116
  Int    bexpl, bexpr;             /* left and right biased exponents */
1117
  uByte  *ub, *us, *ut;            /* work */
1118
  uInt   uiwork;                   /* for macros */
1119
  #if QUAD
1120
  uShort uswork;                   /* .. */
1121
  #endif
1122
 
1123
  uInt sourhil, sourhir;           /* top words from source decFloats */
1124
                                   /* [valid only through end of */
1125
                                   /* fastpath code -- before swap] */
1126
  uInt diffsign;                   /* non-zero if signs differ */
1127
  uInt carry;                      /* carry: 0 or 1 before add loop */
1128
  Int  overlap;                    /* coefficient overlap (if full) */
1129
  Int  summ;                       /* sum of the MSDs */
1130
  /* the following buffers hold coefficients with various alignments */
1131
  /* (see commentary and diagrams below) */
1132
  uByte acc[4+2+DECPMAX*3+8];
1133
  uByte buf[4+2+DECPMAX*2];
1134
  uByte *umsd, *ulsd;              /* local MSD and LSD pointers */
1135
 
1136
  #if DECLITEND
1137
    #define CARRYPAT 0x01000000    /* carry=1 pattern */
1138
  #else
1139
    #define CARRYPAT 0x00000001    /* carry=1 pattern */
1140
  #endif
1141
 
1142
  /* Start decoding the arguments */
1143
  /* The initial exponents are placed into the opposite Ints to */
1144
  /* that which might be expected; there are two sets of data to */
1145
  /* keep track of (each decFloat and the corresponding exponent), */
1146
  /* and this scheme means that at the swap point (after comparing */
1147
  /* exponents) only one pair of words needs to be swapped */
1148
  /* whichever path is taken (thereby minimising worst-case path). */
1149
  /* The calculated exponents will be nonsense when the arguments are */
1150
  /* Special, but are not used in that path */
1151
  sourhil=DFWORD(dfl, 0);           /* LHS top word */
1152
  summ=DECTESTMSD[sourhil>>26];    /* get first MSD for testing */
1153
  bexpr=DECCOMBEXP[sourhil>>26];   /* get exponent high bits (in place) */
1154
  bexpr+=GETECON(dfl);             /* .. + continuation */
1155
 
1156
  sourhir=DFWORD(dfr, 0);           /* RHS top word */
1157
  summ+=DECTESTMSD[sourhir>>26];   /* sum MSDs for testing */
1158
  bexpl=DECCOMBEXP[sourhir>>26];
1159
  bexpl+=GETECON(dfr);
1160
 
1161
  /* here bexpr has biased exponent from lhs, and vice versa */
1162
 
1163
  diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
1164
 
1165
  /* now determine whether to take a fast path or the full-function */
1166
  /* slow path.  The slow path must be taken when: */
1167
  /*   -- both numbers are finite, and: */
1168
  /*         the exponents are different, or */
1169
  /*         the signs are different, or */
1170
  /*         the sum of the MSDs is >8 (hence might overflow) */
1171
  /* specialness and the sum of the MSDs can be tested at once using */
1172
  /* the summ value just calculated, so the test for specials is no */
1173
  /* longer on the worst-case path (as of 3.60) */
1174
 
1175
  if (summ<=8) {                   /* MSD+MSD is good, or there is a special */
1176
    if (summ<0) {                   /* there is a special */
1177
      /* Inf+Inf would give -64; Inf+finite is -32 or higher */
1178
      if (summ<-64) return decNaNs(result, dfl, dfr, set);  /* one or two NaNs */
1179
      /* two infinities with different signs is invalid */
1180
      if (summ==-64 && diffsign) return decInvalid(result, set);
1181
      if (DFISINF(dfl)) return decInfinity(result, dfl);    /* LHS is infinite */
1182
      return decInfinity(result, dfr);                      /* RHS must be Inf */
1183
      }
1184
    /* Here when both arguments are finite; fast path is possible */
1185
    /* (currently only for aligned and same-sign) */
1186
    if (bexpr==bexpl && !diffsign) {
1187
      uInt tac[DECLETS+1];              /* base-1000 coefficient */
1188
      uInt encode;                      /* work */
1189
 
1190
      /* Get one coefficient as base-1000 and add the other */
1191
      GETCOEFFTHOU(dfl, tac);           /* least-significant goes to [0] */
1192
      ADDCOEFFTHOU(dfr, tac);
1193
      /* here the sum of the MSDs (plus any carry) will be <10 due to */
1194
      /* the fastpath test earlier */
1195
 
1196
      /* construct the result; low word is the same for both formats */
1197
      encode =BIN2DPD[tac[0]];
1198
      encode|=BIN2DPD[tac[1]]<<10;
1199
      encode|=BIN2DPD[tac[2]]<<20;
1200
      encode|=BIN2DPD[tac[3]]<<30;
1201
      DFWORD(result, (DECBYTES/4)-1)=encode;
1202
 
1203
      /* collect next two declets (all that remains, for Double) */
1204
      encode =BIN2DPD[tac[3]]>>2;
1205
      encode|=BIN2DPD[tac[4]]<<8;
1206
 
1207
      #if QUAD
1208
      /* complete and lay out middling words */
1209
      encode|=BIN2DPD[tac[5]]<<18;
1210
      encode|=BIN2DPD[tac[6]]<<28;
1211
      DFWORD(result, 2)=encode;
1212
 
1213
      encode =BIN2DPD[tac[6]]>>4;
1214
      encode|=BIN2DPD[tac[7]]<<6;
1215
      encode|=BIN2DPD[tac[8]]<<16;
1216
      encode|=BIN2DPD[tac[9]]<<26;
1217
      DFWORD(result, 1)=encode;
1218
 
1219
      /* and final two declets */
1220
      encode =BIN2DPD[tac[9]]>>6;
1221
      encode|=BIN2DPD[tac[10]]<<4;
1222
      #endif
1223
 
1224
      /* add exponent continuation and sign (from either argument) */
1225
      encode|=sourhil & (ECONMASK | DECFLOAT_Sign);
1226
 
1227
      /* create lookup index = MSD + top two bits of biased exponent <<4 */
1228
      tac[DECLETS]|=(bexpl>>DECECONL)<<4;
1229
      encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */
1230
      DFWORD(result, 0)=encode;   /* complete */
1231
 
1232
      /* decFloatShow(result, ">"); */
1233
      return result;
1234
      } /* fast path OK */
1235
    /* drop through to slow path */
1236
    } /* low sum or Special(s) */
1237
 
1238
  /* Slow path required -- arguments are finite and might overflow,   */
1239
  /* or require alignment, or might have different signs              */
1240
 
1241
  /* now swap either exponents or argument pointers */
1242
  if (bexpl<=bexpr) {
1243
    /* original left is bigger */
1244
    Int bexpswap=bexpl;
1245
    bexpl=bexpr;
1246
    bexpr=bexpswap;
1247
    /* printf("left bigger\n"); */
1248
    }
1249
   else {
1250
    const decFloat *dfswap=dfl;
1251
    dfl=dfr;
1252
    dfr=dfswap;
1253
    /* printf("right bigger\n"); */
1254
    }
1255
  /* [here dfl and bexpl refer to the datum with the larger exponent, */
1256
  /* of if the exponents are equal then the original LHS argument] */
1257
 
1258
  /* if lhs is zero then result will be the rhs (now known to have */
1259
  /* the smaller exponent), which also may need to be tested for zero */
1260
  /* for the weird IEEE 754 sign rules */
1261
  if (DFISZERO(dfl)) {
1262
    decCanonical(result, dfr);               /* clean copy */
1263
    /* "When the sum of two operands with opposite signs is */
1264
    /* exactly zero, the sign of that sum shall be '+' in all */
1265
    /* rounding modes except round toward -Infinity, in which */
1266
    /* mode that sign shall be '-'." */
1267
    if (diffsign && DFISZERO(result)) {
1268
      DFWORD(result, 0)&=~DECFLOAT_Sign;     /* assume sign 0 */
1269
      if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
1270
      }
1271
    return result;
1272
    } /* numfl is zero */
1273
  /* [here, LHS is non-zero; code below assumes that] */
1274
 
1275
  /* Coefficients layout during the calculations to follow: */
1276
  /* */
1277
  /*       Overlap case: */
1278
  /*       +------------------------------------------------+ */
1279
  /* acc:  |0000|      coeffa      | tail B |               | */
1280
  /*       +------------------------------------------------+ */
1281
  /* buf:  |0000| pad0s |      coeffb       |               | */
1282
  /*       +------------------------------------------------+ */
1283
  /* */
1284
  /*       Touching coefficients or gap: */
1285
  /*       +------------------------------------------------+ */
1286
  /* acc:  |0000|      coeffa      | gap |      coeffb      | */
1287
  /*       +------------------------------------------------+ */
1288
  /*       [buf not used or needed; gap clamped to Pmax] */
1289
 
1290
  /* lay out lhs coefficient into accumulator; this starts at acc+4 */
1291
  /* for decDouble or acc+6 for decQuad so the LSD is word- */
1292
  /* aligned; the top word gap is there only in case a carry digit */
1293
  /* is prefixed after the add -- it does not need to be zeroed */
1294
  #if DOUBLE
1295
    #define COFF 4                      /* offset into acc */
1296
  #elif QUAD
1297
    UBFROMUS(acc+4, 0);          /* prefix 00 */
1298
    #define COFF 6                      /* offset into acc */
1299
  #endif
1300
 
1301
  GETCOEFF(dfl, acc+COFF);              /* decode from decFloat */
1302
  ulsd=acc+COFF+DECPMAX-1;
1303
  umsd=acc+4;                           /* [having this here avoids */
1304
 
1305
  #if DECTRACE
1306
  {bcdnum tum;
1307
  tum.msd=umsd;
1308
  tum.lsd=ulsd;
1309
  tum.exponent=bexpl-DECBIAS;
1310
  tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
1311
  decShowNum(&tum, "dflx");}
1312
  #endif
1313
 
1314
  /* if signs differ, take ten's complement of lhs (here the */
1315
  /* coefficient is subtracted from all-nines; the 1 is added during */
1316
  /* the later add cycle -- zeros to the right do not matter because */
1317
  /* the complement of zero is zero); these are fixed-length inverts */
1318
  /* where the lsd is known to be at a 4-byte boundary (so no borrow */
1319
  /* possible) */
1320
  carry=0;                               /* assume no carry */
1321
  if (diffsign) {
1322
    carry=CARRYPAT;                     /* for +1 during add */
1323
    UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4));
1324
    UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8));
1325
    UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12));
1326
    UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16));
1327
    #if QUAD
1328
    UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20));
1329
    UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24));
1330
    UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28));
1331
    UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32));
1332
    UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36));
1333
    #endif
1334
    } /* diffsign */
1335
 
1336
  /* now process the rhs coefficient; if it cannot overlap lhs then */
1337
  /* it can be put straight into acc (with an appropriate gap, if */
1338
  /* needed) because no actual addition will be needed (except */
1339
  /* possibly to complete ten's complement) */
1340
  overlap=DECPMAX-(bexpl-bexpr);
1341
  #if DECTRACE
1342
  printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS));
1343
  printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
1344
  #endif
1345
 
1346
  if (overlap<=0) {                      /* no overlap possible */
1347
    uInt gap;                           /* local work */
1348
    /* since a full addition is not needed, a ten's complement */
1349
    /* calculation started above may need to be completed */
1350
    if (carry) {
1351
      for (ub=ulsd; *ub==9; ub--) *ub=0;
1352
      *ub+=1;
1353
      carry=0;                           /* taken care of */
1354
      }
1355
    /* up to DECPMAX-1 digits of the final result can extend down */
1356
    /* below the LSD of the lhs, so if the gap is >DECPMAX then the */
1357
    /* rhs will be simply sticky bits.  In this case the gap is */
1358
    /* clamped to DECPMAX and the exponent adjusted to suit [this is */
1359
    /* safe because the lhs is non-zero]. */
1360
    gap=-overlap;
1361
    if (gap>DECPMAX) {
1362
      bexpr+=gap-1;
1363
      gap=DECPMAX;
1364
      }
1365
    ub=ulsd+gap+1;                      /* where MSD will go */
1366
    /* Fill the gap with 0s; note that there is no addition to do */
1367
    ut=acc+COFF+DECPMAX;                /* start of gap */
1368
    for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */
1369
    if (overlap<-DECPMAX) {             /* gap was > DECPMAX */
1370
      *ub=(uByte)(!DFISZERO(dfr));      /* make sticky digit */
1371
      }
1372
     else {                             /* need full coefficient */
1373
      GETCOEFF(dfr, ub);                /* decode from decFloat */
1374
      ub+=DECPMAX-1;                    /* new LSD... */
1375
      }
1376
    ulsd=ub;                            /* save new LSD */
1377
    } /* no overlap possible */
1378
 
1379
   else {                               /* overlap>0 */
1380
    /* coefficients overlap (perhaps completely, although also */
1381
    /* perhaps only where zeros) */
1382
    if (overlap==DECPMAX) {             /* aligned */
1383
      ub=buf+COFF;                      /* where msd will go */
1384
      #if QUAD
1385
      UBFROMUS(buf+4, 0);                /* clear quad's 00 */
1386
      #endif
1387
      GETCOEFF(dfr, ub);                /* decode from decFloat */
1388
      }
1389
     else {                             /* unaligned */
1390
      ub=buf+COFF+DECPMAX-overlap;      /* where MSD will go */
1391
      /* Fill the prefix gap with 0s; 8 will cover most common */
1392
      /* unalignments, so start with direct assignments (a loop is */
1393
      /* then used for any remaining -- the loop (and the one in a */
1394
      /* moment) is not then on the critical path because the number */
1395
      /* of additions is reduced by (at least) two in this case) */
1396
      UBFROMUI(buf+4, 0);                /* [clears decQuad 00 too] */
1397
      UBFROMUI(buf+8, 0);
1398
      if (ub>buf+12) {
1399
        ut=buf+12;                      /* start any remaining */
1400
        for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */
1401
        }
1402
      GETCOEFF(dfr, ub);                /* decode from decFloat */
1403
 
1404
      /* now move tail of rhs across to main acc; again use direct */
1405
      /* copies for 8 digits-worth */
1406
      UBFROMUI(acc+COFF+DECPMAX,   UBTOUI(buf+COFF+DECPMAX));
1407
      UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4));
1408
      if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
1409
        us=buf+COFF+DECPMAX+8;          /* source */
1410
        ut=acc+COFF+DECPMAX+8;          /* target */
1411
        for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us));
1412
        }
1413
      } /* unaligned */
1414
 
1415
    ulsd=acc+(ub-buf+DECPMAX-1);        /* update LSD pointer */
1416
 
1417
    /* Now do the add of the non-tail; this is all nicely aligned, */
1418
    /* and is over a multiple of four digits (because for Quad two */
1419
    /* zero digits were added on the left); words in both acc and */
1420
    /* buf (buf especially) will often be zero */
1421
    /* [byte-by-byte add, here, is about 15% slower total effect than */
1422
    /* the by-fours] */
1423
 
1424
    /* Now effect the add; this is harder on a little-endian */
1425
    /* machine as the inter-digit carry cannot use the usual BCD */
1426
    /* addition trick because the bytes are loaded in the wrong order */
1427
    /* [this loop could be unrolled, but probably scarcely worth it] */
1428
 
1429
    ut=acc+COFF+DECPMAX-4;              /* target LSW (acc) */
1430
    us=buf+COFF+DECPMAX-4;              /* source LSW (buf, to add to acc) */
1431
 
1432
    #if !DECLITEND
1433
    for (; ut>=acc+4; ut-=4, us-=4) {   /* big-endian add loop */
1434
      /* bcd8 add */
1435
      carry+=UBTOUI(us);                /* rhs + carry */
1436
      if (carry==0) continue;            /* no-op */
1437
      carry+=UBTOUI(ut);                /* lhs */
1438
      /* Big-endian BCD adjust (uses internal carry) */
1439
      carry+=0x76f6f6f6;                /* note top nibble not all bits */
1440
      /* apply BCD adjust and save */
1441
      UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4));
1442
      carry>>=31;                       /* true carry was at far left */
1443
      } /* add loop */
1444
    #else
1445
    for (; ut>=acc+4; ut-=4, us-=4) {   /* little-endian add loop */
1446
      /* bcd8 add */
1447
      carry+=UBTOUI(us);                /* rhs + carry */
1448
      if (carry==0) continue;            /* no-op [common if unaligned] */
1449
      carry+=UBTOUI(ut);                /* lhs */
1450
      /* Little-endian BCD adjust; inter-digit carry must be manual */
1451
      /* because the lsb from the array will be in the most-significant */
1452
      /* byte of carry */
1453
      carry+=0x76767676;                /* note no inter-byte carries */
1454
      carry+=(carry & 0x80000000)>>15;
1455
      carry+=(carry & 0x00800000)>>15;
1456
      carry+=(carry & 0x00008000)>>15;
1457
      carry-=(carry & 0x60606060)>>4;   /* BCD adjust back */
1458
      UBFROMUI(ut, carry & 0x0f0f0f0f); /* clear debris and save */
1459
      /* here, final carry-out bit is at 0x00000080; move it ready */
1460
      /* for next word-add (i.e., to 0x01000000) */
1461
      carry=(carry & 0x00000080)<<17;
1462
      } /* add loop */
1463
    #endif
1464
 
1465
    #if DECTRACE
1466
    {bcdnum tum;
1467
    printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
1468
    tum.msd=umsd;  /* acc+4; */
1469
    tum.lsd=ulsd;
1470
    tum.exponent=0;
1471
    tum.sign=0;
1472
    decShowNum(&tum, "dfadd");}
1473
    #endif
1474
    } /* overlap possible */
1475
 
1476
  /* ordering here is a little strange in order to have slowest path */
1477
  /* first in GCC asm listing */
1478
  if (diffsign) {                  /* subtraction */
1479
    if (!carry) {                  /* no carry out means RHS<LHS */
1480
      /* borrowed -- take ten's complement */
1481
      /* sign is lhs sign */
1482
      num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
1483
 
1484
      /* invert the coefficient first by fours, then add one; space */
1485
      /* at the end of the buffer ensures the by-fours is always */
1486
      /* safe, but lsd+1 must be cleared to prevent a borrow */
1487
      /* if big-endian */
1488
      #if !DECLITEND
1489
      *(ulsd+1)=0;
1490
      #endif
1491
      /* there are always at least four coefficient words */
1492
      UBFROMUI(umsd,    0x09090909-UBTOUI(umsd));
1493
      UBFROMUI(umsd+4,  0x09090909-UBTOUI(umsd+4));
1494
      UBFROMUI(umsd+8,  0x09090909-UBTOUI(umsd+8));
1495
      UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12));
1496
      #if DOUBLE
1497
        #define BNEXT 16
1498
      #elif QUAD
1499
        UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16));
1500
        UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20));
1501
        UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24));
1502
        UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28));
1503
        UBFROMUI(umsd+32, 0x09090909-UBTOUI(umsd+32));
1504
        #define BNEXT 36
1505
      #endif
1506
      if (ulsd>=umsd+BNEXT) {           /* unaligned */
1507
        /* eight will handle most unaligments for Double; 16 for Quad */
1508
        UBFROMUI(umsd+BNEXT,   0x09090909-UBTOUI(umsd+BNEXT));
1509
        UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4));
1510
        #if DOUBLE
1511
        #define BNEXTY (BNEXT+8)
1512
        #elif QUAD
1513
        UBFROMUI(umsd+BNEXT+8,  0x09090909-UBTOUI(umsd+BNEXT+8));
1514
        UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12));
1515
        #define BNEXTY (BNEXT+16)
1516
        #endif
1517
        if (ulsd>=umsd+BNEXTY) {        /* very unaligned */
1518
          ut=umsd+BNEXTY;               /* -> continue */
1519
          for (;;ut+=4) {
1520
            UBFROMUI(ut, 0x09090909-UBTOUI(ut)); /* invert four digits */
1521
            if (ut>=ulsd-3) break;      /* all done */
1522
            }
1523
          }
1524
        }
1525
      /* complete the ten's complement by adding 1 */
1526
      for (ub=ulsd; *ub==9; ub--) *ub=0;
1527
      *ub+=1;
1528
      } /* borrowed */
1529
 
1530
     else {                        /* carry out means RHS>=LHS */
1531
      num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
1532
      /* all done except for the special IEEE 754 exact-zero-result */
1533
      /* rule (see above); while testing for zero, strip leading */
1534
      /* zeros (which will save decFinalize doing it) (this is in */
1535
      /* diffsign path, so carry impossible and true umsd is */
1536
      /* acc+COFF) */
1537
 
1538
      /* Check the initial coefficient area using the fast macro; */
1539
      /* this will often be all that needs to be done (as on the */
1540
      /* worst-case path when the subtraction was aligned and */
1541
      /* full-length) */
1542
      if (ISCOEFFZERO(acc+COFF)) {
1543
        umsd=acc+COFF+DECPMAX-1;   /* so far, so zero */
1544
        if (ulsd>umsd) {           /* more to check */
1545
          umsd++;                  /* to align after checked area */
1546
          for (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4;
1547
          for (; *umsd==0 && umsd<ulsd;) umsd++;
1548
          }
1549
        if (*umsd==0) {     /* must be true zero (and diffsign) */
1550
          num.sign=0;               /* assume + */
1551
          if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
1552
          }
1553
        }
1554
      /* [else was not zero, might still have leading zeros] */
1555
      } /* subtraction gave positive result */
1556
    } /* diffsign */
1557
 
1558
   else { /* same-sign addition */
1559
    num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
1560
    #if DOUBLE
1561
    if (carry) {                   /* only possible with decDouble */
1562
      *(acc+3)=1;                  /* [Quad has leading 00] */
1563
      umsd=acc+3;
1564
      }
1565
    #endif
1566
    } /* same sign */
1567
 
1568
  num.msd=umsd;                    /* set MSD .. */
1569
  num.lsd=ulsd;                    /* .. and LSD */
1570
  num.exponent=bexpr-DECBIAS;      /* set exponent to smaller, unbiassed */
1571
 
1572
  #if DECTRACE
1573
  decFloatShow(dfl, "dfl");
1574
  decFloatShow(dfr, "dfr");
1575
  decShowNum(&num, "postadd");
1576
  #endif
1577
  return decFinalize(result, &num, set); /* round, check, and lay out */
1578
  } /* decFloatAdd */
1579
 
1580
/* ------------------------------------------------------------------ */
1581
/* decFloatAnd -- logical digitwise AND of two decFloats              */
1582
/*                                                                    */
1583
/*   result gets the result of ANDing dfl and dfr                     */
1584
/*   dfl    is the first decFloat (lhs)                               */
1585
/*   dfr    is the second decFloat (rhs)                              */
1586
/*   set    is the context                                            */
1587
/*   returns result, which will be canonical with sign=0              */
1588
/*                                                                    */
1589
/* The operands must be positive, finite with exponent q=0, and       */
1590
/* comprise just zeros and ones; if not, Invalid operation results.   */
1591
/* ------------------------------------------------------------------ */
1592
decFloat * decFloatAnd(decFloat *result,
1593
                       const decFloat *dfl, const decFloat *dfr,
1594
                       decContext *set) {
1595
  if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
1596
   || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
1597
  /* the operands are positive finite integers (q=0) with just 0s and 1s */
1598
  #if DOUBLE
1599
   DFWORD(result, 0)=ZEROWORD
1600
                   |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
1601
   DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
1602
  #elif QUAD
1603
   DFWORD(result, 0)=ZEROWORD
1604
                   |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
1605
   DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
1606
   DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
1607
   DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
1608
  #endif
1609
  return result;
1610
  } /* decFloatAnd */
1611
 
1612
/* ------------------------------------------------------------------ */
1613
/* decFloatCanonical -- copy a decFloat, making canonical             */
1614
/*                                                                    */
1615
/*   result gets the canonicalized df                                 */
1616
/*   df     is the decFloat to copy and make canonical                */
1617
/*   returns result                                                   */
1618
/*                                                                    */
1619
/* This works on specials, too; no error or exception is possible.    */
1620
/* ------------------------------------------------------------------ */
1621
decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
1622
  return decCanonical(result, df);
1623
  } /* decFloatCanonical */
1624
 
1625
/* ------------------------------------------------------------------ */
1626
/* decFloatClass -- return the class of a decFloat                    */
1627
/*                                                                    */
1628
/*   df is the decFloat to test                                       */
1629
/*   returns the decClass that df falls into                          */
1630
/* ------------------------------------------------------------------ */
1631
enum decClass decFloatClass(const decFloat *df) {
1632
  Int exp;                         /* exponent */
1633
  if (DFISSPECIAL(df)) {
1634
    if (DFISQNAN(df)) return DEC_CLASS_QNAN;
1635
    if (DFISSNAN(df)) return DEC_CLASS_SNAN;
1636
    /* must be an infinity */
1637
    if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
1638
    return DEC_CLASS_POS_INF;
1639
    }
1640
  if (DFISZERO(df)) {              /* quite common */
1641
    if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
1642
    return DEC_CLASS_POS_ZERO;
1643
    }
1644
  /* is finite and non-zero; similar code to decFloatIsNormal, here */
1645
  /* [this could be speeded up slightly by in-lining decFloatDigits] */
1646
  exp=GETEXPUN(df)                 /* get unbiased exponent .. */
1647
     +decFloatDigits(df)-1;        /* .. and make adjusted exponent */
1648
  if (exp>=DECEMIN) {              /* is normal */
1649
    if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
1650
    return DEC_CLASS_POS_NORMAL;
1651
    }
1652
  /* is subnormal */
1653
  if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
1654
  return DEC_CLASS_POS_SUBNORMAL;
1655
  } /* decFloatClass */
1656
 
1657
/* ------------------------------------------------------------------ */
1658
/* decFloatClassString -- return the class of a decFloat as a string  */
1659
/*                                                                    */
1660
/*   df is the decFloat to test                                       */
1661
/*   returns a constant string describing the class df falls into     */
1662
/* ------------------------------------------------------------------ */
1663
const char *decFloatClassString(const decFloat *df) {
1664
  enum decClass eclass=decFloatClass(df);
1665
  if (eclass==DEC_CLASS_POS_NORMAL)    return DEC_ClassString_PN;
1666
  if (eclass==DEC_CLASS_NEG_NORMAL)    return DEC_ClassString_NN;
1667
  if (eclass==DEC_CLASS_POS_ZERO)      return DEC_ClassString_PZ;
1668
  if (eclass==DEC_CLASS_NEG_ZERO)      return DEC_ClassString_NZ;
1669
  if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
1670
  if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
1671
  if (eclass==DEC_CLASS_POS_INF)       return DEC_ClassString_PI;
1672
  if (eclass==DEC_CLASS_NEG_INF)       return DEC_ClassString_NI;
1673
  if (eclass==DEC_CLASS_QNAN)          return DEC_ClassString_QN;
1674
  if (eclass==DEC_CLASS_SNAN)          return DEC_ClassString_SN;
1675
  return DEC_ClassString_UN;           /* Unknown */
1676
  } /* decFloatClassString */
1677
 
1678
/* ------------------------------------------------------------------ */
1679
/* decFloatCompare -- compare two decFloats; quiet NaNs allowed       */
1680
/*                                                                    */
1681
/*   result gets the result of comparing dfl and dfr                  */
1682
/*   dfl    is the first decFloat (lhs)                               */
1683
/*   dfr    is the second decFloat (rhs)                              */
1684
/*   set    is the context                                            */
1685
/*   returns result, which may be -1, 0, 1, or NaN (Unordered)        */
1686
/* ------------------------------------------------------------------ */
1687
decFloat * decFloatCompare(decFloat *result,
1688
                           const decFloat *dfl, const decFloat *dfr,
1689
                           decContext *set) {
1690
  Int comp;                                  /* work */
1691
  /* NaNs are handled as usual */
1692
  if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
1693
  /* numeric comparison needed */
1694
  comp=decNumCompare(dfl, dfr, 0);
1695
  decFloatZero(result);
1696
  if (comp==0) return result;
1697
  DFBYTE(result, DECBYTES-1)=0x01;      /* LSD=1 */
1698
  if (comp<0) DFBYTE(result, 0)|=0x80;    /* set sign bit */
1699
  return result;
1700
  } /* decFloatCompare */
1701
 
1702
/* ------------------------------------------------------------------ */
1703
/* decFloatCompareSignal -- compare two decFloats; all NaNs signal    */
1704
/*                                                                    */
1705
/*   result gets the result of comparing dfl and dfr                  */
1706
/*   dfl    is the first decFloat (lhs)                               */
1707
/*   dfr    is the second decFloat (rhs)                              */
1708
/*   set    is the context                                            */
1709
/*   returns result, which may be -1, 0, 1, or NaN (Unordered)        */
1710
/* ------------------------------------------------------------------ */
1711
decFloat * decFloatCompareSignal(decFloat *result,
1712
                                 const decFloat *dfl, const decFloat *dfr,
1713
                                 decContext *set) {
1714
  Int comp;                                  /* work */
1715
  /* NaNs are handled as usual, except that all NaNs signal */
1716
  if (DFISNAN(dfl) || DFISNAN(dfr)) {
1717
    set->status|=DEC_Invalid_operation;
1718
    return decNaNs(result, dfl, dfr, set);
1719
    }
1720
  /* numeric comparison needed */
1721
  comp=decNumCompare(dfl, dfr, 0);
1722
  decFloatZero(result);
1723
  if (comp==0) return result;
1724
  DFBYTE(result, DECBYTES-1)=0x01;      /* LSD=1 */
1725
  if (comp<0) DFBYTE(result, 0)|=0x80;    /* set sign bit */
1726
  return result;
1727
  } /* decFloatCompareSignal */
1728
 
1729
/* ------------------------------------------------------------------ */
1730
/* decFloatCompareTotal -- compare two decFloats with total ordering  */
1731
/*                                                                    */
1732
/*   result gets the result of comparing dfl and dfr                  */
1733
/*   dfl    is the first decFloat (lhs)                               */
1734
/*   dfr    is the second decFloat (rhs)                              */
1735
/*   returns result, which may be -1, 0, or 1                         */
1736
/* ------------------------------------------------------------------ */
1737
decFloat * decFloatCompareTotal(decFloat *result,
1738
                                const decFloat *dfl, const decFloat *dfr) {
1739
  Int  comp;                                 /* work */
1740
  uInt uiwork;                               /* for macros */
1741
  #if QUAD
1742
  uShort uswork;                             /* .. */
1743
  #endif
1744
  if (DFISNAN(dfl) || DFISNAN(dfr)) {
1745
    Int nanl, nanr;                          /* work */
1746
    /* morph NaNs to +/- 1 or 2, leave numbers as 0 */
1747
    nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2;      /* quiet > signalling */
1748
    if (DFISSIGNED(dfl)) nanl=-nanl;
1749
    nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
1750
    if (DFISSIGNED(dfr)) nanr=-nanr;
1751
    if (nanl>nanr) comp=+1;
1752
     else if (nanl<nanr) comp=-1;
1753
     else { /* NaNs are the same type and sign .. must compare payload */
1754
      /* buffers need +2 for QUAD */
1755
      uByte bufl[DECPMAX+4];                 /* for LHS coefficient + foot */
1756
      uByte bufr[DECPMAX+4];                 /* for RHS coefficient + foot */
1757
      uByte *ub, *uc;                        /* work */
1758
      Int sigl;                              /* signum of LHS */
1759
      sigl=(DFISSIGNED(dfl) ? -1 : +1);
1760
 
1761
      /* decode the coefficients */
1762
      /* (shift both right two if Quad to make a multiple of four) */
1763
      #if QUAD
1764
        UBFROMUS(bufl, 0);
1765
        UBFROMUS(bufr, 0);
1766
      #endif
1767
      GETCOEFF(dfl, bufl+QUAD*2);            /* decode from decFloat */
1768
      GETCOEFF(dfr, bufr+QUAD*2);            /* .. */
1769
      /* all multiples of four, here */
1770
      comp=0;                                 /* assume equal */
1771
      for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
1772
        uInt ui=UBTOUI(ub);
1773
        if (ui==UBTOUI(uc)) continue; /* so far so same */
1774
        /* about to find a winner; go by bytes in case little-endian */
1775
        for (;; ub++, uc++) {
1776
          if (*ub==*uc) continue;
1777
          if (*ub>*uc) comp=sigl;            /* difference found */
1778
           else comp=-sigl;                  /* .. */
1779
           break;
1780
          }
1781
        }
1782
      } /* same NaN type and sign */
1783
    }
1784
   else {
1785
    /* numeric comparison needed */
1786
    comp=decNumCompare(dfl, dfr, 1);    /* total ordering */
1787
    }
1788
  decFloatZero(result);
1789
  if (comp==0) return result;
1790
  DFBYTE(result, DECBYTES-1)=0x01;      /* LSD=1 */
1791
  if (comp<0) DFBYTE(result, 0)|=0x80;    /* set sign bit */
1792
  return result;
1793
  } /* decFloatCompareTotal */
1794
 
1795
/* ------------------------------------------------------------------ */
1796
/* decFloatCompareTotalMag -- compare magnitudes with total ordering  */
1797
/*                                                                    */
1798
/*   result gets the result of comparing abs(dfl) and abs(dfr)        */
1799
/*   dfl    is the first decFloat (lhs)                               */
1800
/*   dfr    is the second decFloat (rhs)                              */
1801
/*   returns result, which may be -1, 0, or 1                         */
1802
/* ------------------------------------------------------------------ */
1803
decFloat * decFloatCompareTotalMag(decFloat *result,
1804
                                const decFloat *dfl, const decFloat *dfr) {
1805
  decFloat a, b;                        /* for copy if needed */
1806
  /* copy and redirect signed operand(s) */
1807
  if (DFISSIGNED(dfl)) {
1808
    decFloatCopyAbs(&a, dfl);
1809
    dfl=&a;
1810
    }
1811
  if (DFISSIGNED(dfr)) {
1812
    decFloatCopyAbs(&b, dfr);
1813
    dfr=&b;
1814
    }
1815
  return decFloatCompareTotal(result, dfl, dfr);
1816
  } /* decFloatCompareTotalMag */
1817
 
1818
/* ------------------------------------------------------------------ */
1819
/* decFloatCopy -- copy a decFloat as-is                              */
1820
/*                                                                    */
1821
/*   result gets the copy of dfl                                      */
1822
/*   dfl    is the decFloat to copy                                   */
1823
/*   returns result                                                   */
1824
/*                                                                    */
1825
/* This is a bitwise operation; no errors or exceptions are possible. */
1826
/* ------------------------------------------------------------------ */
1827
decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
1828
  if (dfl!=result) *result=*dfl;             /* copy needed */
1829
  return result;
1830
  } /* decFloatCopy */
1831
 
1832
/* ------------------------------------------------------------------ */
1833
/* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0     */
1834
/*                                                                    */
1835
/*   result gets the copy of dfl with sign bit 0                      */
1836
/*   dfl    is the decFloat to copy                                   */
1837
/*   returns result                                                   */
1838
/*                                                                    */
1839
/* This is a bitwise operation; no errors or exceptions are possible. */
1840
/* ------------------------------------------------------------------ */
1841
decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
1842
  if (dfl!=result) *result=*dfl;        /* copy needed */
1843
  DFBYTE(result, 0)&=~0x80;              /* zero sign bit */
1844
  return result;
1845
  } /* decFloatCopyAbs */
1846
 
1847
/* ------------------------------------------------------------------ */
1848
/* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
1849
/*                                                                    */
1850
/*   result gets the copy of dfl with sign bit inverted               */
1851
/*   dfl    is the decFloat to copy                                   */
1852
/*   returns result                                                   */
1853
/*                                                                    */
1854
/* This is a bitwise operation; no errors or exceptions are possible. */
1855
/* ------------------------------------------------------------------ */
1856
decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
1857
  if (dfl!=result) *result=*dfl;        /* copy needed */
1858
  DFBYTE(result, 0)^=0x80;               /* invert sign bit */
1859
  return result;
1860
  } /* decFloatCopyNegate */
1861
 
1862
/* ------------------------------------------------------------------ */
1863
/* decFloatCopySign -- copy a decFloat with the sign of another       */
1864
/*                                                                    */
1865
/*   result gets the result of copying dfl with the sign of dfr       */
1866
/*   dfl    is the first decFloat (lhs)                               */
1867
/*   dfr    is the second decFloat (rhs)                              */
1868
/*   returns result                                                   */
1869
/*                                                                    */
1870
/* This is a bitwise operation; no errors or exceptions are possible. */
1871
/* ------------------------------------------------------------------ */
1872
decFloat * decFloatCopySign(decFloat *result,
1873
                            const decFloat *dfl, const decFloat *dfr) {
1874
  uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80);   /* save sign bit */
1875
  if (dfl!=result) *result=*dfl;             /* copy needed */
1876
  DFBYTE(result, 0)&=~0x80;                   /* clear sign .. */
1877
  DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
1878
  return result;
1879
  } /* decFloatCopySign */
1880
 
1881
/* ------------------------------------------------------------------ */
1882
/* decFloatDigits -- return the number of digits in a decFloat        */
1883
/*                                                                    */
1884
/*   df is the decFloat to investigate                                */
1885
/*   returns the number of significant digits in the decFloat; a      */
1886
/*     zero coefficient returns 1 as does an infinity (a NaN returns  */
1887
/*     the number of digits in the payload)                           */
1888
/* ------------------------------------------------------------------ */
1889
/* private macro to extract a declet according to provided formula */
1890
/* (form), and if it is non-zero then return the calculated digits */
1891
/* depending on the declet number (n), where n=0 for the most */
1892
/* significant declet; uses uInt dpd for work */
1893
#define dpdlenchk(n, form) {dpd=(form)&0x3ff;     \
1894
  if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1895
/* next one is used when it is known that the declet must be */
1896
/* non-zero, or is the final zero declet */
1897
#define dpdlendun(n, form) {dpd=(form)&0x3ff;     \
1898
  if (dpd==0) return 1;                    \
1899
  return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
1900
 
1901
uInt decFloatDigits(const decFloat *df) {
1902
  uInt dpd;                        /* work */
1903
  uInt sourhi=DFWORD(df, 0);        /* top word from source decFloat */
1904
  #if QUAD
1905
  uInt sourmh, sourml;
1906
  #endif
1907
  uInt sourlo;
1908
 
1909
  if (DFISINF(df)) return 1;
1910
  /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
1911
  /* then the coefficient is full-length */
1912
  if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
1913
 
1914
  #if DOUBLE
1915
    if (sourhi&0x0003ffff) {     /* ends in first */
1916
      dpdlenchk(0, sourhi>>8);
1917
      sourlo=DFWORD(df, 1);
1918
      dpdlendun(1, (sourhi<<2) | (sourlo>>30));
1919
      } /* [cannot drop through] */
1920
    sourlo=DFWORD(df, 1);  /* sourhi not involved now */
1921
    if (sourlo&0xfff00000) {     /* in one of first two */
1922
      dpdlenchk(1, sourlo>>30);  /* very rare */
1923
      dpdlendun(2, sourlo>>20);
1924
      } /* [cannot drop through] */
1925
    dpdlenchk(3, sourlo>>10);
1926
    dpdlendun(4, sourlo);
1927
    /* [cannot drop through] */
1928
 
1929
  #elif QUAD
1930
    if (sourhi&0x00003fff) {     /* ends in first */
1931
      dpdlenchk(0, sourhi>>4);
1932
      sourmh=DFWORD(df, 1);
1933
      dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
1934
      } /* [cannot drop through] */
1935
    sourmh=DFWORD(df, 1);
1936
    if (sourmh) {
1937
      dpdlenchk(1, sourmh>>26);
1938
      dpdlenchk(2, sourmh>>16);
1939
      dpdlenchk(3, sourmh>>6);
1940
      sourml=DFWORD(df, 2);
1941
      dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
1942
      } /* [cannot drop through] */
1943
    sourml=DFWORD(df, 2);
1944
    if (sourml) {
1945
      dpdlenchk(4, sourml>>28);
1946
      dpdlenchk(5, sourml>>18);
1947
      dpdlenchk(6, sourml>>8);
1948
      sourlo=DFWORD(df, 3);
1949
      dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
1950
      } /* [cannot drop through] */
1951
    sourlo=DFWORD(df, 3);
1952
    if (sourlo&0xfff00000) {     /* in one of first two */
1953
      dpdlenchk(7, sourlo>>30);  /* very rare */
1954
      dpdlendun(8, sourlo>>20);
1955
      } /* [cannot drop through] */
1956
    dpdlenchk(9, sourlo>>10);
1957
    dpdlendun(10, sourlo);
1958
    /* [cannot drop through] */
1959
  #endif
1960
  } /* decFloatDigits */
1961
 
1962
/* ------------------------------------------------------------------ */
1963
/* decFloatDivide -- divide a decFloat by another                     */
1964
/*                                                                    */
1965
/*   result gets the result of dividing dfl by dfr:                   */
1966
/*   dfl    is the first decFloat (lhs)                               */
1967
/*   dfr    is the second decFloat (rhs)                              */
1968
/*   set    is the context                                            */
1969
/*   returns result                                                   */
1970
/*                                                                    */
1971
/* ------------------------------------------------------------------ */
1972
/* This is just a wrapper. */
1973
decFloat * decFloatDivide(decFloat *result,
1974
                          const decFloat *dfl, const decFloat *dfr,
1975
                          decContext *set) {
1976
  return decDivide(result, dfl, dfr, set, DIVIDE);
1977
  } /* decFloatDivide */
1978
 
1979
/* ------------------------------------------------------------------ */
1980
/* decFloatDivideInteger -- integer divide a decFloat by another      */
1981
/*                                                                    */
1982
/*   result gets the result of dividing dfl by dfr:                   */
1983
/*   dfl    is the first decFloat (lhs)                               */
1984
/*   dfr    is the second decFloat (rhs)                              */
1985
/*   set    is the context                                            */
1986
/*   returns result                                                   */
1987
/*                                                                    */
1988
/* ------------------------------------------------------------------ */
1989
decFloat * decFloatDivideInteger(decFloat *result,
1990
                             const decFloat *dfl, const decFloat *dfr,
1991
                             decContext *set) {
1992
  return decDivide(result, dfl, dfr, set, DIVIDEINT);
1993
  } /* decFloatDivideInteger */
1994
 
1995
/* ------------------------------------------------------------------ */
1996
/* decFloatFMA -- multiply and add three decFloats, fused             */
1997
/*                                                                    */
1998
/*   result gets the result of (dfl*dfr)+dff with a single rounding   */
1999
/*   dfl    is the first decFloat (lhs)                               */
2000
/*   dfr    is the second decFloat (rhs)                              */
2001
/*   dff    is the final decFloat (fhs)                               */
2002
/*   set    is the context                                            */
2003
/*   returns result                                                   */
2004
/*                                                                    */
2005
/* ------------------------------------------------------------------ */
2006
decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
2007
                       const decFloat *dfr, const decFloat *dff,
2008
                       decContext *set) {
2009
 
2010
  /* The accumulator has the bytes needed for FiniteMultiply, plus */
2011
  /* one byte to the left in case of carry, plus DECPMAX+2 to the */
2012
  /* right for the final addition (up to full fhs + round & sticky) */
2013
  #define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2))
2014
  uByte  acc[FMALEN];              /* for multiplied coefficient in BCD */
2015
                                   /* .. and for final result */
2016
  bcdnum mul;                      /* for multiplication result */
2017
  bcdnum fin;                      /* for final operand, expanded */
2018
  uByte  coe[ROUNDUP4(DECPMAX)];   /* dff coefficient in BCD */
2019
  bcdnum *hi, *lo;                 /* bcdnum with higher/lower exponent */
2020
  uInt   diffsign;                 /* non-zero if signs differ */
2021
  uInt   hipad;                    /* pad digit for hi if needed */
2022
  Int    padding;                  /* excess exponent */
2023
  uInt   carry;                    /* +1 for ten's complement and during add */
2024
  uByte  *ub, *uh, *ul;            /* work */
2025
  uInt   uiwork;                   /* for macros */
2026
 
2027
  /* handle all the special values [any special operand leads to a */
2028
  /* special result] */
2029
  if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
2030
    decFloat proxy;                /* multiplication result proxy */
2031
    /* NaNs are handled as usual, giving priority to sNaNs */
2032
    if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2033
    if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
2034
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2035
    if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
2036
    /* One or more of the three is infinite */
2037
    /* infinity times zero is bad */
2038
    decFloatZero(&proxy);
2039
    if (DFISINF(dfl)) {
2040
      if (DFISZERO(dfr)) return decInvalid(result, set);
2041
      decInfinity(&proxy, &proxy);
2042
      }
2043
     else if (DFISINF(dfr)) {
2044
      if (DFISZERO(dfl)) return decInvalid(result, set);
2045
      decInfinity(&proxy, &proxy);
2046
      }
2047
    /* compute sign of multiplication and place in proxy */
2048
    DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
2049
    if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
2050
    /* dff is Infinite */
2051
    if (!DFISINF(&proxy)) return decInfinity(result, dff);
2052
    /* both sides of addition are infinite; different sign is bad */
2053
    if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
2054
      return decInvalid(result, set);
2055
    return decFloatCopy(result, &proxy);
2056
    }
2057
 
2058
  /* Here when all operands are finite */
2059
 
2060
  /* First multiply dfl*dfr */
2061
  decFiniteMultiply(&mul, acc+1, dfl, dfr);
2062
  /* The multiply is complete, exact and unbounded, and described in */
2063
  /* mul with the coefficient held in acc[1...] */
2064
 
2065
  /* now add in dff; the algorithm is essentially the same as */
2066
  /* decFloatAdd, but the code is different because the code there */
2067
  /* is highly optimized for adding two numbers of the same size */
2068
  fin.exponent=GETEXPUN(dff);           /* get dff exponent and sign */
2069
  fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
2070
  diffsign=mul.sign^fin.sign;           /* note if signs differ */
2071
  fin.msd=coe;
2072
  fin.lsd=coe+DECPMAX-1;
2073
  GETCOEFF(dff, coe);                   /* extract the coefficient */
2074
 
2075
  /* now set hi and lo so that hi points to whichever of mul and fin */
2076
  /* has the higher exponent and lo points to the other [don't care, */
2077
  /* if the same].  One coefficient will be in acc, the other in coe. */
2078
  if (mul.exponent>=fin.exponent) {
2079
    hi=&mul;
2080
    lo=&fin;
2081
    }
2082
   else {
2083
    hi=&fin;
2084
    lo=&mul;
2085
    }
2086
 
2087
  /* remove leading zeros on both operands; this will save time later */
2088
  /* and make testing for zero trivial (tests are safe because acc */
2089
  /* and coe are rounded up to uInts) */
2090
  for (; UBTOUI(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
2091
  for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
2092
  for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
2093
  for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
2094
 
2095
  /* if hi is zero then result will be lo (which has the smaller */
2096
  /* exponent), which also may need to be tested for zero for the */
2097
  /* weird IEEE 754 sign rules */
2098
  if (*hi->msd==0) {                          /* hi is zero */
2099
    /* "When the sum of two operands with opposite signs is */
2100
    /* exactly zero, the sign of that sum shall be '+' in all */
2101
    /* rounding modes except round toward -Infinity, in which */
2102
    /* mode that sign shall be '-'." */
2103
    if (diffsign) {
2104
      if (*lo->msd==0) {              /* lo is zero */
2105
        lo->sign=0;
2106
        if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2107
        } /* diffsign && lo=0 */
2108
      } /* diffsign */
2109
    return decFinalize(result, lo, set);     /* may need clamping */
2110
    } /* numfl is zero */
2111
  /* [here, both are minimal length and hi is non-zero] */
2112
  /* (if lo is zero then padding with zeros may be needed, below) */
2113
 
2114
  /* if signs differ, take the ten's complement of hi (zeros to the */
2115
  /* right do not matter because the complement of zero is zero); the */
2116
  /* +1 is done later, as part of the addition, inserted at the */
2117
  /* correct digit */
2118
  hipad=0;
2119
  carry=0;
2120
  if (diffsign) {
2121
    hipad=9;
2122
    carry=1;
2123
    /* exactly the correct number of digits must be inverted */
2124
    for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(uh));
2125
    for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
2126
    }
2127
 
2128
  /* ready to add; note that hi has no leading zeros so gap */
2129
  /* calculation does not have to be as pessimistic as in decFloatAdd */
2130
  /* (this is much more like the arbitrary-precision algorithm in */
2131
  /* Rexx and decNumber) */
2132
 
2133
  /* padding is the number of zeros that would need to be added to hi */
2134
  /* for its lsd to be aligned with the lsd of lo */
2135
  padding=hi->exponent-lo->exponent;
2136
  /* printf("FMA pad %ld\n", (LI)padding); */
2137
 
2138
  /* the result of the addition will be built into the accumulator, */
2139
  /* starting from the far right; this could be either hi or lo, and */
2140
  /* will be aligned */
2141
  ub=acc+FMALEN-1;                 /* where lsd of result will go */
2142
  ul=lo->lsd;                      /* lsd of rhs */
2143
 
2144
  if (padding!=0) {                 /* unaligned */
2145
    /* if the msd of lo is more than DECPMAX+2 digits to the right of */
2146
    /* the original msd of hi then it can be reduced to a single */
2147
    /* digit at the right place, as it stays clear of hi digits */
2148
    /* [it must be DECPMAX+2 because during a subtraction the msd */
2149
    /* could become 0 after a borrow from 1.000 to 0.9999...] */
2150
 
2151
    Int hilen=(Int)(hi->lsd-hi->msd+1); /* length of hi */
2152
    Int lolen=(Int)(lo->lsd-lo->msd+1); /* and of lo */
2153
 
2154
    if (hilen+padding-lolen > DECPMAX+2) {   /* can reduce lo to single */
2155
      /* make sure it is virtually at least DECPMAX from hi->msd, at */
2156
      /* least to right of hi->lsd (in case of destructive subtract), */
2157
      /* and separated by at least two digits from either of those */
2158
      /* (the tricky DOUBLE case is when hi is a 1 that will become a */
2159
      /* 0.9999... by subtraction: */
2160
      /*   hi:   1                                   E+16 */
2161
      /*   lo:    .................1000000000000000  E-16 */
2162
      /* which for the addition pads to: */
2163
      /*   hi:   1000000000000000000                 E-16 */
2164
      /*   lo:    .................1000000000000000  E-16 */
2165
      Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
2166
 
2167
      /* printf("FMA reduce: %ld\n", (LI)reduce); */
2168
      lo->lsd=lo->msd;                       /* to single digit [maybe 0] */
2169
      lo->exponent=newexp;                   /* new lowest exponent */
2170
      padding=hi->exponent-lo->exponent;     /* recalculate */
2171
      ul=lo->lsd;                            /* .. and repoint */
2172
      }
2173
 
2174
    /* padding is still > 0, but will fit in acc (less leading carry slot) */
2175
    #if DECCHECK
2176
      if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
2177
      if (hilen+padding+1>FMALEN)
2178
        printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (LI)padding);
2179
      /* printf("FMA padding: %ld\n", (LI)padding); */
2180
    #endif
2181
 
2182
    /* padding digits can now be set in the result; one or more of */
2183
    /* these will come from lo; others will be zeros in the gap */
2184
    for (; ul-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) {
2185
      UBFROMUI(ub-3, UBTOUI(ul-3));          /* [cannot overlap] */
2186
      }
2187
    for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
2188
    for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
2189
    }
2190
 
2191
  /* addition now complete to the right of the rightmost digit of hi */
2192
  uh=hi->lsd;
2193
 
2194
  /* dow do the add from hi->lsd to the left */
2195
  /* [bytewise, because either operand can run out at any time] */
2196
  /* carry was set up depending on ten's complement above */
2197
  /* first assume both operands have some digits */
2198
  for (;; ub--) {
2199
    if (uh<hi->msd || ul<lo->msd) break;
2200
    *ub=(uByte)(carry+(*uh--)+(*ul--));
2201
    carry=0;
2202
    if (*ub<10) continue;
2203
    *ub-=10;
2204
    carry=1;
2205
    } /* both loop */
2206
 
2207
  if (ul<lo->msd) {           /* to left of lo */
2208
    for (;; ub--) {
2209
      if (uh<hi->msd) break;
2210
      *ub=(uByte)(carry+(*uh--));  /* [+0] */
2211
      carry=0;
2212
      if (*ub<10) continue;
2213
      *ub-=10;
2214
      carry=1;
2215
      } /* hi loop */
2216
    }
2217
   else {                     /* to left of hi */
2218
    for (;; ub--) {
2219
      if (ul<lo->msd) break;
2220
      *ub=(uByte)(carry+hipad+(*ul--));
2221
      carry=0;
2222
      if (*ub<10) continue;
2223
      *ub-=10;
2224
      carry=1;
2225
      } /* lo loop */
2226
    }
2227
 
2228
  /* addition complete -- now handle carry, borrow, etc. */
2229
  /* use lo to set up the num (its exponent is already correct, and */
2230
  /* sign usually is) */
2231
  lo->msd=ub+1;
2232
  lo->lsd=acc+FMALEN-1;
2233
  /* decShowNum(lo, "lo"); */
2234
  if (!diffsign) {                 /* same-sign addition */
2235
    if (carry) {                   /* carry out */
2236
      *ub=1;                       /* place the 1 .. */
2237
      lo->msd--;                   /* .. and update */
2238
      }
2239
    } /* same sign */
2240
   else {                          /* signs differed (subtraction) */
2241
    if (!carry) {                  /* no carry out means hi<lo */
2242
      /* borrowed -- take ten's complement of the right digits */
2243
      lo->sign=hi->sign;           /* sign is lhs sign */
2244
      for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(ul));
2245
      for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
2246
      /* complete the ten's complement by adding 1 [cannot overrun] */
2247
      for (ul--; *ul==9; ul--) *ul=0;
2248
      *ul+=1;
2249
      } /* borrowed */
2250
     else {                        /* carry out means hi>=lo */
2251
      /* sign to use is lo->sign */
2252
      /* all done except for the special IEEE 754 exact-zero-result */
2253
      /* rule (see above); while testing for zero, strip leading */
2254
      /* zeros (which will save decFinalize doing it) */
2255
      for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
2256
      for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
2257
      if (*lo->msd==0) {    /* must be true zero (and diffsign) */
2258
        lo->sign=0;                 /* assume + */
2259
        if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
2260
        }
2261
      /* [else was not zero, might still have leading zeros] */
2262
      } /* subtraction gave positive result */
2263
    } /* diffsign */
2264
 
2265
  #if DECCHECK
2266
  /* assert no left underrun */
2267
  if (lo->msd<acc) {
2268
    printf("FMA underrun by %ld \n", (LI)(acc-lo->msd));
2269
    }
2270
  #endif
2271
 
2272
  return decFinalize(result, lo, set);  /* round, check, and lay out */
2273
  } /* decFloatFMA */
2274
 
2275
/* ------------------------------------------------------------------ */
2276
/* decFloatFromInt -- initialise a decFloat from an Int               */
2277
/*                                                                    */
2278
/*   result gets the converted Int                                    */
2279
/*   n      is the Int to convert                                     */
2280
/*   returns result                                                   */
2281
/*                                                                    */
2282
/* The result is Exact; no errors or exceptions are possible.         */
2283
/* ------------------------------------------------------------------ */
2284
decFloat * decFloatFromInt32(decFloat *result, Int n) {
2285
  uInt u=(uInt)n;                       /* copy as bits */
2286
  uInt encode;                          /* work */
2287
  DFWORD(result, 0)=ZEROWORD;            /* always */
2288
  #if QUAD
2289
    DFWORD(result, 1)=0;
2290
    DFWORD(result, 2)=0;
2291
  #endif
2292
  if (n<0) {                             /* handle -n with care */
2293
    /* [This can be done without the test, but is then slightly slower] */
2294
    u=(~u)+1;
2295
    DFWORD(result, 0)|=DECFLOAT_Sign;
2296
    }
2297
  /* Since the maximum value of u now is 2**31, only the low word of */
2298
  /* result is affected */
2299
  encode=BIN2DPD[u%1000];
2300
  u/=1000;
2301
  encode|=BIN2DPD[u%1000]<<10;
2302
  u/=1000;
2303
  encode|=BIN2DPD[u%1000]<<20;
2304
  u/=1000;                              /* now 0, 1, or 2 */
2305
  encode|=u<<30;
2306
  DFWORD(result, DECWORDS-1)=encode;
2307
  return result;
2308
  } /* decFloatFromInt32 */
2309
 
2310
/* ------------------------------------------------------------------ */
2311
/* decFloatFromUInt -- initialise a decFloat from a uInt              */
2312
/*                                                                    */
2313
/*   result gets the converted uInt                                   */
2314
/*   n      is the uInt to convert                                    */
2315
/*   returns result                                                   */
2316
/*                                                                    */
2317
/* The result is Exact; no errors or exceptions are possible.         */
2318
/* ------------------------------------------------------------------ */
2319
decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
2320
  uInt encode;                          /* work */
2321
  DFWORD(result, 0)=ZEROWORD;            /* always */
2322
  #if QUAD
2323
    DFWORD(result, 1)=0;
2324
    DFWORD(result, 2)=0;
2325
  #endif
2326
  encode=BIN2DPD[u%1000];
2327
  u/=1000;
2328
  encode|=BIN2DPD[u%1000]<<10;
2329
  u/=1000;
2330
  encode|=BIN2DPD[u%1000]<<20;
2331
  u/=1000;                              /* now 0 -> 4 */
2332
  encode|=u<<30;
2333
  DFWORD(result, DECWORDS-1)=encode;
2334
  DFWORD(result, DECWORDS-2)|=u>>2;     /* rarely non-zero */
2335
  return result;
2336
  } /* decFloatFromUInt32 */
2337
 
2338
/* ------------------------------------------------------------------ */
2339
/* decFloatInvert -- logical digitwise INVERT of a decFloat           */
2340
/*                                                                    */
2341
/*   result gets the result of INVERTing df                           */
2342
/*   df     is the decFloat to invert                                 */
2343
/*   set    is the context                                            */
2344
/*   returns result, which will be canonical with sign=0              */
2345
/*                                                                    */
2346
/* The operand must be positive, finite with exponent q=0, and        */
2347
/* comprise just zeros and ones; if not, Invalid operation results.   */
2348
/* ------------------------------------------------------------------ */
2349
decFloat * decFloatInvert(decFloat *result, const decFloat *df,
2350
                          decContext *set) {
2351
  uInt sourhi=DFWORD(df, 0);             /* top word of dfs */
2352
 
2353
  if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
2354
  /* the operand is a finite integer (q=0) */
2355
  #if DOUBLE
2356
   DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
2357
   DFWORD(result, 1)=(~DFWORD(df, 1))   &0x49124491;
2358
  #elif QUAD
2359
   DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
2360
   DFWORD(result, 1)=(~DFWORD(df, 1))   &0x44912449;
2361
   DFWORD(result, 2)=(~DFWORD(df, 2))   &0x12449124;
2362
   DFWORD(result, 3)=(~DFWORD(df, 3))   &0x49124491;
2363
  #endif
2364
  return result;
2365
  } /* decFloatInvert */
2366
 
2367
/* ------------------------------------------------------------------ */
2368
/* decFloatIs -- decFloat tests (IsSigned, etc.)                      */
2369
/*                                                                    */
2370
/*   df is the decFloat to test                                       */
2371
/*   returns 0 or 1 in a uInt                                         */
2372
/*                                                                    */
2373
/* Many of these could be macros, but having them as real functions   */
2374
/* is a little cleaner (and they can be referred to here by the       */
2375
/* generic names)                                                     */
2376
/* ------------------------------------------------------------------ */
2377
uInt decFloatIsCanonical(const decFloat *df) {
2378
  if (DFISSPECIAL(df)) {
2379
    if (DFISINF(df)) {
2380
      if (DFWORD(df, 0)&ECONMASK) return 0;  /* exponent continuation */
2381
      if (!DFISCCZERO(df)) return 0;          /* coefficient continuation */
2382
      return 1;
2383
      }
2384
    /* is a NaN */
2385
    if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
2386
    if (DFISCCZERO(df)) return 1;            /* coefficient continuation */
2387
    /* drop through to check payload */
2388
    }
2389
  { /* declare block */
2390
  #if DOUBLE
2391
    uInt sourhi=DFWORD(df, 0);
2392
    uInt sourlo=DFWORD(df, 1);
2393
    if (CANONDPDOFF(sourhi, 8)
2394
     && CANONDPDTWO(sourhi, sourlo, 30)
2395
     && CANONDPDOFF(sourlo, 20)
2396
     && CANONDPDOFF(sourlo, 10)
2397
     && CANONDPDOFF(sourlo, 0)) return 1;
2398
  #elif QUAD
2399
    uInt sourhi=DFWORD(df, 0);
2400
    uInt sourmh=DFWORD(df, 1);
2401
    uInt sourml=DFWORD(df, 2);
2402
    uInt sourlo=DFWORD(df, 3);
2403
    if (CANONDPDOFF(sourhi, 4)
2404
     && CANONDPDTWO(sourhi, sourmh, 26)
2405
     && CANONDPDOFF(sourmh, 16)
2406
     && CANONDPDOFF(sourmh, 6)
2407
     && CANONDPDTWO(sourmh, sourml, 28)
2408
     && CANONDPDOFF(sourml, 18)
2409
     && CANONDPDOFF(sourml, 8)
2410
     && CANONDPDTWO(sourml, sourlo, 30)
2411
     && CANONDPDOFF(sourlo, 20)
2412
     && CANONDPDOFF(sourlo, 10)
2413
     && CANONDPDOFF(sourlo, 0)) return 1;
2414
  #endif
2415
  } /* block */
2416
  return 0;    /* a declet is non-canonical */
2417
  }
2418
 
2419
uInt decFloatIsFinite(const decFloat *df) {
2420
  return !DFISSPECIAL(df);
2421
  }
2422
uInt decFloatIsInfinite(const decFloat *df) {
2423
  return DFISINF(df);
2424
  }
2425
uInt decFloatIsInteger(const decFloat *df) {
2426
  return DFISINT(df);
2427
  }
2428
uInt decFloatIsNaN(const decFloat *df) {
2429
  return DFISNAN(df);
2430
  }
2431
uInt decFloatIsNormal(const decFloat *df) {
2432
  Int exp;                         /* exponent */
2433
  if (DFISSPECIAL(df)) return 0;
2434
  if (DFISZERO(df)) return 0;
2435
  /* is finite and non-zero */
2436
  exp=GETEXPUN(df)                 /* get unbiased exponent .. */
2437
     +decFloatDigits(df)-1;        /* .. and make adjusted exponent */
2438
  return (exp>=DECEMIN);           /* < DECEMIN is subnormal */
2439
  }
2440
uInt decFloatIsSignaling(const decFloat *df) {
2441
  return DFISSNAN(df);
2442
  }
2443
uInt decFloatIsSignalling(const decFloat *df) {
2444
  return DFISSNAN(df);
2445
  }
2446
uInt decFloatIsSigned(const decFloat *df) {
2447
  return DFISSIGNED(df);
2448
  }
2449
uInt decFloatIsSubnormal(const decFloat *df) {
2450
  if (DFISSPECIAL(df)) return 0;
2451
  /* is finite */
2452
  if (decFloatIsNormal(df)) return 0;
2453
  /* it is <Nmin, but could be zero */
2454
  if (DFISZERO(df)) return 0;
2455
  return 1;                                  /* is subnormal */
2456
  }
2457
uInt decFloatIsZero(const decFloat *df) {
2458
  return DFISZERO(df);
2459
  } /* decFloatIs... */
2460
 
2461
/* ------------------------------------------------------------------ */
2462
/* decFloatLogB -- return adjusted exponent, by 754 rules             */
2463
/*                                                                    */
2464
/*   result gets the adjusted exponent as an integer, or a NaN etc.   */
2465
/*   df     is the decFloat to be examined                            */
2466
/*   set    is the context                                            */
2467
/*   returns result                                                   */
2468
/*                                                                    */
2469
/* Notable cases:                                                     */
2470
/*   A<0 -> Use |A|                                                   */
2471
/*   A=0 -> -Infinity (Division by zero)                              */
2472
/*   A=Infinite -> +Infinity (Exact)                                  */
2473
/*   A=1 exactly -> 0 (Exact)                                         */
2474
/*   NaNs are propagated as usual                                     */
2475
/* ------------------------------------------------------------------ */
2476
decFloat * decFloatLogB(decFloat *result, const decFloat *df,
2477
                        decContext *set) {
2478
  Int ae;                                    /* adjusted exponent */
2479
  if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2480
  if (DFISINF(df)) {
2481
    DFWORD(result, 0)=0;               /* need +ve */
2482
    return decInfinity(result, result);      /* canonical +Infinity */
2483
    }
2484
  if (DFISZERO(df)) {
2485
    set->status|=DEC_Division_by_zero;       /* as per 754 */
2486
    DFWORD(result, 0)=DECFLOAT_Sign;          /* make negative */
2487
    return decInfinity(result, result);      /* canonical -Infinity */
2488
    }
2489
  ae=GETEXPUN(df)                       /* get unbiased exponent .. */
2490
    +decFloatDigits(df)-1;              /* .. and make adjusted exponent */
2491
  /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
2492
  /* it is worth using a special case of decFloatFromInt32 */
2493
  DFWORD(result, 0)=ZEROWORD;            /* always */
2494
  if (ae<0) {
2495
    DFWORD(result, 0)|=DECFLOAT_Sign;    /* -0 so far */
2496
    ae=-ae;
2497
    }
2498
  #if DOUBLE
2499
    DFWORD(result, 1)=BIN2DPD[ae];      /* a single declet */
2500
  #elif QUAD
2501
    DFWORD(result, 1)=0;
2502
    DFWORD(result, 2)=0;
2503
    DFWORD(result, 3)=(ae/1000)<<10;    /* is <10, so need no DPD encode */
2504
    DFWORD(result, 3)|=BIN2DPD[ae%1000];
2505
  #endif
2506
  return result;
2507
  } /* decFloatLogB */
2508
 
2509
/* ------------------------------------------------------------------ */
2510
/* decFloatMax -- return maxnum of two operands                       */
2511
/*                                                                    */
2512
/*   result gets the chosen decFloat                                  */
2513
/*   dfl    is the first decFloat (lhs)                               */
2514
/*   dfr    is the second decFloat (rhs)                              */
2515
/*   set    is the context                                            */
2516
/*   returns result                                                   */
2517
/*                                                                    */
2518
/* If just one operand is a quiet NaN it is ignored.                  */
2519
/* ------------------------------------------------------------------ */
2520
decFloat * decFloatMax(decFloat *result,
2521
                       const decFloat *dfl, const decFloat *dfr,
2522
                       decContext *set) {
2523
  Int comp;
2524
  if (DFISNAN(dfl)) {
2525
    /* sNaN or both NaNs leads to normal NaN processing */
2526
    if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2527
    return decCanonical(result, dfr);        /* RHS is numeric */
2528
    }
2529
  if (DFISNAN(dfr)) {
2530
    /* sNaN leads to normal NaN processing (both NaN handled above) */
2531
    if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2532
    return decCanonical(result, dfl);        /* LHS is numeric */
2533
    }
2534
  /* Both operands are numeric; numeric comparison needed -- use */
2535
  /* total order for a well-defined choice (and +0 > -0) */
2536
  comp=decNumCompare(dfl, dfr, 1);
2537
  if (comp>=0) return decCanonical(result, dfl);
2538
  return decCanonical(result, dfr);
2539
  } /* decFloatMax */
2540
 
2541
/* ------------------------------------------------------------------ */
2542
/* decFloatMaxMag -- return maxnummag of two operands                 */
2543
/*                                                                    */
2544
/*   result gets the chosen decFloat                                  */
2545
/*   dfl    is the first decFloat (lhs)                               */
2546
/*   dfr    is the second decFloat (rhs)                              */
2547
/*   set    is the context                                            */
2548
/*   returns result                                                   */
2549
/*                                                                    */
2550
/* Returns according to the magnitude comparisons if both numeric and */
2551
/* unequal, otherwise returns maxnum                                  */
2552
/* ------------------------------------------------------------------ */
2553
decFloat * decFloatMaxMag(decFloat *result,
2554
                       const decFloat *dfl, const decFloat *dfr,
2555
                       decContext *set) {
2556
  Int comp;
2557
  decFloat absl, absr;
2558
  if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
2559
 
2560
  decFloatCopyAbs(&absl, dfl);
2561
  decFloatCopyAbs(&absr, dfr);
2562
  comp=decNumCompare(&absl, &absr, 0);
2563
  if (comp>0) return decCanonical(result, dfl);
2564
  if (comp<0) return decCanonical(result, dfr);
2565
  return decFloatMax(result, dfl, dfr, set);
2566
  } /* decFloatMaxMag */
2567
 
2568
/* ------------------------------------------------------------------ */
2569
/* decFloatMin -- return minnum of two operands                       */
2570
/*                                                                    */
2571
/*   result gets the chosen decFloat                                  */
2572
/*   dfl    is the first decFloat (lhs)                               */
2573
/*   dfr    is the second decFloat (rhs)                              */
2574
/*   set    is the context                                            */
2575
/*   returns result                                                   */
2576
/*                                                                    */
2577
/* If just one operand is a quiet NaN it is ignored.                  */
2578
/* ------------------------------------------------------------------ */
2579
decFloat * decFloatMin(decFloat *result,
2580
                       const decFloat *dfl, const decFloat *dfr,
2581
                       decContext *set) {
2582
  Int comp;
2583
  if (DFISNAN(dfl)) {
2584
    /* sNaN or both NaNs leads to normal NaN processing */
2585
    if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
2586
    return decCanonical(result, dfr);        /* RHS is numeric */
2587
    }
2588
  if (DFISNAN(dfr)) {
2589
    /* sNaN leads to normal NaN processing (both NaN handled above) */
2590
    if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2591
    return decCanonical(result, dfl);        /* LHS is numeric */
2592
    }
2593
  /* Both operands are numeric; numeric comparison needed -- use */
2594
  /* total order for a well-defined choice (and +0 > -0) */
2595
  comp=decNumCompare(dfl, dfr, 1);
2596
  if (comp<=0) return decCanonical(result, dfl);
2597
  return decCanonical(result, dfr);
2598
  } /* decFloatMin */
2599
 
2600
/* ------------------------------------------------------------------ */
2601
/* decFloatMinMag -- return minnummag of two operands                 */
2602
/*                                                                    */
2603
/*   result gets the chosen decFloat                                  */
2604
/*   dfl    is the first decFloat (lhs)                               */
2605
/*   dfr    is the second decFloat (rhs)                              */
2606
/*   set    is the context                                            */
2607
/*   returns result                                                   */
2608
/*                                                                    */
2609
/* Returns according to the magnitude comparisons if both numeric and */
2610
/* unequal, otherwise returns minnum                                  */
2611
/* ------------------------------------------------------------------ */
2612
decFloat * decFloatMinMag(decFloat *result,
2613
                       const decFloat *dfl, const decFloat *dfr,
2614
                       decContext *set) {
2615
  Int comp;
2616
  decFloat absl, absr;
2617
  if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
2618
 
2619
  decFloatCopyAbs(&absl, dfl);
2620
  decFloatCopyAbs(&absr, dfr);
2621
  comp=decNumCompare(&absl, &absr, 0);
2622
  if (comp<0) return decCanonical(result, dfl);
2623
  if (comp>0) return decCanonical(result, dfr);
2624
  return decFloatMin(result, dfl, dfr, set);
2625
  } /* decFloatMinMag */
2626
 
2627
/* ------------------------------------------------------------------ */
2628
/* decFloatMinus -- negate value, heeding NaNs, etc.                  */
2629
/*                                                                    */
2630
/*   result gets the canonicalized 0-df                               */
2631
/*   df     is the decFloat to minus                                  */
2632
/*   set    is the context                                            */
2633
/*   returns result                                                   */
2634
/*                                                                    */
2635
/* This has the same effect as 0-df where the exponent of the zero is */
2636
/* the same as that of df (if df is finite).                          */
2637
/* The effect is also the same as decFloatCopyNegate except that NaNs */
2638
/* are handled normally (the sign of a NaN is not affected, and an    */
2639
/* sNaN will signal), the result is canonical, and zero gets sign 0.  */
2640
/* ------------------------------------------------------------------ */
2641
decFloat * decFloatMinus(decFloat *result, const decFloat *df,
2642
                         decContext *set) {
2643
  if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2644
  decCanonical(result, df);                       /* copy and check */
2645
  if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80;      /* turn off sign bit */
2646
   else DFBYTE(result, 0)^=0x80;           /* flip sign bit */
2647
  return result;
2648
  } /* decFloatMinus */
2649
 
2650
/* ------------------------------------------------------------------ */
2651
/* decFloatMultiply -- multiply two decFloats                         */
2652
/*                                                                    */
2653
/*   result gets the result of multiplying dfl and dfr:               */
2654
/*   dfl    is the first decFloat (lhs)                               */
2655
/*   dfr    is the second decFloat (rhs)                              */
2656
/*   set    is the context                                            */
2657
/*   returns result                                                   */
2658
/*                                                                    */
2659
/* ------------------------------------------------------------------ */
2660
decFloat * decFloatMultiply(decFloat *result,
2661
                            const decFloat *dfl, const decFloat *dfr,
2662
                            decContext *set) {
2663
  bcdnum num;                      /* for final conversion */
2664
  uByte  bcdacc[DECPMAX9*18+1];    /* for coefficent in BCD */
2665
 
2666
  if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
2667
    /* NaNs are handled as usual */
2668
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2669
    /* infinity times zero is bad */
2670
    if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
2671
    if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
2672
    /* both infinite; return canonical infinity with computed sign */
2673
    DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
2674
    return decInfinity(result, result);
2675
    }
2676
 
2677
  /* Here when both operands are finite */
2678
  decFiniteMultiply(&num, bcdacc, dfl, dfr);
2679
  return decFinalize(result, &num, set); /* round, check, and lay out */
2680
  } /* decFloatMultiply */
2681
 
2682
/* ------------------------------------------------------------------ */
2683
/* decFloatNextMinus -- next towards -Infinity                        */
2684
/*                                                                    */
2685
/*   result gets the next lesser decFloat                             */
2686
/*   dfl    is the decFloat to start with                             */
2687
/*   set    is the context                                            */
2688
/*   returns result                                                   */
2689
/*                                                                    */
2690
/* This is 754 nextdown; Invalid is the only status possible (from    */
2691
/* an sNaN).                                                          */
2692
/* ------------------------------------------------------------------ */
2693
decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
2694
                             decContext *set) {
2695
  decFloat delta;                       /* tiny increment */
2696
  uInt savestat;                        /* saves status */
2697
  enum rounding saveround;              /* .. and mode */
2698
 
2699
  /* +Infinity is the special case */
2700
  if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
2701
    DFSETNMAX(result);
2702
    return result;                      /* [no status to set] */
2703
    }
2704
  /* other cases are effected by sutracting a tiny delta -- this */
2705
  /* should be done in a wider format as the delta is unrepresentable */
2706
  /* here (but can be done with normal add if the sign of zero is */
2707
  /* treated carefully, because no Inexactitude is interesting); */
2708
  /* rounding to -Infinity then pushes the result to next below */
2709
  decFloatZero(&delta);                 /* set up tiny delta */
2710
  DFWORD(&delta, DECWORDS-1)=1;         /* coefficient=1 */
2711
  DFWORD(&delta, 0)=DECFLOAT_Sign;       /* Sign=1 + biased exponent=0 */
2712
  /* set up for the directional round */
2713
  saveround=set->round;                 /* save mode */
2714
  set->round=DEC_ROUND_FLOOR;           /* .. round towards -Infinity */
2715
  savestat=set->status;                 /* save status */
2716
  decFloatAdd(result, dfl, &delta, set);
2717
  /* Add rules mess up the sign when going from +Ntiny to 0 */
2718
  if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2719
  set->status&=DEC_Invalid_operation;   /* preserve only sNaN status */
2720
  set->status|=savestat;                /* restore pending flags */
2721
  set->round=saveround;                 /* .. and mode */
2722
  return result;
2723
  } /* decFloatNextMinus */
2724
 
2725
/* ------------------------------------------------------------------ */
2726
/* decFloatNextPlus -- next towards +Infinity                         */
2727
/*                                                                    */
2728
/*   result gets the next larger decFloat                             */
2729
/*   dfl    is the decFloat to start with                             */
2730
/*   set    is the context                                            */
2731
/*   returns result                                                   */
2732
/*                                                                    */
2733
/* This is 754 nextup; Invalid is the only status possible (from      */
2734
/* an sNaN).                                                          */
2735
/* ------------------------------------------------------------------ */
2736
decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
2737
                            decContext *set) {
2738
  uInt savestat;                        /* saves status */
2739
  enum rounding saveround;              /* .. and mode */
2740
  decFloat delta;                       /* tiny increment */
2741
 
2742
  /* -Infinity is the special case */
2743
  if (DFISINF(dfl) && DFISSIGNED(dfl)) {
2744
    DFSETNMAX(result);
2745
    DFWORD(result, 0)|=DECFLOAT_Sign;    /* make negative */
2746
    return result;                      /* [no status to set] */
2747
    }
2748
  /* other cases are effected by sutracting a tiny delta -- this */
2749
  /* should be done in a wider format as the delta is unrepresentable */
2750
  /* here (but can be done with normal add if the sign of zero is */
2751
  /* treated carefully, because no Inexactitude is interesting); */
2752
  /* rounding to +Infinity then pushes the result to next above */
2753
  decFloatZero(&delta);                 /* set up tiny delta */
2754
  DFWORD(&delta, DECWORDS-1)=1;         /* coefficient=1 */
2755
  DFWORD(&delta, 0)=0;                    /* Sign=0 + biased exponent=0 */
2756
  /* set up for the directional round */
2757
  saveround=set->round;                 /* save mode */
2758
  set->round=DEC_ROUND_CEILING;         /* .. round towards +Infinity */
2759
  savestat=set->status;                 /* save status */
2760
  decFloatAdd(result, dfl, &delta, set);
2761
  /* Add rules mess up the sign when going from -Ntiny to -0 */
2762
  if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
2763
  set->status&=DEC_Invalid_operation;   /* preserve only sNaN status */
2764
  set->status|=savestat;                /* restore pending flags */
2765
  set->round=saveround;                 /* .. and mode */
2766
  return result;
2767
  } /* decFloatNextPlus */
2768
 
2769
/* ------------------------------------------------------------------ */
2770
/* decFloatNextToward -- next towards a decFloat                      */
2771
/*                                                                    */
2772
/*   result gets the next decFloat                                    */
2773
/*   dfl    is the decFloat to start with                             */
2774
/*   dfr    is the decFloat to move toward                            */
2775
/*   set    is the context                                            */
2776
/*   returns result                                                   */
2777
/*                                                                    */
2778
/* This is 754-1985 nextafter, as modified during revision (dropped   */
2779
/* from 754-2008); status may be set unless the result is a normal    */
2780
/* number.                                                            */
2781
/* ------------------------------------------------------------------ */
2782
decFloat * decFloatNextToward(decFloat *result,
2783
                              const decFloat *dfl, const decFloat *dfr,
2784
                              decContext *set) {
2785
  decFloat delta;                       /* tiny increment or decrement */
2786
  decFloat pointone;                    /* 1e-1 */
2787
  uInt  savestat;                       /* saves status */
2788
  enum  rounding saveround;             /* .. and mode */
2789
  uInt  deltatop;                       /* top word for delta */
2790
  Int   comp;                           /* work */
2791
 
2792
  if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2793
  /* Both are numeric, so Invalid no longer a possibility */
2794
  comp=decNumCompare(dfl, dfr, 0);
2795
  if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
2796
  /* unequal; do NextPlus or NextMinus but with different status rules */
2797
 
2798
  if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
2799
    if (DFISINF(dfl) && DFISSIGNED(dfl)) {   /* -Infinity special case */
2800
      DFSETNMAX(result);
2801
      DFWORD(result, 0)|=DECFLOAT_Sign;
2802
      return result;
2803
      }
2804
    saveround=set->round;                    /* save mode */
2805
    set->round=DEC_ROUND_CEILING;            /* .. round towards +Infinity */
2806
    deltatop=0;                       /* positive delta */
2807
    }
2808
   else { /* lhs>rhs, do NextMinus, see above for commentary */
2809
    if (DFISINF(dfl) && !DFISSIGNED(dfl)) {  /* +Infinity special case */
2810
      DFSETNMAX(result);
2811
      return result;
2812
      }
2813
    saveround=set->round;                    /* save mode */
2814
    set->round=DEC_ROUND_FLOOR;              /* .. round towards -Infinity */
2815
    deltatop=DECFLOAT_Sign;                  /* negative delta */
2816
    }
2817
  savestat=set->status;                      /* save status */
2818
  /* Here, Inexact is needed where appropriate (and hence Underflow, */
2819
  /* etc.).  Therefore the tiny delta which is otherwise */
2820
  /* unrepresentable (see NextPlus and NextMinus) is constructed */
2821
  /* using the multiplication of FMA. */
2822
  decFloatZero(&delta);                 /* set up tiny delta */
2823
  DFWORD(&delta, DECWORDS-1)=1;         /* coefficient=1 */
2824
  DFWORD(&delta, 0)=deltatop;            /* Sign + biased exponent=0 */
2825
  decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
2826
  decFloatFMA(result, &delta, &pointone, dfl, set);
2827
  /* [Delta is truly tiny, so no need to correct sign of zero] */
2828
  /* use new status unless the result is normal */
2829
  if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
2830
  set->round=saveround;                 /* restore mode */
2831
  return result;
2832
  } /* decFloatNextToward */
2833
 
2834
/* ------------------------------------------------------------------ */
2835
/* decFloatOr -- logical digitwise OR of two decFloats                */
2836
/*                                                                    */
2837
/*   result gets the result of ORing dfl and dfr                      */
2838
/*   dfl    is the first decFloat (lhs)                               */
2839
/*   dfr    is the second decFloat (rhs)                              */
2840
/*   set    is the context                                            */
2841
/*   returns result, which will be canonical with sign=0              */
2842
/*                                                                    */
2843
/* The operands must be positive, finite with exponent q=0, and       */
2844
/* comprise just zeros and ones; if not, Invalid operation results.   */
2845
/* ------------------------------------------------------------------ */
2846
decFloat * decFloatOr(decFloat *result,
2847
                       const decFloat *dfl, const decFloat *dfr,
2848
                       decContext *set) {
2849
  if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
2850
   || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
2851
  /* the operands are positive finite integers (q=0) with just 0s and 1s */
2852
  #if DOUBLE
2853
   DFWORD(result, 0)=ZEROWORD
2854
                   |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
2855
   DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
2856
  #elif QUAD
2857
   DFWORD(result, 0)=ZEROWORD
2858
                   |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
2859
   DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
2860
   DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
2861
   DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
2862
  #endif
2863
  return result;
2864
  } /* decFloatOr */
2865
 
2866
/* ------------------------------------------------------------------ */
2867
/* decFloatPlus -- add value to 0, heeding NaNs, etc.                 */
2868
/*                                                                    */
2869
/*   result gets the canonicalized 0+df                               */
2870
/*   df     is the decFloat to plus                                   */
2871
/*   set    is the context                                            */
2872
/*   returns result                                                   */
2873
/*                                                                    */
2874
/* This has the same effect as 0+df where the exponent of the zero is */
2875
/* the same as that of df (if df is finite).                          */
2876
/* The effect is also the same as decFloatCopy except that NaNs       */
2877
/* are handled normally (the sign of a NaN is not affected, and an    */
2878
/* sNaN will signal), the result is canonical, and zero gets sign 0.  */
2879
/* ------------------------------------------------------------------ */
2880
decFloat * decFloatPlus(decFloat *result, const decFloat *df,
2881
                        decContext *set) {
2882
  if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
2883
  decCanonical(result, df);                       /* copy and check */
2884
  if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80;      /* turn off sign bit */
2885
  return result;
2886
  } /* decFloatPlus */
2887
 
2888
/* ------------------------------------------------------------------ */
2889
/* decFloatQuantize -- quantize a decFloat                            */
2890
/*                                                                    */
2891
/*   result gets the result of quantizing dfl to match dfr            */
2892
/*   dfl    is the first decFloat (lhs)                               */
2893
/*   dfr    is the second decFloat (rhs), which sets the exponent     */
2894
/*   set    is the context                                            */
2895
/*   returns result                                                   */
2896
/*                                                                    */
2897
/* Unless there is an error or the result is infinite, the exponent   */
2898
/* of result is guaranteed to be the same as that of dfr.             */
2899
/* ------------------------------------------------------------------ */
2900
decFloat * decFloatQuantize(decFloat *result,
2901
                            const decFloat *dfl, const decFloat *dfr,
2902
                            decContext *set) {
2903
  Int   explb, exprb;         /* left and right biased exponents */
2904
  uByte *ulsd;                /* local LSD pointer */
2905
  uByte *ub, *uc;             /* work */
2906
  Int   drop;                 /* .. */
2907
  uInt  dpd;                  /* .. */
2908
  uInt  encode;               /* encoding accumulator */
2909
  uInt  sourhil, sourhir;     /* top words from source decFloats */
2910
  uInt  uiwork;               /* for macros */
2911
  #if QUAD
2912
  uShort uswork;              /* .. */
2913
  #endif
2914
  /* the following buffer holds the coefficient for manipulation */
2915
  uByte buf[4+DECPMAX*3+2*QUAD];   /* + space for zeros to left or right */
2916
  #if DECTRACE
2917
  bcdnum num;                      /* for trace displays */
2918
  #endif
2919
 
2920
  /* Start decoding the arguments */
2921
  sourhil=DFWORD(dfl, 0);           /* LHS top word */
2922
  explb=DECCOMBEXP[sourhil>>26];   /* get exponent high bits (in place) */
2923
  sourhir=DFWORD(dfr, 0);           /* RHS top word */
2924
  exprb=DECCOMBEXP[sourhir>>26];
2925
 
2926
  if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
2927
    /* NaNs are handled as usual */
2928
    if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
2929
    /* one infinity but not both is bad */
2930
    if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
2931
    /* both infinite; return canonical infinity with sign of LHS */
2932
    return decInfinity(result, dfl);
2933
    }
2934
 
2935
  /* Here when both arguments are finite */
2936
  /* complete extraction of the exponents [no need to unbias] */
2937
  explb+=GETECON(dfl);             /* + continuation */
2938
  exprb+=GETECON(dfr);             /* .. */
2939
 
2940
  /* calculate the number of digits to drop from the coefficient */
2941
  drop=exprb-explb;                /* 0 if nothing to do */
2942
  if (drop==0) return decCanonical(result, dfl); /* return canonical */
2943
 
2944
  /* the coefficient is needed; lay it out into buf, offset so zeros */
2945
  /* can be added before or after as needed -- an extra heading is */
2946
  /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
2947
  /* fours */
2948
  #define BUFOFF (buf+4+DECPMAX)
2949
  GETCOEFF(dfl, BUFOFF);           /* decode from decFloat */
2950
  /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
2951
 
2952
  #if DECTRACE
2953
  num.msd=BUFOFF;
2954
  num.lsd=BUFOFF+DECPMAX-1;
2955
  num.exponent=explb-DECBIAS;
2956
  num.sign=sourhil & DECFLOAT_Sign;
2957
  decShowNum(&num, "dfl");
2958
  #endif
2959
 
2960
  if (drop>0) {                  /* [most common case] */
2961
    /* (this code is very similar to that in decFloatFinalize, but */
2962
    /* has many differences so is duplicated here -- so any changes */
2963
    /* may need to be made there, too) */
2964
    uByte *roundat;                          /* -> re-round digit */
2965
    uByte reround;                           /* reround value */
2966
    /* printf("Rounding; drop=%ld\n", (LI)drop); */
2967
 
2968
    /* there is at least one zero needed to the left, in all but one */
2969
    /* exceptional (all-nines) case, so place four zeros now; this is */
2970
    /* needed almost always and makes rounding all-nines by fours safe */
2971
    UBFROMUI(BUFOFF-4, 0);
2972
 
2973
    /* Three cases here: */
2974
    /*   1. new LSD is in coefficient (almost always) */
2975
    /*   2. new LSD is digit to left of coefficient (so MSD is */
2976
    /*      round-for-reround digit) */
2977
    /*   3. new LSD is to left of case 2 (whole coefficient is sticky) */
2978
    /* Note that leading zeros can safely be treated as useful digits */
2979
 
2980
    /* [duplicate check-stickies code to save a test] */
2981
    /* [by-digit check for stickies as runs of zeros are rare] */
2982
    if (drop<DECPMAX) {                      /* NB lengths not addresses */
2983
      roundat=BUFOFF+DECPMAX-drop;
2984
      reround=*roundat;
2985
      for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
2986
        if (*ub!=0) {                         /* non-zero to be discarded */
2987
          reround=DECSTICKYTAB[reround];     /* apply sticky bit */
2988
          break;                             /* [remainder don't-care] */
2989
          }
2990
        } /* check stickies */
2991
      ulsd=roundat-1;                        /* set LSD */
2992
      }
2993
     else {                                  /* edge case */
2994
      if (drop==DECPMAX) {
2995
        roundat=BUFOFF;
2996
        reround=*roundat;
2997
        }
2998
       else {
2999
        roundat=BUFOFF-1;
3000
        reround=0;
3001
        }
3002
      for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
3003
        if (*ub!=0) {                         /* non-zero to be discarded */
3004
          reround=DECSTICKYTAB[reround];     /* apply sticky bit */
3005
          break;                             /* [remainder don't-care] */
3006
          }
3007
        } /* check stickies */
3008
      *BUFOFF=0;                      /* make a coefficient of 0 */
3009
      ulsd=BUFOFF;                           /* .. at the MSD place */
3010
      }
3011
 
3012
    if (reround!=0) {                         /* discarding non-zero */
3013
      uInt bump=0;
3014
      set->status|=DEC_Inexact;
3015
 
3016
      /* next decide whether to increment the coefficient */
3017
      if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
3018
        if (reround>5) bump=1;               /* >0.5 goes up */
3019
         else if (reround==5)                /* exactly 0.5000 .. */
3020
          bump=*ulsd & 0x01;                 /* .. up iff [new] lsd is odd */
3021
        } /* r-h-e */
3022
       else switch (set->round) {
3023
        case DEC_ROUND_DOWN: {
3024
          /* no change */
3025
          break;} /* r-d */
3026
        case DEC_ROUND_HALF_DOWN: {
3027
          if (reround>5) bump=1;
3028
          break;} /* r-h-d */
3029
        case DEC_ROUND_HALF_UP: {
3030
          if (reround>=5) bump=1;
3031
          break;} /* r-h-u */
3032
        case DEC_ROUND_UP: {
3033
          if (reround>0) bump=1;
3034
          break;} /* r-u */
3035
        case DEC_ROUND_CEILING: {
3036
          /* same as _UP for positive numbers, and as _DOWN for negatives */
3037
          if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
3038
          break;} /* r-c */
3039
        case DEC_ROUND_FLOOR: {
3040
          /* same as _UP for negative numbers, and as _DOWN for positive */
3041
          /* [negative reround cannot occur on 0] */
3042
          if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
3043
          break;} /* r-f */
3044
        case DEC_ROUND_05UP: {
3045
          if (reround>0) { /* anything out there is 'sticky' */
3046
            /* bump iff lsd=0 or 5; this cannot carry so it could be */
3047
            /* effected immediately with no bump -- but the code */
3048
            /* is clearer if this is done the same way as the others */
3049
            if (*ulsd==0 || *ulsd==5) bump=1;
3050
            }
3051
          break;} /* r-r */
3052
        default: {      /* e.g., DEC_ROUND_MAX */
3053
          set->status|=DEC_Invalid_context;
3054
          #if DECCHECK
3055
          printf("Unknown rounding mode: %ld\n", (LI)set->round);
3056
          #endif
3057
          break;}
3058
        } /* switch (not r-h-e) */
3059
      /* printf("ReRound: %ld  bump: %ld\n", (LI)reround, (LI)bump); */
3060
 
3061
      if (bump!=0) {                          /* need increment */
3062
        /* increment the coefficient; this could give 1000... (after */
3063
        /* the all nines case) */
3064
        ub=ulsd;
3065
        for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
3066
        /* now at most 3 digits left to non-9 (usually just the one) */
3067
        for (; *ub==9; ub--) *ub=0;
3068
        *ub+=1;
3069
        /* [the all-nines case will have carried one digit to the */
3070
        /* left of the original MSD -- just where it is needed] */
3071
        } /* bump needed */
3072
      } /* inexact rounding */
3073
 
3074
    /* now clear zeros to the left so exactly DECPMAX digits will be */
3075
    /* available in the coefficent -- the first word to the left was */
3076
    /* cleared earlier for safe carry; now add any more needed */
3077
    if (drop>4) {
3078
      UBFROMUI(BUFOFF-8, 0);                  /* must be at least 5 */
3079
      for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 0);
3080
      }
3081
    } /* need round (drop>0) */
3082
 
3083
   else { /* drop<0; padding with -drop digits is needed */
3084
    /* This is the case where an error can occur if the padded */
3085
    /* coefficient will not fit; checking for this can be done in the */
3086
    /* same loop as padding for zeros if the no-hope and zero cases */
3087
    /* are checked first */
3088
    if (-drop>DECPMAX-1) {                   /* cannot fit unless 0 */
3089
      if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
3090
      /* a zero can have any exponent; just drop through and use it */
3091
      ulsd=BUFOFF+DECPMAX-1;
3092
      }
3093
     else { /* padding will fit (but may still be too long) */
3094
      /* final-word mask depends on endianess */
3095
      #if DECLITEND
3096
      static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
3097
      #else
3098
      static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
3099
      #endif
3100
      /* note that here zeros to the right are added by fours, so in */
3101
      /* the Quad case this could write 36 zeros if the coefficient has */
3102
      /* fewer than three significant digits (hence the +2*QUAD for buf) */
3103
      for (uc=BUFOFF+DECPMAX;; uc+=4) {
3104
        UBFROMUI(uc, 0);
3105
        if (UBTOUI(uc-DECPMAX)!=0) {               /* could be bad */
3106
          /* if all four digits should be zero, definitely bad */
3107
          if (uc<=BUFOFF+DECPMAX+(-drop)-4)
3108
            return decInvalid(result, set);
3109
          /* must be a 1- to 3-digit sequence; check more carefully */
3110
          if ((UBTOUI(uc-DECPMAX)&dmask[(-drop)%4])!=0)
3111
            return decInvalid(result, set);
3112
          break;    /* no need for loop end test */
3113
          }
3114
        if (uc>=BUFOFF+DECPMAX+(-drop)-4) break;  /* done */
3115
        }
3116
      ulsd=BUFOFF+DECPMAX+(-drop)-1;
3117
      } /* pad and check leading zeros */
3118
    } /* drop<0 */
3119
 
3120
  #if DECTRACE
3121
  num.msd=ulsd-DECPMAX+1;
3122
  num.lsd=ulsd;
3123
  num.exponent=explb-DECBIAS;
3124
  num.sign=sourhil & DECFLOAT_Sign;
3125
  decShowNum(&num, "res");
3126
  #endif
3127
 
3128
  /*------------------------------------------------------------------*/
3129
  /* At this point the result is DECPMAX digits, ending at ulsd, so   */
3130
  /* fits the encoding exactly; there is no possibility of error      */
3131
  /*------------------------------------------------------------------*/
3132
  encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
3133
  encode=DECCOMBFROM[encode];                /* indexed by (0-2)*16+msd */
3134
  /* the exponent continuation can be extracted from the original RHS */
3135
  encode|=sourhir & ECONMASK;
3136
  encode|=sourhil&DECFLOAT_Sign;             /* add the sign from LHS */
3137
 
3138
  /* finally encode the coefficient */
3139
  /* private macro to encode a declet; this version can be used */
3140
  /* because all coefficient digits exist */
3141
  #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2;                   \
3142
    dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
3143
 
3144
  #if DOUBLE
3145
    getDPD3q(dpd, 4); encode|=dpd<<8;
3146
    getDPD3q(dpd, 3); encode|=dpd>>2;
3147
    DFWORD(result, 0)=encode;
3148
    encode=dpd<<30;
3149
    getDPD3q(dpd, 2); encode|=dpd<<20;
3150
    getDPD3q(dpd, 1); encode|=dpd<<10;
3151
    getDPD3q(dpd, 0); encode|=dpd;
3152
    DFWORD(result, 1)=encode;
3153
 
3154
  #elif QUAD
3155
    getDPD3q(dpd,10); encode|=dpd<<4;
3156
    getDPD3q(dpd, 9); encode|=dpd>>6;
3157
    DFWORD(result, 0)=encode;
3158
    encode=dpd<<26;
3159
    getDPD3q(dpd, 8); encode|=dpd<<16;
3160
    getDPD3q(dpd, 7); encode|=dpd<<6;
3161
    getDPD3q(dpd, 6); encode|=dpd>>4;
3162
    DFWORD(result, 1)=encode;
3163
    encode=dpd<<28;
3164
    getDPD3q(dpd, 5); encode|=dpd<<18;
3165
    getDPD3q(dpd, 4); encode|=dpd<<8;
3166
    getDPD3q(dpd, 3); encode|=dpd>>2;
3167
    DFWORD(result, 2)=encode;
3168
    encode=dpd<<30;
3169
    getDPD3q(dpd, 2); encode|=dpd<<20;
3170
    getDPD3q(dpd, 1); encode|=dpd<<10;
3171
    getDPD3q(dpd, 0); encode|=dpd;
3172
    DFWORD(result, 3)=encode;
3173
  #endif
3174
  return result;
3175
  } /* decFloatQuantize */
3176
 
3177
/* ------------------------------------------------------------------ */
3178
/* decFloatReduce -- reduce finite coefficient to minimum length      */
3179
/*                                                                    */
3180
/*   result gets the reduced decFloat                                 */
3181
/*   df     is the source decFloat                                    */
3182
/*   set    is the context                                            */
3183
/*   returns result, which will be canonical                          */
3184
/*                                                                    */
3185
/* This removes all possible trailing zeros from the coefficient;     */
3186
/* some may remain when the number is very close to Nmax.             */
3187
/* Special values are unchanged and no status is set unless df=sNaN.  */
3188
/* Reduced zero has an exponent q=0.                                  */
3189
/* ------------------------------------------------------------------ */
3190
decFloat * decFloatReduce(decFloat *result, const decFloat *df,
3191
                          decContext *set) {
3192
  bcdnum num;                           /* work */
3193
  uByte buf[DECPMAX], *ub;              /* coefficient and pointer */
3194
  if (df!=result) *result=*df;          /* copy, if needed */
3195
  if (DFISNAN(df)) return decNaNs(result, df, NULL, set);   /* sNaN */
3196
  /* zeros and infinites propagate too */
3197
  if (DFISINF(df)) return decInfinity(result, df);     /* canonical */
3198
  if (DFISZERO(df)) {
3199
    uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
3200
    decFloatZero(result);
3201
    DFWORD(result, 0)|=sign;
3202
    return result;                      /* exponent dropped, sign OK */
3203
    }
3204
  /* non-zero finite */
3205
  GETCOEFF(df, buf);
3206
  ub=buf+DECPMAX-1;                     /* -> lsd */
3207
  if (*ub) return result;               /* no trailing zeros */
3208
  for (ub--; *ub==0;) ub--;              /* terminates because non-zero */
3209
  /* *ub is the first non-zero from the right */
3210
  num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
3211
  num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
3212
  num.msd=buf;
3213
  num.lsd=ub;
3214
  return decFinalize(result, &num, set);
3215
  } /* decFloatReduce */
3216
 
3217
/* ------------------------------------------------------------------ */
3218
/* decFloatRemainder -- integer divide and return remainder           */
3219
/*                                                                    */
3220
/*   result gets the remainder of dividing dfl by dfr:                */
3221
/*   dfl    is the first decFloat (lhs)                               */
3222
/*   dfr    is the second decFloat (rhs)                              */
3223
/*   set    is the context                                            */
3224
/*   returns result                                                   */
3225
/*                                                                    */
3226
/* ------------------------------------------------------------------ */
3227
decFloat * decFloatRemainder(decFloat *result,
3228
                             const decFloat *dfl, const decFloat *dfr,
3229
                             decContext *set) {
3230
  return decDivide(result, dfl, dfr, set, REMAINDER);
3231
  } /* decFloatRemainder */
3232
 
3233
/* ------------------------------------------------------------------ */
3234
/* decFloatRemainderNear -- integer divide to nearest and remainder   */
3235
/*                                                                    */
3236
/*   result gets the remainder of dividing dfl by dfr:                */
3237
/*   dfl    is the first decFloat (lhs)                               */
3238
/*   dfr    is the second decFloat (rhs)                              */
3239
/*   set    is the context                                            */
3240
/*   returns result                                                   */
3241
/*                                                                    */
3242
/* This is the IEEE remainder, where the nearest integer is used.     */
3243
/* ------------------------------------------------------------------ */
3244
decFloat * decFloatRemainderNear(decFloat *result,
3245
                             const decFloat *dfl, const decFloat *dfr,
3246
                             decContext *set) {
3247
  return decDivide(result, dfl, dfr, set, REMNEAR);
3248
  } /* decFloatRemainderNear */
3249
 
3250
/* ------------------------------------------------------------------ */
3251
/* decFloatRotate -- rotate the coefficient of a decFloat left/right  */
3252
/*                                                                    */
3253
/*   result gets the result of rotating dfl                           */
3254
/*   dfl    is the source decFloat to rotate                          */
3255
/*   dfr    is the count of digits to rotate, an integer (with q=0)   */
3256
/*   set    is the context                                            */
3257
/*   returns result                                                   */
3258
/*                                                                    */
3259
/* The digits of the coefficient of dfl are rotated to the left (if   */
3260
/* dfr is positive) or to the right (if dfr is negative) without      */
3261
/* adjusting the exponent or the sign of dfl.                         */
3262
/*                                                                    */
3263
/* dfr must be in the range -DECPMAX through +DECPMAX.                */
3264
/* NaNs are propagated as usual.  An infinite dfl is unaffected (but  */
3265
/* dfr must be valid).  No status is set unless dfr is invalid or an  */
3266
/* operand is an sNaN.  The result is canonical.                      */
3267
/* ------------------------------------------------------------------ */
3268
#define PHALF (ROUNDUP(DECPMAX/2, 4))   /* half length, rounded up */
3269
decFloat * decFloatRotate(decFloat *result,
3270
                         const decFloat *dfl, const decFloat *dfr,
3271
                         decContext *set) {
3272
  Int rotate;                           /* dfr as an Int */
3273
  uByte buf[DECPMAX+PHALF];             /* coefficient + half */
3274
  uInt digits, savestat;                /* work */
3275
  bcdnum num;                           /* .. */
3276
  uByte *ub;                            /* .. */
3277
 
3278
  if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3279
  if (!DFISINT(dfr)) return decInvalid(result, set);
3280
  digits=decFloatDigits(dfr);                    /* calculate digits */
3281
  if (digits>2) return decInvalid(result, set);  /* definitely out of range */
3282
  rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
3283
  if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
3284
  /* [from here on no error or status change is possible] */
3285
  if (DFISINF(dfl)) return decInfinity(result, dfl);  /* canonical */
3286
  /* handle no-rotate cases */
3287
  if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
3288
  /* a real rotate is needed: 0 < rotate < DECPMAX */
3289
  /* reduce the rotation to no more than half to reduce copying later */
3290
  /* (for QUAD in fact half + 2 digits) */
3291
  if (DFISSIGNED(dfr)) rotate=-rotate;
3292
  if (abs(rotate)>PHALF) {
3293
    if (rotate<0) rotate=DECPMAX+rotate;
3294
     else rotate=rotate-DECPMAX;
3295
    }
3296
  /* now lay out the coefficient, leaving room to the right or the */
3297
  /* left depending on the direction of rotation */
3298
  ub=buf;
3299
  if (rotate<0) ub+=PHALF;    /* rotate right, so space to left */
3300
  GETCOEFF(dfl, ub);
3301
  /* copy half the digits to left or right, and set num.msd */
3302
  if (rotate<0) {
3303
    memcpy(buf, buf+DECPMAX, PHALF);
3304
    num.msd=buf+PHALF+rotate;
3305
    }
3306
   else {
3307
    memcpy(buf+DECPMAX, buf, PHALF);
3308
    num.msd=buf+rotate;
3309
    }
3310
  /* fill in rest of num */
3311
  num.lsd=num.msd+DECPMAX-1;
3312
  num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3313
  num.exponent=GETEXPUN(dfl);
3314
  savestat=set->status;                 /* record */
3315
  decFinalize(result, &num, set);
3316
  set->status=savestat;                 /* restore */
3317
  return result;
3318
  } /* decFloatRotate */
3319
 
3320
/* ------------------------------------------------------------------ */
3321
/* decFloatSameQuantum -- test decFloats for same quantum             */
3322
/*                                                                    */
3323
/*   dfl    is the first decFloat (lhs)                               */
3324
/*   dfr    is the second decFloat (rhs)                              */
3325
/*   returns 1 if the operands have the same quantum, 0 otherwise     */
3326
/*                                                                    */
3327
/* No error is possible and no status results.                        */
3328
/* ------------------------------------------------------------------ */
3329
uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
3330
  if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
3331
    if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
3332
    if (DFISINF(dfl) && DFISINF(dfr)) return 1;
3333
    return 0;  /* any other special mixture gives false */
3334
    }
3335
  if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
3336
  return 0;
3337
  } /* decFloatSameQuantum */
3338
 
3339
/* ------------------------------------------------------------------ */
3340
/* decFloatScaleB -- multiply by a power of 10, as per 754            */
3341
/*                                                                    */
3342
/*   result gets the result of the operation                          */
3343
/*   dfl    is the first decFloat (lhs)                               */
3344
/*   dfr    is the second decFloat (rhs), am integer (with q=0)       */
3345
/*   set    is the context                                            */
3346
/*   returns result                                                   */
3347
/*                                                                    */
3348
/* This computes result=dfl x 10**dfr where dfr is an integer in the  */
3349
/* range +/-2*(emax+pmax), typically resulting from LogB.             */
3350
/* Underflow and Overflow (with Inexact) may occur.  NaNs propagate   */
3351
/* as usual.                                                          */
3352
/* ------------------------------------------------------------------ */
3353
#define SCALEBMAX 2*(DECEMAX+DECPMAX)   /* D=800, Q=12356 */
3354
decFloat * decFloatScaleB(decFloat *result,
3355
                          const decFloat *dfl, const decFloat *dfr,
3356
                          decContext *set) {
3357
  uInt digits;                          /* work */
3358
  Int  expr;                            /* dfr as an Int */
3359
 
3360
  if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3361
  if (!DFISINT(dfr)) return decInvalid(result, set);
3362
  digits=decFloatDigits(dfr);                /* calculate digits */
3363
 
3364
  #if DOUBLE
3365
  if (digits>3) return decInvalid(result, set);   /* definitely out of range */
3366
  expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff];             /* must be in bottom declet */
3367
  #elif QUAD
3368
  if (digits>5) return decInvalid(result, set);   /* definitely out of range */
3369
  expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff]              /* in bottom 2 declets .. */
3370
      +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000;  /* .. */
3371
  #endif
3372
  if (expr>SCALEBMAX) return decInvalid(result, set);  /* oops */
3373
  /* [from now on no error possible] */
3374
  if (DFISINF(dfl)) return decInfinity(result, dfl);   /* canonical */
3375
  if (DFISSIGNED(dfr)) expr=-expr;
3376
  /* dfl is finite and expr is valid */
3377
  *result=*dfl;                              /* copy to target */
3378
  return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
3379
  } /* decFloatScaleB */
3380
 
3381
/* ------------------------------------------------------------------ */
3382
/* decFloatShift -- shift the coefficient of a decFloat left or right */
3383
/*                                                                    */
3384
/*   result gets the result of shifting dfl                           */
3385
/*   dfl    is the source decFloat to shift                           */
3386
/*   dfr    is the count of digits to shift, an integer (with q=0)    */
3387
/*   set    is the context                                            */
3388
/*   returns result                                                   */
3389
/*                                                                    */
3390
/* The digits of the coefficient of dfl are shifted to the left (if   */
3391
/* dfr is positive) or to the right (if dfr is negative) without      */
3392
/* adjusting the exponent or the sign of dfl.                         */
3393
/*                                                                    */
3394
/* dfr must be in the range -DECPMAX through +DECPMAX.                */
3395
/* NaNs are propagated as usual.  An infinite dfl is unaffected (but  */
3396
/* dfr must be valid).  No status is set unless dfr is invalid or an  */
3397
/* operand is an sNaN.  The result is canonical.                      */
3398
/* ------------------------------------------------------------------ */
3399
decFloat * decFloatShift(decFloat *result,
3400
                         const decFloat *dfl, const decFloat *dfr,
3401
                         decContext *set) {
3402
  Int    shift;                         /* dfr as an Int */
3403
  uByte  buf[DECPMAX*2];                /* coefficient + padding */
3404
  uInt   digits, savestat;              /* work */
3405
  bcdnum num;                           /* .. */
3406
  uInt   uiwork;                        /* for macros */
3407
 
3408
  if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
3409
  if (!DFISINT(dfr)) return decInvalid(result, set);
3410
  digits=decFloatDigits(dfr);                     /* calculate digits */
3411
  if (digits>2) return decInvalid(result, set);   /* definitely out of range */
3412
  shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff];   /* is in bottom declet */
3413
  if (shift>DECPMAX) return decInvalid(result, set);   /* too big */
3414
  /* [from here on no error or status change is possible] */
3415
 
3416
  if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
3417
  /* handle no-shift and all-shift (clear to zero) cases */
3418
  if (shift==0) return decCanonical(result, dfl);
3419
  if (shift==DECPMAX) {                      /* zero with sign */
3420
    uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
3421
    decFloatZero(result);                    /* make +0 */
3422
    DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
3423
    /* [cannot safely use CopySign] */
3424
    return result;
3425
    }
3426
  /* a real shift is needed: 0 < shift < DECPMAX */
3427
  num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
3428
  num.exponent=GETEXPUN(dfl);
3429
  num.msd=buf;
3430
  GETCOEFF(dfl, buf);
3431
  if (DFISSIGNED(dfr)) { /* shift right */
3432
    /* edge cases are taken care of, so this is easy */
3433
    num.lsd=buf+DECPMAX-shift-1;
3434
    }
3435
   else { /* shift left -- zero padding needed to right */
3436
    UBFROMUI(buf+DECPMAX, 0);            /* 8 will handle most cases */
3437
    UBFROMUI(buf+DECPMAX+4, 0);  /* .. */
3438
    if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
3439
    num.msd+=shift;
3440
    num.lsd=num.msd+DECPMAX-1;
3441
    }
3442
  savestat=set->status;                 /* record */
3443
  decFinalize(result, &num, set);
3444
  set->status=savestat;                 /* restore */
3445
  return result;
3446
  } /* decFloatShift */
3447
 
3448
/* ------------------------------------------------------------------ */
3449
/* decFloatSubtract -- subtract a decFloat from another               */
3450
/*                                                                    */
3451
/*   result gets the result of subtracting dfr from dfl:              */
3452
/*   dfl    is the first decFloat (lhs)                               */
3453
/*   dfr    is the second decFloat (rhs)                              */
3454
/*   set    is the context                                            */
3455
/*   returns result                                                   */
3456
/*                                                                    */
3457
/* ------------------------------------------------------------------ */
3458
decFloat * decFloatSubtract(decFloat *result,
3459
                            const decFloat *dfl, const decFloat *dfr,
3460
                            decContext *set) {
3461
  decFloat temp;
3462
  /* NaNs must propagate without sign change */
3463
  if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
3464
  temp=*dfr;                                   /* make a copy */
3465
  DFBYTE(&temp, 0)^=0x80;                       /* flip sign */
3466
  return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
3467
  } /* decFloatSubtract */
3468
 
3469
/* ------------------------------------------------------------------ */
3470
/* decFloatToInt -- round to 32-bit binary integer (4 flavours)       */
3471
/*                                                                    */
3472
/*   df    is the decFloat to round                                   */
3473
/*   set   is the context                                             */
3474
/*   round is the rounding mode to use                                */
3475
/*   returns a uInt or an Int, rounded according to the name          */
3476
/*                                                                    */
3477
/* Invalid will always be signaled if df is a NaN, is Infinite, or is */
3478
/* outside the range of the target; Inexact will not be signaled for  */
3479
/* simple rounding unless 'Exact' appears in the name.                */
3480
/* ------------------------------------------------------------------ */
3481
uInt decFloatToUInt32(const decFloat *df, decContext *set,
3482
                      enum rounding round) {
3483
  return decToInt32(df, set, round, 0, 1);}
3484
 
3485
uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
3486
                           enum rounding round) {
3487
  return decToInt32(df, set, round, 1, 1);}
3488
 
3489
Int decFloatToInt32(const decFloat *df, decContext *set,
3490
                    enum rounding round) {
3491
  return (Int)decToInt32(df, set, round, 0, 0);}
3492
 
3493
Int decFloatToInt32Exact(const decFloat *df, decContext *set,
3494
                         enum rounding round) {
3495
  return (Int)decToInt32(df, set, round, 1, 0);}
3496
 
3497
/* ------------------------------------------------------------------ */
3498
/* decFloatToIntegral -- round to integral value (two flavours)       */
3499
/*                                                                    */
3500
/*   result gets the result                                           */
3501
/*   df     is the decFloat to round                                  */
3502
/*   set    is the context                                            */
3503
/*   round  is the rounding mode to use                               */
3504
/*   returns result                                                   */
3505
/*                                                                    */
3506
/* No exceptions, even Inexact, are raised except for sNaN input, or  */
3507
/* if 'Exact' appears in the name.                                    */
3508
/* ------------------------------------------------------------------ */
3509
decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
3510
                                   decContext *set, enum rounding round) {
3511
  return decToIntegral(result, df, set, round, 0);}
3512
 
3513
decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
3514
                                   decContext *set) {
3515
  return decToIntegral(result, df, set, set->round, 1);}
3516
 
3517
/* ------------------------------------------------------------------ */
3518
/* decFloatXor -- logical digitwise XOR of two decFloats              */
3519
/*                                                                    */
3520
/*   result gets the result of XORing dfl and dfr                     */
3521
/*   dfl    is the first decFloat (lhs)                               */
3522
/*   dfr    is the second decFloat (rhs)                              */
3523
/*   set    is the context                                            */
3524
/*   returns result, which will be canonical with sign=0              */
3525
/*                                                                    */
3526
/* The operands must be positive, finite with exponent q=0, and       */
3527
/* comprise just zeros and ones; if not, Invalid operation results.   */
3528
/* ------------------------------------------------------------------ */
3529
decFloat * decFloatXor(decFloat *result,
3530
                       const decFloat *dfl, const decFloat *dfr,
3531
                       decContext *set) {
3532
  if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
3533
   || !DFISCC01(dfl)   || !DFISCC01(dfr)) return decInvalid(result, set);
3534
  /* the operands are positive finite integers (q=0) with just 0s and 1s */
3535
  #if DOUBLE
3536
   DFWORD(result, 0)=ZEROWORD
3537
                   |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
3538
   DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
3539
  #elif QUAD
3540
   DFWORD(result, 0)=ZEROWORD
3541
                   |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
3542
   DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
3543
   DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
3544
   DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
3545
  #endif
3546
  return result;
3547
  } /* decFloatXor */
3548
 
3549
/* ------------------------------------------------------------------ */
3550
/* decInvalid -- set Invalid_operation result                         */
3551
/*                                                                    */
3552
/*   result gets a canonical NaN                                      */
3553
/*   set    is the context                                            */
3554
/*   returns result                                                   */
3555
/*                                                                    */
3556
/* status has Invalid_operation added                                 */
3557
/* ------------------------------------------------------------------ */
3558
static decFloat *decInvalid(decFloat *result, decContext *set) {
3559
  decFloatZero(result);
3560
  DFWORD(result, 0)=DECFLOAT_qNaN;
3561
  set->status|=DEC_Invalid_operation;
3562
  return result;
3563
  } /* decInvalid */
3564
 
3565
/* ------------------------------------------------------------------ */
3566
/* decInfinity -- set canonical Infinity with sign from a decFloat    */
3567
/*                                                                    */
3568
/*   result gets a canonical Infinity                                 */
3569
/*   df     is source decFloat (only the sign is used)                */
3570
/*   returns result                                                   */
3571
/*                                                                    */
3572
/* df may be the same as result                                       */
3573
/* ------------------------------------------------------------------ */
3574
static decFloat *decInfinity(decFloat *result, const decFloat *df) {
3575
  uInt sign=DFWORD(df, 0);          /* save source signword */
3576
  decFloatZero(result);            /* clear everything */
3577
  DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
3578
  return result;
3579
  } /* decInfinity */
3580
 
3581
/* ------------------------------------------------------------------ */
3582
/* decNaNs -- handle NaN argument(s)                                  */
3583
/*                                                                    */
3584
/*   result gets the result of handling dfl and dfr, one or both of   */
3585
/*          which is a NaN                                            */
3586
/*   dfl    is the first decFloat (lhs)                               */
3587
/*   dfr    is the second decFloat (rhs) -- may be NULL for a single- */
3588
/*          operand operation                                         */
3589
/*   set    is the context                                            */
3590
/*   returns result                                                   */
3591
/*                                                                    */
3592
/* Called when one or both operands is a NaN, and propagates the      */
3593
/* appropriate result to res.  When an sNaN is found, it is changed   */
3594
/* to a qNaN and Invalid operation is set.                            */
3595
/* ------------------------------------------------------------------ */
3596
static decFloat *decNaNs(decFloat *result,
3597
                         const decFloat *dfl, const decFloat *dfr,
3598
                         decContext *set) {
3599
  /* handle sNaNs first */
3600
  if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
3601
  if (DFISSNAN(dfl)) {
3602
    decCanonical(result, dfl);          /* propagate canonical sNaN */
3603
    DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
3604
    set->status|=DEC_Invalid_operation;
3605
    return result;
3606
    }
3607
  /* one or both is a quiet NaN */
3608
  if (!DFISNAN(dfl)) dfl=dfr;           /* RHS must be NaN, use it */
3609
  return decCanonical(result, dfl);     /* propagate canonical qNaN */
3610
  } /* decNaNs */
3611
 
3612
/* ------------------------------------------------------------------ */
3613
/* decNumCompare -- numeric comparison of two decFloats               */
3614
/*                                                                    */
3615
/*   dfl    is the left-hand decFloat, which is not a NaN             */
3616
/*   dfr    is the right-hand decFloat, which is not a NaN            */
3617
/*   tot    is 1 for total order compare, 0 for simple numeric        */
3618
/*   returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr               */
3619
/*                                                                    */
3620
/* No error is possible; status and mode are unchanged.               */
3621
/* ------------------------------------------------------------------ */
3622
static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
3623
  Int   sigl, sigr;                     /* LHS and RHS non-0 signums */
3624
  Int   shift;                          /* shift needed to align operands */
3625
  uByte *ub, *uc;                       /* work */
3626
  uInt  uiwork;                         /* for macros */
3627
  /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
3628
  uByte bufl[DECPMAX*2+QUAD*2+4];       /* for LHS coefficient + padding */
3629
  uByte bufr[DECPMAX*2+QUAD*2+4];       /* for RHS coefficient + padding */
3630
 
3631
  sigl=1;
3632
  if (DFISSIGNED(dfl)) {
3633
    if (!DFISSIGNED(dfr)) {             /* -LHS +RHS */
3634
      if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3635
      return -1;                        /* RHS wins */
3636
      }
3637
    sigl=-1;
3638
    }
3639
  if (DFISSIGNED(dfr)) {
3640
    if (!DFISSIGNED(dfl)) {             /* +LHS -RHS */
3641
      if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
3642
      return +1;                        /* LHS wins */
3643
      }
3644
    }
3645
 
3646
  /* signs are the same; operand(s) could be zero */
3647
  sigr=-sigl;                           /* sign to return if abs(RHS) wins */
3648
 
3649
  if (DFISINF(dfl)) {
3650
    if (DFISINF(dfr)) return 0;  /* both infinite & same sign */
3651
    return sigl;                        /* inf > n */
3652
    }
3653
  if (DFISINF(dfr)) return sigr;        /* n < inf [dfl is finite] */
3654
 
3655
  /* here, both are same sign and finite; calculate their offset */
3656
  shift=GETEXP(dfl)-GETEXP(dfr);        /* [0 means aligned] */
3657
  /* [bias can be ignored -- the absolute exponent is not relevant] */
3658
 
3659
  if (DFISZERO(dfl)) {
3660
    if (!DFISZERO(dfr)) return sigr;    /* LHS=0, RHS!=0 */
3661
    /* both are zero, return 0 if both same exponent or numeric compare */
3662
    if (shift==0 || !tot) return 0;
3663
    if (shift>0) return sigl;
3664
    return sigr;                        /* [shift<0] */
3665
    }
3666
   else {                               /* LHS!=0 */
3667
    if (DFISZERO(dfr)) return sigl;     /* LHS!=0, RHS=0 */
3668
    }
3669
  /* both are known to be non-zero at this point */
3670
 
3671
  /* if the exponents are so different that the coefficients do not */
3672
  /* overlap (by even one digit) then a full comparison is not needed */
3673
  if (abs(shift)>=DECPMAX) {            /* no overlap */
3674
    /* coefficients are known to be non-zero */
3675
    if (shift>0) return sigl;
3676
    return sigr;                        /* [shift<0] */
3677
    }
3678
 
3679
  /* decode the coefficients */
3680
  /* (shift both right two if Quad to make a multiple of four) */
3681
  #if QUAD
3682
    UBFROMUI(bufl, 0);
3683
    UBFROMUI(bufr, 0);
3684
  #endif
3685
  GETCOEFF(dfl, bufl+QUAD*2);           /* decode from decFloat */
3686
  GETCOEFF(dfr, bufr+QUAD*2);           /* .. */
3687
  if (shift==0) {                        /* aligned; common and easy */
3688
    /* all multiples of four, here */
3689
    for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
3690
      uInt ui=UBTOUI(ub);
3691
      if (ui==UBTOUI(uc)) continue;     /* so far so same */
3692
      /* about to find a winner; go by bytes in case little-endian */
3693
      for (;; ub++, uc++) {
3694
        if (*ub>*uc) return sigl;       /* difference found */
3695
        if (*ub<*uc) return sigr;       /* .. */
3696
        }
3697
      }
3698
    } /* aligned */
3699
   else if (shift>0) {                   /* lhs to left */
3700
    ub=bufl;                            /* RHS pointer */
3701
    /* pad bufl so right-aligned; most shifts will fit in 8 */
3702
    UBFROMUI(bufl+DECPMAX+QUAD*2, 0);    /* add eight zeros */
3703
    UBFROMUI(bufl+DECPMAX+QUAD*2+4, 0); /* .. */
3704
    if (shift>8) {
3705
      /* more than eight; fill the rest, and also worth doing the */
3706
      /* lead-in by fours */
3707
      uByte *up;                        /* work */
3708
      uByte *upend=bufl+DECPMAX+QUAD*2+shift;
3709
      for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
3710
      /* [pads up to 36 in all for Quad] */
3711
      for (;; ub+=4) {
3712
        if (UBTOUI(ub)!=0) return sigl;
3713
        if (ub+4>bufl+shift-4) break;
3714
        }
3715
      }
3716
    /* check remaining leading digits */
3717
    for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
3718
    /* now start the overlapped part; bufl has been padded, so the */
3719
    /* comparison can go for the full length of bufr, which is a */
3720
    /* multiple of 4 bytes */
3721
    for (uc=bufr; ; uc+=4, ub+=4) {
3722
      uInt ui=UBTOUI(ub);
3723
      if (ui!=UBTOUI(uc)) {             /* mismatch found */
3724
        for (;; uc++, ub++) {           /* check from left [little-endian?] */
3725
          if (*ub>*uc) return sigl;     /* difference found */
3726
          if (*ub<*uc) return sigr;     /* .. */
3727
          }
3728
        } /* mismatch */
3729
      if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
3730
      }
3731
    } /* shift>0 */
3732
 
3733
   else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
3734
    uc=bufr;                            /* RHS pointer */
3735
    /* pad bufr so right-aligned; most shifts will fit in 8 */
3736
    UBFROMUI(bufr+DECPMAX+QUAD*2, 0);    /* add eight zeros */
3737
    UBFROMUI(bufr+DECPMAX+QUAD*2+4, 0); /* .. */
3738
    if (shift<-8) {
3739
      /* more than eight; fill the rest, and also worth doing the */
3740
      /* lead-in by fours */
3741
      uByte *up;                        /* work */
3742
      uByte *upend=bufr+DECPMAX+QUAD*2-shift;
3743
      for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
3744
      /* [pads up to 36 in all for Quad] */
3745
      for (;; uc+=4) {
3746
        if (UBTOUI(uc)!=0) return sigr;
3747
        if (uc+4>bufr-shift-4) break;
3748
        }
3749
      }
3750
    /* check remaining leading digits */
3751
    for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
3752
    /* now start the overlapped part; bufr has been padded, so the */
3753
    /* comparison can go for the full length of bufl, which is a */
3754
    /* multiple of 4 bytes */
3755
    for (ub=bufl; ; ub+=4, uc+=4) {
3756
      uInt ui=UBTOUI(ub);
3757
      if (ui!=UBTOUI(uc)) {             /* mismatch found */
3758
        for (;; ub++, uc++) {           /* check from left [little-endian?] */
3759
          if (*ub>*uc) return sigl;     /* difference found */
3760
          if (*ub<*uc) return sigr;     /* .. */
3761
          }
3762
        } /* mismatch */
3763
      if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
3764
      }
3765
    } /* shift<0 */
3766
 
3767
  /* Here when compare equal */
3768
  if (!tot) return 0;                    /* numerically equal */
3769
  /* total ordering .. exponent matters */
3770
  if (shift>0) return sigl;              /* total order by exponent */
3771
  if (shift<0) return sigr;              /* .. */
3772
  return 0;
3773
  } /* decNumCompare */
3774
 
3775
/* ------------------------------------------------------------------ */
3776
/* decToInt32 -- local routine to effect ToInteger conversions        */
3777
/*                                                                    */
3778
/*   df     is the decFloat to convert                                */
3779
/*   set    is the context                                            */
3780
/*   rmode  is the rounding mode to use                               */
3781
/*   exact  is 1 if Inexact should be signalled                       */
3782
/*   unsign is 1 if the result a uInt, 0 if an Int (cast to uInt)     */
3783
/*   returns 32-bit result as a uInt                                  */
3784
/*                                                                    */
3785
/* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
3786
/* these cases 0 is returned.                                         */
3787
/* ------------------------------------------------------------------ */
3788
static uInt decToInt32(const decFloat *df, decContext *set,
3789
                       enum rounding rmode, Flag exact, Flag unsign) {
3790
  Int  exp;                        /* exponent */
3791
  uInt sourhi, sourpen, sourlo;    /* top word from source decFloat .. */
3792
  uInt hi, lo;                     /* .. penultimate, least, etc. */
3793
  decFloat zero, result;           /* work */
3794
  Int  i;                          /* .. */
3795
 
3796
  /* Start decoding the argument */
3797
  sourhi=DFWORD(df, 0);          /* top word */
3798
  exp=DECCOMBEXP[sourhi>>26];           /* get exponent high bits (in place) */
3799
  if (EXPISSPECIAL(exp)) {              /* is special? */
3800
    set->status|=DEC_Invalid_operation; /* signal */
3801
    return 0;
3802
    }
3803
 
3804
  /* Here when the argument is finite */
3805
  if (GETEXPUN(df)==0) result=*df;       /* already a true integer */
3806
   else {                               /* need to round to integer */
3807
    enum rounding saveround;            /* saver */
3808
    uInt savestatus;                    /* .. */
3809
    saveround=set->round;               /* save rounding mode .. */
3810
    savestatus=set->status;             /* .. and status */
3811
    set->round=rmode;                   /* set mode */
3812
    decFloatZero(&zero);                /* make 0E+0 */
3813
    set->status=0;                       /* clear */
3814
    decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
3815
    set->round=saveround;               /* restore rounding mode .. */
3816
    if (exact) set->status|=savestatus; /* include Inexact */
3817
     else set->status=savestatus;       /* .. or just original status */
3818
    }
3819
 
3820
  /* only the last four declets of the coefficient can contain */
3821
  /* non-zero; check for others (and also NaN or Infinity from the */
3822
  /* Quantize) first (see DFISZERO for explanation): */
3823
  /* decFloatShow(&result, "sofar"); */
3824
  #if DOUBLE
3825
  if ((DFWORD(&result, 0)&0x1c03ff00)!=0
3826
   || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3827
  #elif QUAD
3828
  if ((DFWORD(&result, 2)&0xffffff00)!=0
3829
   ||  DFWORD(&result, 1)!=0
3830
   || (DFWORD(&result, 0)&0x1c003fff)!=0
3831
   || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
3832
  #endif
3833
    set->status|=DEC_Invalid_operation; /* Invalid or out of range */
3834
    return 0;
3835
    }
3836
  /* get last twelve digits of the coefficent into hi & ho, base */
3837
  /* 10**9 (see GETCOEFFBILL): */
3838
  sourlo=DFWORD(&result, DECWORDS-1);
3839
  lo=DPD2BIN0[sourlo&0x3ff]
3840
    +DPD2BINK[(sourlo>>10)&0x3ff]
3841
    +DPD2BINM[(sourlo>>20)&0x3ff];
3842
  sourpen=DFWORD(&result, DECWORDS-2);
3843
  hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
3844
 
3845
  /* according to request, check range carefully */
3846
  if (unsign) {
3847
    if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
3848
      set->status|=DEC_Invalid_operation; /* out of range */
3849
      return 0;
3850
      }
3851
    return hi*BILLION+lo;
3852
    }
3853
  /* signed */
3854
  if (hi>2 || (hi==2 && lo>147483647)) {
3855
    /* handle the usual edge case */
3856
    if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
3857
    set->status|=DEC_Invalid_operation; /* truly out of range */
3858
    return 0;
3859
    }
3860
  i=hi*BILLION+lo;
3861
  if (DFISSIGNED(&result)) i=-i;
3862
  return (uInt)i;
3863
  } /* decToInt32 */
3864
 
3865
/* ------------------------------------------------------------------ */
3866
/* decToIntegral -- local routine to effect ToIntegral value          */
3867
/*                                                                    */
3868
/*   result gets the result                                           */
3869
/*   df     is the decFloat to round                                  */
3870
/*   set    is the context                                            */
3871
/*   rmode  is the rounding mode to use                               */
3872
/*   exact  is 1 if Inexact should be signalled                       */
3873
/*   returns result                                                   */
3874
/* ------------------------------------------------------------------ */
3875
static decFloat * decToIntegral(decFloat *result, const decFloat *df,
3876
                                decContext *set, enum rounding rmode,
3877
                                Flag exact) {
3878
  Int  exp;                        /* exponent */
3879
  uInt sourhi;                     /* top word from source decFloat */
3880
  enum rounding saveround;         /* saver */
3881
  uInt savestatus;                 /* .. */
3882
  decFloat zero;                   /* work */
3883
 
3884
  /* Start decoding the argument */
3885
  sourhi=DFWORD(df, 0);     /* top word */
3886
  exp=DECCOMBEXP[sourhi>>26];      /* get exponent high bits (in place) */
3887
 
3888
  if (EXPISSPECIAL(exp)) {         /* is special? */
3889
    /* NaNs are handled as usual */
3890
    if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
3891
    /* must be infinite; return canonical infinity with sign of df */
3892
    return decInfinity(result, df);
3893
    }
3894
 
3895
  /* Here when the argument is finite */
3896
  /* complete extraction of the exponent */
3897
  exp+=GETECON(df)-DECBIAS;             /* .. + continuation and unbias */
3898
 
3899
  if (exp>=0) return decCanonical(result, df); /* already integral */
3900
 
3901
  saveround=set->round;                 /* save rounding mode .. */
3902
  savestatus=set->status;               /* .. and status */
3903
  set->round=rmode;                     /* set mode */
3904
  decFloatZero(&zero);                  /* make 0E+0 */
3905
  decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
3906
  set->round=saveround;                 /* restore rounding mode .. */
3907
  if (!exact) set->status=savestatus;   /* .. and status, unless exact */
3908
  return result;
3909
  } /* decToIntegral */

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.