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

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [gcc-4.2.2/] [libdecnumber/] [decNumber.c] - Blame information for rev 455

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

Line No. Rev Author Line
1 38 julius
/* Decimal Number module for the decNumber C Library
2
   Copyright (C) 2005 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 2, or (at your option) any later
10
   version.
11
 
12
   In addition to the permissions in the GNU General Public License,
13
   the Free Software Foundation gives you unlimited permission to link
14
   the compiled version of this file into combinations with other
15
   programs, and to distribute those combinations without any
16
   restriction coming from the use of this file.  (The General Public
17
   License restrictions do apply in other respects; for example, they
18
   cover modification of the file, and distribution when not linked
19
   into a combine executable.)
20
 
21
   GCC is distributed in the hope that it will be useful, but WITHOUT ANY
22
   WARRANTY; without even the implied warranty of MERCHANTABILITY or
23
   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
24
   for more details.
25
 
26
   You should have received a copy of the GNU General Public License
27
   along with GCC; see the file COPYING.  If not, write to the Free
28
   Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
29
   02110-1301, USA.  */
30
 
31
/* ------------------------------------------------------------------ */
32
/* This module comprises the routines for Standard Decimal Arithmetic */
33
/* as defined in the specification which may be found on the          */
34
/* http://www2.hursley.ibm.com/decimal web pages.  It implements both */
35
/* the full ('extended') arithmetic and the simpler ('subset')        */
36
/* arithmetic.                                                        */
37
/*                                                                    */
38
/* Usage notes:                                                       */
39
/*                                                                    */
40
/* 1. This code is ANSI C89 except:                                   */
41
/*                                                                    */
42
/*    a) Line comments (double forward slash) are used.  (Most C      */
43
/*       compilers accept these.  If yours does not, a simple script  */
44
/*       can be used to convert them to ANSI C comments.)             */
45
/*                                                                    */
46
/*    b) Types from C99 stdint.h are used.  If you do not have this   */
47
/*       header file, see the User's Guide section of the decNumber   */
48
/*       documentation; this lists the necessary definitions.         */
49
/*                                                                    */
50
/*    c) If DECDPUN>4, non-ANSI 64-bit 'long long' types are used.    */
51
/*       To avoid these, set DECDPUN <= 4 (see documentation).        */
52
/*                                                                    */
53
/* 2. The decNumber format which this library uses is optimized for   */
54
/*    efficient processing of relatively short numbers; in particular */
55
/*    it allows the use of fixed sized structures and minimizes copy  */
56
/*    and move operations.  It does, however, support arbitrary       */
57
/*    precision (up to 999,999,999 digits) and arbitrary exponent     */
58
/*    range (Emax in the range 0 through 999,999,999 and Emin in the  */
59
/*    range -999,999,999 through 0).                                  */
60
/*                                                                    */
61
/* 3. Operands to operator functions are never modified unless they   */
62
/*    are also specified to be the result number (which is always     */
63
/*    permitted).  Other than that case, operands may not overlap.    */
64
/*                                                                    */
65
/* 4. Error handling: the type of the error is ORed into the status   */
66
/*    flags in the current context (decContext structure).  The       */
67
/*    SIGFPE signal is then raised if the corresponding trap-enabler  */
68
/*    flag in the decContext is set (is 1).                           */
69
/*                                                                    */
70
/*    It is the responsibility of the caller to clear the status      */
71
/*    flags as required.                                              */
72
/*                                                                    */
73
/*    The result of any routine which returns a number will always    */
74
/*    be a valid number (which may be a special value, such as an     */
75
/*    Infinity or NaN).                                               */
76
/*                                                                    */
77
/* 5. The decNumber format is not an exchangeable concrete            */
78
/*    representation as it comprises fields which may be machine-     */
79
/*    dependent (big-endian or little-endian, for example).           */
80
/*    Canonical conversions to and from strings are provided; other   */
81
/*    conversions are available in separate modules.                  */
82
/*                                                                    */
83
/* 6. Normally, input operands are assumed to be valid.  Set DECCHECK */
84
/*    to 1 for extended operand checking (including NULL operands).   */
85
/*    Results are undefined if a badly-formed structure (or a NULL    */
86
/*    NULL pointer to a structure) is provided, though with DECCHECK  */
87
/*    enabled the operator routines are protected against exceptions. */
88
/*    (Except if the result pointer is NULL, which is unrecoverable.) */
89
/*                                                                    */
90
/*    However, the routines will never cause exceptions if they are   */
91
/*    given well-formed operands, even if the value of the operands   */
92
/*    is inappropriate for the operation and DECCHECK is not set.     */
93
/*                                                                    */
94
/* 7. Subset arithmetic is available only if DECSUBSET is set to 1.   */
95
/* ------------------------------------------------------------------ */
96
/* Implementation notes for maintenance of this module:               */
97
/*                                                                    */
98
/* 1. Storage leak protection:  Routines which use malloc are not     */
99
/*    permitted to use return for fastpath or error exits (i.e.,      */
100
/*    they follow strict structured programming conventions).         */
101
/*    Instead they have a do{}while(0); construct surrounding the     */
102
/*    code which is protected -- break may be used from this.         */
103
/*    Other routines are allowed to use the return statement inline.  */
104
/*                                                                    */
105
/*    Storage leak accounting can be enabled using DECALLOC.          */
106
/*                                                                    */
107
/* 2. All loops use the for(;;) construct.  Any do construct is for   */
108
/*    protection as just described.                                   */
109
/*                                                                    */
110
/* 3. Setting status in the context must always be the very last      */
111
/*    action in a routine, as non-0 status may raise a trap and hence */
112
/*    the call to set status may not return (if the handler uses long */
113
/*    jump).  Therefore all cleanup must be done first.  In general,  */
114
/*    to achieve this we accumulate status and only finally apply it  */
115
/*    by calling decContextSetStatus (via decStatus).                 */
116
/*                                                                    */
117
/*    Routines which allocate storage cannot, therefore, use the      */
118
/*    'top level' routines which could cause a non-returning          */
119
/*    transfer of control.  The decXxxxOp routines are safe (do not   */
120
/*    call decStatus even if traps are set in the context) and should */
121
/*    be used instead (they are also a little faster).                */
122
/*                                                                    */
123
/* 4. Exponent checking is minimized by allowing the exponent to      */
124
/*    grow outside its limits during calculations, provided that      */
125
/*    the decFinalize function is called later.  Multiplication and   */
126
/*    division, and intermediate calculations in exponentiation,      */
127
/*    require more careful checks because of the risk of 31-bit       */
128
/*    overflow (the most negative valid exponent is -1999999997, for  */
129
/*    a 999999999-digit number with adjusted exponent of -999999999). */
130
/*                                                                    */
131
/* 5. Rounding is deferred until finalization of results, with any    */
132
/*    'off to the right' data being represented as a single digit     */
133
/*    residue (in the range -1 through 9).  This avoids any double-   */
134
/*    rounding when more than one shortening takes place (for         */
135
/*    example, when a result is subnormal).                           */
136
/*                                                                    */
137
/* 6. The digits count is allowed to rise to a multiple of DECDPUN    */
138
/*    during many operations, so whole Units are handled and exact    */
139
/*    accounting of digits is not needed.  The correct digits value   */
140
/*    is found by decGetDigits, which accounts for leading zeros.     */
141
/*    This must be called before any rounding if the number of digits */
142
/*    is not known exactly.                                           */
143
/*                                                                    */
144
/* 7. We use the multiply-by-reciprocal 'trick' for partitioning      */
145
/*    numbers up to four digits, using appropriate constants.  This   */
146
/*    is not useful for longer numbers because overflow of 32 bits    */
147
/*    would lead to 4 multiplies, which is almost as expensive as     */
148
/*    a divide (unless we assumed floating-point multiply available). */
149
/*                                                                    */
150
/* 8. Unusual abbreviations possibly used in the commentary:          */
151
/*      lhs -- left hand side (operand, of an operation)              */
152
/*      lsd -- least significant digit (of coefficient)               */
153
/*      lsu -- least significant Unit (of coefficient)                */
154
/*      msd -- most significant digit (of coefficient)                */
155
/*      msu -- most significant Unit (of coefficient)                 */
156
/*      rhs -- right hand side (operand, of an operation)             */
157
/*      +ve -- positive                                               */
158
/*      -ve -- negative                                               */
159
/* ------------------------------------------------------------------ */
160
 
161
/* Some of glibc's string inlines cause warnings.  Plus we'd rather
162
   rely on (and therefore test) GCC's string builtins.  */
163
#define __NO_STRING_INLINES
164
 
165
#include <stdlib.h>             /* for malloc, free, etc. */
166
#include <stdio.h>              /* for printf [if needed] */
167
#include <string.h>             /* for strcpy */
168
#include <ctype.h>              /* for lower */
169
#include "config.h"
170
#include "decNumber.h"          /* base number library */
171
#include "decNumberLocal.h"     /* decNumber local types, etc. */
172
 
173
/* Constants */
174
/* Public constant array: powers of ten (powers[n]==10**n) */
175
const uInt powers[] = { 1, 10, 100, 1000, 10000, 100000, 1000000,
176
  10000000, 100000000, 1000000000
177
};
178
 
179
/* Local constants */
180
#define DIVIDE    0x80          /* Divide operators */
181
#define REMAINDER 0x40          /* .. */
182
#define DIVIDEINT 0x20          /* .. */
183
#define REMNEAR   0x10          /* .. */
184
#define COMPARE   0x01          /* Compare operators */
185
#define COMPMAX   0x02          /* .. */
186
#define COMPMIN   0x03          /* .. */
187
#define COMPNAN   0x04          /* .. [NaN processing] */
188
 
189
#define DEC_sNaN 0x40000000     /* local status: sNaN signal */
190
#define BADINT (Int)0x80000000  /* most-negative Int; error indicator */
191
 
192
static Unit one[] = { 1 };      /* Unit array of 1, used for incrementing */
193
 
194
/* Granularity-dependent code */
195
#if DECDPUN<=4
196
#define eInt  Int               /* extended integer */
197
#define ueInt uInt              /* unsigned extended integer */
198
  /* Constant multipliers for divide-by-power-of five using reciprocal */
199
  /* multiply, after removing powers of 2 by shifting, and final shift */
200
  /* of 17 [we only need up to **4] */
201
static const uInt multies[] = { 131073, 26215, 5243, 1049, 210 };
202
 
203
  /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
204
#define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
205
#else
206
  /* For DECDPUN>4 we currently use non-ANSI 64-bit types.  These could */
207
  /* be replaced by subroutine calls later. */
208
#ifdef long
209
#undef long
210
#endif
211
typedef signed long long Long;
212
typedef unsigned long long uLong;
213
#define eInt  Long              /* extended integer */
214
#define ueInt uLong             /* unsigned extended integer */
215
#endif
216
 
217
/* Local routines */
218
static decNumber *decAddOp (decNumber *, const decNumber *,
219
                            const decNumber *, decContext *,
220
                            uByte, uInt *);
221
static void decApplyRound (decNumber *, decContext *, Int, uInt *);
222
static Int decCompare (const decNumber * lhs, const decNumber * rhs);
223
static decNumber *decCompareOp (decNumber *, const decNumber *, const decNumber *,
224
                                decContext *, Flag, uInt *);
225
static void decCopyFit (decNumber *, const decNumber *, decContext *,
226
                        Int *, uInt *);
227
static decNumber *decDivideOp (decNumber *, const decNumber *, const decNumber *,
228
                               decContext *, Flag, uInt *);
229
static void decFinalize (decNumber *, decContext *, Int *, uInt *);
230
static Int decGetDigits (const Unit *, Int);
231
#if DECSUBSET
232
static Int decGetInt (const decNumber *, decContext *);
233
#else
234
static Int decGetInt (const decNumber *);
235
#endif
236
static decNumber *decMultiplyOp (decNumber *, const decNumber *,
237
                                 const decNumber *, decContext *, uInt *);
238
static decNumber *decNaNs (decNumber *, const decNumber *, const decNumber *, uInt *);
239
static decNumber *decQuantizeOp (decNumber *, const decNumber *,
240
                                 const decNumber *, decContext *, Flag, uInt *);
241
static void decSetCoeff (decNumber *, decContext *, const Unit *,
242
                         Int, Int *, uInt *);
243
static void decSetOverflow (decNumber *, decContext *, uInt *);
244
static void decSetSubnormal (decNumber *, decContext *, Int *, uInt *);
245
static Int decShiftToLeast (Unit *, Int, Int);
246
static Int decShiftToMost (Unit *, Int, Int);
247
static void decStatus (decNumber *, uInt, decContext *);
248
static Flag decStrEq (const char *, const char *);
249
static void decToString (const decNumber *, char[], Flag);
250
static decNumber *decTrim (decNumber *, Flag, Int *);
251
static Int decUnitAddSub (const Unit *, Int, const Unit *, Int, Int, Unit *, Int);
252
static Int decUnitCompare (const Unit *, Int, const Unit *, Int, Int);
253
 
254
#if !DECSUBSET
255
/* decFinish == decFinalize when no subset arithmetic needed */
256
#define decFinish(a,b,c,d) decFinalize(a,b,c,d)
257
#else
258
static void decFinish (decNumber *, decContext *, Int *, uInt *);
259
static decNumber *decRoundOperand (const decNumber *, decContext *, uInt *);
260
#endif
261
 
262
/* Diagnostic macros, etc. */
263
#if DECALLOC
264
/* Handle malloc/free accounting.  If enabled, our accountable routines */
265
/* are used; otherwise the code just goes straight to the system malloc */
266
/* and free routines. */
267
#define malloc(a) decMalloc(a)
268
#define free(a) decFree(a)
269
#define DECFENCE 0x5a           /* corruption detector */
270
/* 'Our' malloc and free: */
271
static void *decMalloc (size_t);
272
static void decFree (void *);
273
uInt decAllocBytes = 0;          /* count of bytes allocated */
274
/* Note that DECALLOC code only checks for storage buffer overflow. */
275
/* To check for memory leaks, the decAllocBytes variable should be */
276
/* checked to be 0 at appropriate times (e.g., after the test */
277
/* harness completes a set of tests).  This checking may be unreliable */
278
/* if the testing is done in a multi-thread environment. */
279
#endif
280
 
281
#if DECCHECK
282
/* Optional operand checking routines.  Enabling these means that */
283
/* decNumber and decContext operands to operator routines are checked */
284
/* for correctness.  This roughly doubles the execution time of the */
285
/* fastest routines (and adds 600+ bytes), so should not normally be */
286
/* used in 'production'. */
287
#define DECUNUSED (void *)(0xffffffff)
288
static Flag decCheckOperands (decNumber *, const decNumber *,
289
                              const decNumber *, decContext *);
290
static Flag decCheckNumber (const decNumber *, decContext *);
291
#endif
292
 
293
#if DECTRACE || DECCHECK
294
/* Optional trace/debugging routines. */
295
void decNumberShow (const decNumber *); /* displays the components of a number */
296
static void decDumpAr (char, const Unit *, Int);
297
#endif
298
 
299
/* ================================================================== */
300
/* Conversions                                                        */
301
/* ================================================================== */
302
 
303 152 jeremybenn
/* The following two functions are copied from later versions, since  */
304
/* they are needed by GDB which shares this library.                  */
305
 
306 38 julius
/* ------------------------------------------------------------------ */
307 152 jeremybenn
/* from-int32 -- conversion from Int or uInt                          */
308
/*                                                                    */
309
/*  dn is the decNumber to receive the integer                        */
310
/*  in or uin is the integer to be converted                          */
311
/*  returns dn                                                        */
312
/*                                                                    */
313
/* No error is possible.                                              */
314
/* ------------------------------------------------------------------ */
315
decNumber * decNumberFromInt32(decNumber *dn, Int in) {
316
  uInt unsig;
317
  if (in>=0) unsig=in;
318
   else {                               /* negative (possibly BADINT) */
319
    if (in==BADINT) unsig=(uInt)1073741824*2; /* special case */
320
     else unsig=-in;                    /* invert */
321
    }
322
  /* in is now positive */
323
  decNumberFromUInt32(dn, unsig);
324
  if (in<0) dn->bits=DECNEG;             /* sign needed */
325
  return dn;
326
  } /* decNumberFromInt32 */
327
 
328
decNumber * decNumberFromUInt32(decNumber *dn, uInt uin) {
329
  Unit *up;                             /* work pointer */
330
  decNumberZero(dn);                    /* clean */
331
  if (uin==0) return dn;         /* [or decGetDigits bad call] */
332
  for (up=dn->lsu; uin>0; up++) {
333
    *up=(Unit)(uin%(DECDPUNMAX+1));
334
    uin=uin/(DECDPUNMAX+1);
335
    }
336
  dn->digits=decGetDigits(dn->lsu, up-dn->lsu);
337
  return dn;
338
  } /* decNumberFromUInt32 */
339
 
340
 
341
/* ------------------------------------------------------------------ */
342 38 julius
/* to-scientific-string -- conversion to numeric string               */
343
/* to-engineering-string -- conversion to numeric string              */
344
/*                                                                    */
345
/*   decNumberToString(dn, string);                                   */
346
/*   decNumberToEngString(dn, string);                                */
347
/*                                                                    */
348
/*  dn is the decNumber to convert                                    */
349
/*  string is the string where the result will be laid out            */
350
/*                                                                    */
351
/*  string must be at least dn->digits+14 characters long             */
352
/*                                                                    */
353
/*  No error is possible, and no status can be set.                   */
354
/* ------------------------------------------------------------------ */
355
char *
356
decNumberToString (const decNumber * dn, char *string)
357
{
358
  decToString (dn, string, 0);
359
  return string;
360
}
361
 
362
char *
363
decNumberToEngString (const decNumber * dn, char *string)
364
{
365
  decToString (dn, string, 1);
366
  return string;
367
}
368
 
369
/* ------------------------------------------------------------------ */
370
/* to-number -- conversion from numeric string                        */
371
/*                                                                    */
372
/* decNumberFromString -- convert string to decNumber                 */
373
/*   dn        -- the number structure to fill                        */
374
/*   chars[]   -- the string to convert ('\0' terminated)             */
375
/*   set       -- the context used for processing any error,          */
376
/*                determining the maximum precision available         */
377
/*                (set.digits), determining the maximum and minimum   */
378
/*                exponent (set.emax and set.emin), determining if    */
379
/*                extended values are allowed, and checking the       */
380
/*                rounding mode if overflow occurs or rounding is     */
381
/*                needed.                                             */
382
/*                                                                    */
383
/* The length of the coefficient and the size of the exponent are     */
384
/* checked by this routine, so the correct error (Underflow or        */
385
/* Overflow) can be reported or rounding applied, as necessary.       */
386
/*                                                                    */
387
/* If bad syntax is detected, the result will be a quiet NaN.         */
388
/* ------------------------------------------------------------------ */
389
decNumber *
390
decNumberFromString (decNumber * dn, const char chars[], decContext * set)
391
{
392
  Int exponent = 0;              /* working exponent [assume 0] */
393
  uByte bits = 0;                /* working flags [assume +ve] */
394
  Unit *res;                    /* where result will be built */
395
  Unit resbuff[D2U (DECBUFFER + 1)];    /* local buffer in case need temporary */
396
  Unit *allocres = NULL;        /* -> allocated result, iff allocated */
397
  Int need;                     /* units needed for result */
398
  Int d = 0;                     /* count of digits found in decimal part */
399
  const char *dotchar = NULL;   /* where dot was found */
400
  const char *cfirst;           /* -> first character of decimal part */
401
  const char *last = NULL;      /* -> last digit of decimal part */
402
  const char *firstexp;         /* -> first significant exponent digit */
403
  const char *c;                /* work */
404
  Unit *up;                     /* .. */
405
#if DECDPUN>1
406
  Int i;                        /* .. */
407
#endif
408
  Int residue = 0;               /* rounding residue */
409
  uInt status = 0;               /* error code */
410
 
411
#if DECCHECK
412
  if (decCheckOperands (DECUNUSED, DECUNUSED, DECUNUSED, set))
413
    return decNumberZero (dn);
414
#endif
415
 
416
  do
417
    {                           /* status & malloc protection */
418
      c = chars;                /* -> input character */
419
      if (*c == '-')
420
        {                       /* handle leading '-' */
421
          bits = DECNEG;
422
          c++;
423
        }
424
      else if (*c == '+')
425
        c++;                    /* step over leading '+' */
426
      /* We're at the start of the number [we think] */
427
      cfirst = c;               /* save */
428
      for (;; c++)
429
        {
430
          if (*c >= '0' && *c <= '9')
431
            {                   /* test for Arabic digit */
432
              last = c;
433
              d++;              /* count of real digits */
434
              continue;         /* still in decimal part */
435
            }
436
          if (*c != '.')
437
            break;              /* done with decimal part */
438
          /* dot: record, check, and ignore */
439
          if (dotchar != NULL)
440
            {                   /* two dots */
441
              last = NULL;      /* indicate bad */
442
              break;
443
            }                   /* .. and go report */
444
          dotchar = c;          /* offset into decimal part */
445
        }                       /* c */
446
 
447
      if (last == NULL)
448
        {                       /* no decimal digits, or >1 . */
449
#if DECSUBSET
450
          /* If subset then infinities and NaNs are not allowed */
451
          if (!set->extended)
452
            {
453
              status = DEC_Conversion_syntax;
454
              break;            /* all done */
455
            }
456
          else
457
            {
458
#endif
459
              /* Infinities and NaNs are possible, here */
460
              decNumberZero (dn);       /* be optimistic */
461
              if (decStrEq (c, "Infinity") || decStrEq (c, "Inf"))
462
                {
463
                  dn->bits = bits | DECINF;
464
                  break;        /* all done */
465
                }
466
              else
467
                {               /* a NaN expected */
468
                  /* 2003.09.10 NaNs are now permitted to have a sign */
469
                  status = DEC_Conversion_syntax;       /* assume the worst */
470
                  dn->bits = bits | DECNAN;     /* assume simple NaN */
471
                  if (*c == 's' || *c == 'S')
472
                    {           /* looks like an` sNaN */
473
                      c++;
474
                      dn->bits = bits | DECSNAN;
475
                    }
476
                  if (*c != 'n' && *c != 'N')
477
                    break;      /* check caseless "NaN" */
478
                  c++;
479
                  if (*c != 'a' && *c != 'A')
480
                    break;      /* .. */
481
                  c++;
482
                  if (*c != 'n' && *c != 'N')
483
                    break;      /* .. */
484
                  c++;
485
                  /* now nothing, or nnnn, expected */
486
                  /* -> start of integer and skip leading 0s [including plain 0] */
487
                  for (cfirst = c; *cfirst == '0';)
488
                    cfirst++;
489
                  if (*cfirst == '\0')
490
                    {           /* "NaN" or "sNaN", maybe with all 0s */
491
                      status = 0;        /* it's good */
492
                      break;    /* .. */
493
                    }
494
                  /* something other than 0s; setup last and d as usual [no dots] */
495
                  for (c = cfirst;; c++, d++)
496
                    {
497
                      if (*c < '0' || *c > '9')
498
                        break;  /* test for Arabic digit */
499
                      last = c;
500
                    }
501
                  if (*c != '\0')
502
                    break;      /* not all digits */
503
                  if (d > set->digits)
504
                    break;      /* too many digits */
505
                  /* good; drop through and convert the integer */
506
                  status = 0;
507
                  bits = dn->bits;      /* for copy-back */
508
                }               /* NaN expected */
509
#if DECSUBSET
510
            }
511
#endif
512
        }                       /* last==NULL */
513
 
514
      if (*c != '\0')
515
        {                       /* more there; exponent expected... */
516
          Flag nege = 0; /* 1=negative exponent */
517
          if (*c != 'e' && *c != 'E')
518
            {
519
              status = DEC_Conversion_syntax;
520
              break;
521
            }
522
 
523
          /* Found 'e' or 'E' -- now process explicit exponent */
524
          /* 1998.07.11: sign no longer required */
525
          c++;                  /* to (expected) sign */
526
          if (*c == '-')
527
            {
528
              nege = 1;
529
              c++;
530
            }
531
          else if (*c == '+')
532
            c++;
533
          if (*c == '\0')
534
            {
535
              status = DEC_Conversion_syntax;
536
              break;
537
            }
538
 
539
          for (; *c == '0' && *(c + 1) != '\0';)
540
            c++;                /* strip insignificant zeros */
541
          firstexp = c;         /* save exponent digit place */
542
          for (;; c++)
543
            {
544
              if (*c < '0' || *c > '9')
545
                break;          /* not a digit */
546
              exponent = X10 (exponent) + (Int) * c - (Int) '0';
547
            }                   /* c */
548
          /* if we didn't end on '\0' must not be a digit */
549
          if (*c != '\0')
550
            {
551
              status = DEC_Conversion_syntax;
552
              break;
553
            }
554
 
555
          /* (this next test must be after the syntax check) */
556
          /* if it was too long the exponent may have wrapped, so check */
557
          /* carefully and set it to a certain overflow if wrap possible */
558
          if (c >= firstexp + 9 + 1)
559
            {
560
              if (c > firstexp + 9 + 1 || *firstexp > '1')
561
                exponent = DECNUMMAXE * 2;
562
              /* [up to 1999999999 is OK, for example 1E-1000000998] */
563
            }
564
          if (nege)
565
            exponent = -exponent;       /* was negative */
566
        }                       /* had exponent */
567
      /* Here when all inspected; syntax is good */
568
 
569
      /* Handle decimal point... */
570
      if (dotchar != NULL && dotchar < last)    /* embedded . found, so */
571
        exponent = exponent - (last - dotchar); /* .. adjust exponent */
572
      /* [we can now ignore the .] */
573
 
574
      /* strip leading zeros/dot (leave final if all 0's) */
575
      for (c = cfirst; c < last; c++)
576
        {
577
          if (*c == '0')
578
            d--;                /* 0 stripped */
579
          else if (*c != '.')
580
            break;
581
          cfirst++;             /* step past leader */
582
        }                       /* c */
583
 
584
#if DECSUBSET
585
      /* We can now make a rapid exit for zeros if !extended */
586
      if (*cfirst == '0' && !set->extended)
587
        {
588
          decNumberZero (dn);   /* clean result */
589
          break;                /* [could be return] */
590
        }
591
#endif
592
 
593
      /* OK, the digits string is good.  Copy to the decNumber, or to
594
         a temporary decNumber if rounding is needed */
595
      if (d <= set->digits)
596
        res = dn->lsu;          /* fits into given decNumber */
597
      else
598
        {                       /* rounding needed */
599
          need = D2U (d);       /* units needed */
600
          res = resbuff;        /* assume use local buffer */
601
          if (need * sizeof (Unit) > sizeof (resbuff))
602
            {                   /* too big for local */
603
              allocres = (Unit *) malloc (need * sizeof (Unit));
604
              if (allocres == NULL)
605
                {
606
                  status |= DEC_Insufficient_storage;
607
                  break;
608
                }
609
              res = allocres;
610
            }
611
        }
612
      /* res now -> number lsu, buffer, or allocated storage for Unit array */
613
 
614
      /* Place the coefficient into the selected Unit array */
615
#if DECDPUN>1
616
      i = d % DECDPUN;          /* digits in top unit */
617
      if (i == 0)
618
        i = DECDPUN;
619
      up = res + D2U (d) - 1;   /* -> msu */
620
      *up = 0;
621
      for (c = cfirst;; c++)
622
        {                       /* along the digits */
623
          if (*c == '.')
624
            {                   /* ignore . [don't decrement i] */
625
              if (c != last)
626
                continue;
627
              break;
628
            }
629
          *up = (Unit) (X10 (*up) + (Int) * c - (Int) '0');
630
          i--;
631
          if (i > 0)
632
            continue;           /* more for this unit */
633
          if (up == res)
634
            break;              /* just filled the last unit */
635
          i = DECDPUN;
636
          up--;
637
          *up = 0;
638
        }                       /* c */
639
#else
640
      /* DECDPUN==1 */
641
      up = res;                 /* -> lsu */
642
      for (c = last; c >= cfirst; c--)
643
        {                       /* over each character, from least */
644
          if (*c == '.')
645
            continue;           /* ignore . [don't step b] */
646
          *up = (Unit) ((Int) * c - (Int) '0');
647
          up++;
648
        }                       /* c */
649
#endif
650
 
651
      dn->bits = bits;
652
      dn->exponent = exponent;
653
      dn->digits = d;
654
 
655
      /* if not in number (too long) shorten into the number */
656
      if (d > set->digits)
657
        decSetCoeff (dn, set, res, d, &residue, &status);
658
 
659
      /* Finally check for overflow or subnormal and round as needed */
660
      decFinalize (dn, set, &residue, &status);
661
      /* decNumberShow(dn); */
662
    }
663
  while (0);                     /* [for break] */
664
 
665
  if (allocres != NULL)
666
    free (allocres);            /* drop any storage we used */
667
  if (status != 0)
668
    decStatus (dn, status, set);
669
  return dn;
670
}
671
 
672
/* ================================================================== */
673
/* Operators                                                          */
674
/* ================================================================== */
675
 
676
/* ------------------------------------------------------------------ */
677
/* decNumberAbs -- absolute value operator                            */
678
/*                                                                    */
679
/*   This computes C = abs(A)                                         */
680
/*                                                                    */
681
/*   res is C, the result.  C may be A                                */
682
/*   rhs is A                                                         */
683
/*   set is the context                                               */
684
/*                                                                    */
685
/* C must have space for set->digits digits.                          */
686
/* ------------------------------------------------------------------ */
687
/* This has the same effect as decNumberPlus unless A is negative,    */
688
/* in which case it has the same effect as decNumberMinus.            */
689
/* ------------------------------------------------------------------ */
690
decNumber *
691
decNumberAbs (decNumber * res, const decNumber * rhs, decContext * set)
692
{
693
  decNumber dzero;              /* for 0 */
694
  uInt status = 0;               /* accumulator */
695
 
696
#if DECCHECK
697
  if (decCheckOperands (res, DECUNUSED, rhs, set))
698
    return res;
699
#endif
700
 
701
  decNumberZero (&dzero);       /* set 0 */
702
  dzero.exponent = rhs->exponent;       /* [no coefficient expansion] */
703
  decAddOp (res, &dzero, rhs, set, (uByte) (rhs->bits & DECNEG), &status);
704
  if (status != 0)
705
    decStatus (res, status, set);
706
  return res;
707
}
708
 
709
/* ------------------------------------------------------------------ */
710
/* decNumberAdd -- add two Numbers                                    */
711
/*                                                                    */
712
/*   This computes C = A + B                                          */
713
/*                                                                    */
714
/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
715
/*   lhs is A                                                         */
716
/*   rhs is B                                                         */
717
/*   set is the context                                               */
718
/*                                                                    */
719
/* C must have space for set->digits digits.                          */
720
/* ------------------------------------------------------------------ */
721
/* This just calls the routine shared with Subtract                   */
722
decNumber *
723
decNumberAdd (decNumber * res, const decNumber * lhs,
724
              const decNumber * rhs, decContext * set)
725
{
726
  uInt status = 0;               /* accumulator */
727
  decAddOp (res, lhs, rhs, set, 0, &status);
728
  if (status != 0)
729
    decStatus (res, status, set);
730
  return res;
731
}
732
 
733
/* ------------------------------------------------------------------ */
734
/* decNumberCompare -- compare two Numbers                            */
735
/*                                                                    */
736
/*   This computes C = A ? B                                          */
737
/*                                                                    */
738
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
739
/*   lhs is A                                                         */
740
/*   rhs is B                                                         */
741
/*   set is the context                                               */
742
/*                                                                    */
743
/* C must have space for one digit.                                   */
744
/* ------------------------------------------------------------------ */
745
decNumber *
746
decNumberCompare (decNumber * res, const decNumber * lhs,
747
                  const decNumber * rhs, decContext * set)
748
{
749
  uInt status = 0;               /* accumulator */
750
  decCompareOp (res, lhs, rhs, set, COMPARE, &status);
751
  if (status != 0)
752
    decStatus (res, status, set);
753
  return res;
754
}
755
 
756
/* ------------------------------------------------------------------ */
757
/* decNumberDivide -- divide one number by another                    */
758
/*                                                                    */
759
/*   This computes C = A / B                                          */
760
/*                                                                    */
761
/*   res is C, the result.  C may be A and/or B (e.g., X=X/X)         */
762
/*   lhs is A                                                         */
763
/*   rhs is B                                                         */
764
/*   set is the context                                               */
765
/*                                                                    */
766
/* C must have space for set->digits digits.                          */
767
/* ------------------------------------------------------------------ */
768
decNumber *
769
decNumberDivide (decNumber * res, const decNumber * lhs,
770
                 const decNumber * rhs, decContext * set)
771
{
772
  uInt status = 0;               /* accumulator */
773
  decDivideOp (res, lhs, rhs, set, DIVIDE, &status);
774
  if (status != 0)
775
    decStatus (res, status, set);
776
  return res;
777
}
778
 
779
/* ------------------------------------------------------------------ */
780
/* decNumberDivideInteger -- divide and return integer quotient       */
781
/*                                                                    */
782
/*   This computes C = A # B, where # is the integer divide operator  */
783
/*                                                                    */
784
/*   res is C, the result.  C may be A and/or B (e.g., X=X#X)         */
785
/*   lhs is A                                                         */
786
/*   rhs is B                                                         */
787
/*   set is the context                                               */
788
/*                                                                    */
789
/* C must have space for set->digits digits.                          */
790
/* ------------------------------------------------------------------ */
791
decNumber *
792
decNumberDivideInteger (decNumber * res, const decNumber * lhs,
793
                        const decNumber * rhs, decContext * set)
794
{
795
  uInt status = 0;               /* accumulator */
796
  decDivideOp (res, lhs, rhs, set, DIVIDEINT, &status);
797
  if (status != 0)
798
    decStatus (res, status, set);
799
  return res;
800
}
801
 
802
/* ------------------------------------------------------------------ */
803
/* decNumberMax -- compare two Numbers and return the maximum         */
804
/*                                                                    */
805
/*   This computes C = A ? B, returning the maximum or A if equal     */
806
/*                                                                    */
807
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
808
/*   lhs is A                                                         */
809
/*   rhs is B                                                         */
810
/*   set is the context                                               */
811
/*                                                                    */
812
/* C must have space for set->digits digits.                          */
813
/* ------------------------------------------------------------------ */
814
decNumber *
815
decNumberMax (decNumber * res, const decNumber * lhs,
816
              const decNumber * rhs, decContext * set)
817
{
818
  uInt status = 0;               /* accumulator */
819
  decCompareOp (res, lhs, rhs, set, COMPMAX, &status);
820
  if (status != 0)
821
    decStatus (res, status, set);
822
  return res;
823
}
824
 
825
/* ------------------------------------------------------------------ */
826
/* decNumberMin -- compare two Numbers and return the minimum         */
827
/*                                                                    */
828
/*   This computes C = A ? B, returning the minimum or A if equal     */
829
/*                                                                    */
830
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
831
/*   lhs is A                                                         */
832
/*   rhs is B                                                         */
833
/*   set is the context                                               */
834
/*                                                                    */
835
/* C must have space for set->digits digits.                          */
836
/* ------------------------------------------------------------------ */
837
decNumber *
838
decNumberMin (decNumber * res, const decNumber * lhs,
839
              const decNumber * rhs, decContext * set)
840
{
841
  uInt status = 0;               /* accumulator */
842
  decCompareOp (res, lhs, rhs, set, COMPMIN, &status);
843
  if (status != 0)
844
    decStatus (res, status, set);
845
  return res;
846
}
847
 
848
/* ------------------------------------------------------------------ */
849
/* decNumberMinus -- prefix minus operator                            */
850
/*                                                                    */
851
/*   This computes C = 0 - A                                          */
852
/*                                                                    */
853
/*   res is C, the result.  C may be A                                */
854
/*   rhs is A                                                         */
855
/*   set is the context                                               */
856
/*                                                                    */
857
/* C must have space for set->digits digits.                          */
858
/* ------------------------------------------------------------------ */
859
/* We simply use AddOp for the subtract, which will do the necessary. */
860
/* ------------------------------------------------------------------ */
861
decNumber *
862
decNumberMinus (decNumber * res, const decNumber * rhs, decContext * set)
863
{
864
  decNumber dzero;
865
  uInt status = 0;               /* accumulator */
866
 
867
#if DECCHECK
868
  if (decCheckOperands (res, DECUNUSED, rhs, set))
869
    return res;
870
#endif
871
 
872
  decNumberZero (&dzero);       /* make 0 */
873
  dzero.exponent = rhs->exponent;       /* [no coefficient expansion] */
874
  decAddOp (res, &dzero, rhs, set, DECNEG, &status);
875
  if (status != 0)
876
    decStatus (res, status, set);
877
  return res;
878
}
879
 
880
/* ------------------------------------------------------------------ */
881
/* decNumberPlus -- prefix plus operator                              */
882
/*                                                                    */
883
/*   This computes C = 0 + A                                          */
884
/*                                                                    */
885
/*   res is C, the result.  C may be A                                */
886
/*   rhs is A                                                         */
887
/*   set is the context                                               */
888
/*                                                                    */
889
/* C must have space for set->digits digits.                          */
890
/* ------------------------------------------------------------------ */
891
/* We simply use AddOp; Add will take fast path after preparing A.    */
892
/* Performance is a concern here, as this routine is often used to    */
893
/* check operands and apply rounding and overflow/underflow testing.  */
894
/* ------------------------------------------------------------------ */
895
decNumber *
896
decNumberPlus (decNumber * res, const decNumber * rhs, decContext * set)
897
{
898
  decNumber dzero;
899
  uInt status = 0;               /* accumulator */
900
 
901
#if DECCHECK
902
  if (decCheckOperands (res, DECUNUSED, rhs, set))
903
    return res;
904
#endif
905
 
906
  decNumberZero (&dzero);       /* make 0 */
907
  dzero.exponent = rhs->exponent;       /* [no coefficient expansion] */
908
  decAddOp (res, &dzero, rhs, set, 0, &status);
909
  if (status != 0)
910
    decStatus (res, status, set);
911
  return res;
912
}
913
 
914
/* ------------------------------------------------------------------ */
915
/* decNumberMultiply -- multiply two Numbers                          */
916
/*                                                                    */
917
/*   This computes C = A x B                                          */
918
/*                                                                    */
919
/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
920
/*   lhs is A                                                         */
921
/*   rhs is B                                                         */
922
/*   set is the context                                               */
923
/*                                                                    */
924
/* C must have space for set->digits digits.                          */
925
/* ------------------------------------------------------------------ */
926
decNumber *
927
decNumberMultiply (decNumber * res, const decNumber * lhs,
928
                   const decNumber * rhs, decContext * set)
929
{
930
  uInt status = 0;               /* accumulator */
931
  decMultiplyOp (res, lhs, rhs, set, &status);
932
  if (status != 0)
933
    decStatus (res, status, set);
934
  return res;
935
}
936
 
937
/* ------------------------------------------------------------------ */
938
/* decNumberNormalize -- remove trailing zeros                        */
939
/*                                                                    */
940
/*   This computes C = 0 + A, and normalizes the result               */
941
/*                                                                    */
942
/*   res is C, the result.  C may be A                                */
943
/*   rhs is A                                                         */
944
/*   set is the context                                               */
945
/*                                                                    */
946
/* C must have space for set->digits digits.                          */
947
/* ------------------------------------------------------------------ */
948
decNumber *
949
decNumberNormalize (decNumber * res, const decNumber * rhs, decContext * set)
950
{
951
  decNumber *allocrhs = NULL;   /* non-NULL if rounded rhs allocated */
952
  uInt status = 0;               /* as usual */
953
  Int residue = 0;               /* as usual */
954
  Int dropped;                  /* work */
955
 
956
#if DECCHECK
957
  if (decCheckOperands (res, DECUNUSED, rhs, set))
958
    return res;
959
#endif
960
 
961
  do
962
    {                           /* protect allocated storage */
963
#if DECSUBSET
964
      if (!set->extended)
965
        {
966
          /* reduce operand and set lostDigits status, as needed */
967
          if (rhs->digits > set->digits)
968
            {
969
              allocrhs = decRoundOperand (rhs, set, &status);
970
              if (allocrhs == NULL)
971
                break;
972
              rhs = allocrhs;
973
            }
974
        }
975
#endif
976
      /* [following code does not require input rounding] */
977
 
978
      /* specials copy through, except NaNs need care */
979
      if (decNumberIsNaN (rhs))
980
        {
981
          decNaNs (res, rhs, NULL, &status);
982
          break;
983
        }
984
 
985
      /* reduce result to the requested length and copy to result */
986
      decCopyFit (res, rhs, set, &residue, &status);    /* copy & round */
987
      decFinish (res, set, &residue, &status);  /* cleanup/set flags */
988
      decTrim (res, 1, &dropped);       /* normalize in place */
989
    }
990
  while (0);                     /* end protected */
991
 
992
  if (allocrhs != NULL)
993
    free (allocrhs);            /* .. */
994
  if (status != 0)
995
    decStatus (res, status, set);       /* then report status */
996
  return res;
997
}
998
 
999
/* ------------------------------------------------------------------ */
1000
/* decNumberPower -- raise a number to an integer power               */
1001
/*                                                                    */
1002
/*   This computes C = A ** B                                         */
1003
/*                                                                    */
1004
/*   res is C, the result.  C may be A and/or B (e.g., X=X**X)        */
1005
/*   lhs is A                                                         */
1006
/*   rhs is B                                                         */
1007
/*   set is the context                                               */
1008
/*                                                                    */
1009
/* C must have space for set->digits digits.                          */
1010
/*                                                                    */
1011
/* Specification restriction: abs(n) must be <=999999999              */
1012
/* ------------------------------------------------------------------ */
1013
decNumber *
1014
decNumberPower (decNumber * res, const decNumber * lhs,
1015
                const decNumber * rhs, decContext * set)
1016
{
1017
  decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
1018
  decNumber *allocrhs = NULL;   /* .., rhs */
1019
  decNumber *allocdac = NULL;   /* -> allocated acc buffer, iff used */
1020
  const decNumber *inrhs = rhs; /* save original rhs */
1021
  Int reqdigits = set->digits;  /* requested DIGITS */
1022
  Int n;                        /* RHS in binary */
1023
  Int i;                        /* work */
1024
#if DECSUBSET
1025
  Int dropped;                  /* .. */
1026
#endif
1027
  uInt needbytes;               /* buffer size needed */
1028
  Flag seenbit;                 /* seen a bit while powering */
1029
  Int residue = 0;               /* rounding residue */
1030
  uInt status = 0;               /* accumulator */
1031
  uByte bits = 0;                /* result sign if errors */
1032
  decContext workset;           /* working context */
1033
  decNumber dnOne;              /* work value 1... */
1034
  /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
1035
  uByte dacbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
1036
  /* same again for possible 1/lhs calculation */
1037
  uByte lhsbuff[sizeof (decNumber) + D2U (DECBUFFER + 9) * sizeof (Unit)];
1038
  decNumber *dac = (decNumber *) dacbuff;       /* -> result accumulator */
1039
 
1040
#if DECCHECK
1041
  if (decCheckOperands (res, lhs, rhs, set))
1042
    return res;
1043
#endif
1044
 
1045
  do
1046
    {                           /* protect allocated storage */
1047
#if DECSUBSET
1048
      if (!set->extended)
1049
        {
1050
          /* reduce operands and set lostDigits status, as needed */
1051
          if (lhs->digits > reqdigits)
1052
            {
1053
              alloclhs = decRoundOperand (lhs, set, &status);
1054
              if (alloclhs == NULL)
1055
                break;
1056
              lhs = alloclhs;
1057
            }
1058
          /* rounding won't affect the result, but we might signal lostDigits */
1059
          /* as well as the error for non-integer [x**y would need this too] */
1060
          if (rhs->digits > reqdigits)
1061
            {
1062
              allocrhs = decRoundOperand (rhs, set, &status);
1063
              if (allocrhs == NULL)
1064
                break;
1065
              rhs = allocrhs;
1066
            }
1067
        }
1068
#endif
1069
      /* [following code does not require input rounding] */
1070
 
1071
      /* handle rhs Infinity */
1072
      if (decNumberIsInfinite (rhs))
1073
        {
1074
          status |= DEC_Invalid_operation;      /* bad */
1075
          break;
1076
        }
1077
      /* handle NaNs */
1078
      if ((lhs->bits | rhs->bits) & (DECNAN | DECSNAN))
1079
        {
1080
          decNaNs (res, lhs, rhs, &status);
1081
          break;
1082
        }
1083
 
1084
      /* Original rhs must be an integer that fits and is in range */
1085
#if DECSUBSET
1086
      n = decGetInt (inrhs, set);
1087
#else
1088
      n = decGetInt (inrhs);
1089
#endif
1090
      if (n == BADINT || n > 999999999 || n < -999999999)
1091
        {
1092
          status |= DEC_Invalid_operation;
1093
          break;
1094
        }
1095
      if (n < 0)
1096
        {                       /* negative */
1097
          n = -n;               /* use the absolute value */
1098
        }
1099
      if (decNumberIsNegative (lhs)     /* -x .. */
1100
          && (n & 0x00000001))
1101
        bits = DECNEG;          /* .. to an odd power */
1102
 
1103
      /* handle LHS infinity */
1104
      if (decNumberIsInfinite (lhs))
1105
        {                       /* [NaNs already handled] */
1106
          uByte rbits = rhs->bits;      /* save */
1107
          decNumberZero (res);
1108
          if (n == 0)
1109
            *res->lsu = 1;      /* [-]Inf**0 => 1 */
1110
          else
1111
            {
1112
              if (!(rbits & DECNEG))
1113
                bits |= DECINF; /* was not a **-n */
1114
              /* [otherwise will be 0 or -0] */
1115
              res->bits = bits;
1116
            }
1117
          break;
1118
        }
1119
 
1120
      /* clone the context */
1121
      workset = *set;           /* copy all fields */
1122
      /* calculate the working DIGITS */
1123
      workset.digits = reqdigits + (inrhs->digits + inrhs->exponent) + 1;
1124
      /* it's an error if this is more than we can handle */
1125
      if (workset.digits > DECNUMMAXP)
1126
        {
1127
          status |= DEC_Invalid_operation;
1128
          break;
1129
        }
1130
 
1131
      /* workset.digits is the count of digits for the accumulator we need */
1132
      /* if accumulator is too long for local storage, then allocate */
1133
      needbytes =
1134
        sizeof (decNumber) + (D2U (workset.digits) - 1) * sizeof (Unit);
1135
      /* [needbytes also used below if 1/lhs needed] */
1136
      if (needbytes > sizeof (dacbuff))
1137
        {
1138
          allocdac = (decNumber *) malloc (needbytes);
1139
          if (allocdac == NULL)
1140
            {                   /* hopeless -- abandon */
1141
              status |= DEC_Insufficient_storage;
1142
              break;
1143
            }
1144
          dac = allocdac;       /* use the allocated space */
1145
        }
1146
      decNumberZero (dac);      /* acc=1 */
1147
      *dac->lsu = 1;            /* .. */
1148
 
1149
      if (n == 0)
1150
        {                       /* x**0 is usually 1 */
1151
          /* 0**0 is bad unless subset, when it becomes 1 */
1152
          if (ISZERO (lhs)
1153
#if DECSUBSET
1154
              && set->extended
1155
#endif
1156
            )
1157
            status |= DEC_Invalid_operation;
1158
          else
1159
            decNumberCopy (res, dac);   /* copy the 1 */
1160
          break;
1161
        }
1162
 
1163
      /* if a negative power we'll need the constant 1, and if not subset */
1164
      /* we'll invert the lhs now rather than inverting the result later */
1165
      if (decNumberIsNegative (rhs))
1166
        {                       /* was a **-n [hence digits>0] */
1167
          decNumber * newlhs;
1168
          decNumberCopy (&dnOne, dac);  /* dnOne=1;  [needed now or later] */
1169
#if DECSUBSET
1170
          if (set->extended)
1171
            {                   /* need to calculate 1/lhs */
1172
#endif
1173
              /* divide lhs into 1, putting result in dac [dac=1/dac] */
1174
              decDivideOp (dac, &dnOne, lhs, &workset, DIVIDE, &status);
1175
              if (alloclhs != NULL)
1176
                {
1177
                  free (alloclhs);      /* done with intermediate */
1178
                  alloclhs = NULL;      /* indicate freed */
1179
                }
1180
              /* now locate or allocate space for the inverted lhs */
1181
              if (needbytes > sizeof (lhsbuff))
1182
                {
1183
                  alloclhs = (decNumber *) malloc (needbytes);
1184
                  if (alloclhs == NULL)
1185
                    {           /* hopeless -- abandon */
1186
                      status |= DEC_Insufficient_storage;
1187
                      break;
1188
                    }
1189
                  newlhs = alloclhs;    /* use the allocated space */
1190
                }
1191
              else
1192
                newlhs = (decNumber *) lhsbuff; /* use stack storage */
1193
              /* [lhs now points to buffer or allocated storage] */
1194
              decNumberCopy (newlhs, dac);      /* copy the 1/lhs */
1195
              decNumberCopy (dac, &dnOne);      /* restore acc=1 */
1196
              lhs = newlhs;
1197
#if DECSUBSET
1198
            }
1199
#endif
1200
        }
1201
 
1202
      /* Raise-to-the-power loop... */
1203
      seenbit = 0;               /* set once we've seen a 1-bit */
1204
      for (i = 1;; i++)
1205
        {                       /* for each bit [top bit ignored] */
1206
          /* abandon if we have had overflow or terminal underflow */
1207
          if (status & (DEC_Overflow | DEC_Underflow))
1208
            {                   /* interesting? */
1209
              if (status & DEC_Overflow || ISZERO (dac))
1210
                break;
1211
            }
1212
          /* [the following two lines revealed an optimizer bug in a C++ */
1213
          /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
1214
          n = n << 1;           /* move next bit to testable position */
1215
          if (n < 0)
1216
            {                   /* top bit is set */
1217
              seenbit = 1;      /* OK, we're off */
1218
              decMultiplyOp (dac, dac, lhs, &workset, &status); /* dac=dac*x */
1219
            }
1220
          if (i == 31)
1221
            break;              /* that was the last bit */
1222
          if (!seenbit)
1223
            continue;           /* we don't have to square 1 */
1224
          decMultiplyOp (dac, dac, dac, &workset, &status);     /* dac=dac*dac [square] */
1225
        }                       /*i *//* 32 bits */
1226
 
1227
      /* complete internal overflow or underflow processing */
1228
      if (status & (DEC_Overflow | DEC_Subnormal))
1229
        {
1230
#if DECSUBSET
1231
          /* If subset, and power was negative, reverse the kind of -erflow */
1232
          /* [1/x not yet done] */
1233
          if (!set->extended && decNumberIsNegative (rhs))
1234
            {
1235
              if (status & DEC_Overflow)
1236
                status ^= DEC_Overflow | DEC_Underflow | DEC_Subnormal;
1237
              else
1238
                {               /* trickier -- Underflow may or may not be set */
1239
                  status &= ~(DEC_Underflow | DEC_Subnormal);   /* [one or both] */
1240
                  status |= DEC_Overflow;
1241
                }
1242
            }
1243
#endif
1244
          dac->bits = (dac->bits & ~DECNEG) | bits;     /* force correct sign */
1245
          /* round subnormals [to set.digits rather than workset.digits] */
1246
          /* or set overflow result similarly as required */
1247
          decFinalize (dac, set, &residue, &status);
1248
          decNumberCopy (res, dac);     /* copy to result (is now OK length) */
1249
          break;
1250
        }
1251
 
1252
#if DECSUBSET
1253
      if (!set->extended &&     /* subset math */
1254
          decNumberIsNegative (rhs))
1255
        {                       /* was a **-n [hence digits>0] */
1256
          /* so divide result into 1 [dac=1/dac] */
1257
          decDivideOp (dac, &dnOne, dac, &workset, DIVIDE, &status);
1258
        }
1259
#endif
1260
 
1261
      /* reduce result to the requested length and copy to result */
1262
      decCopyFit (res, dac, set, &residue, &status);
1263
      decFinish (res, set, &residue, &status);  /* final cleanup */
1264
#if DECSUBSET
1265
      if (!set->extended)
1266
        decTrim (res, 0, &dropped);      /* trailing zeros */
1267
#endif
1268
    }
1269
  while (0);                     /* end protected */
1270
 
1271
  if (allocdac != NULL)
1272
    free (allocdac);            /* drop any storage we used */
1273
  if (allocrhs != NULL)
1274
    free (allocrhs);            /* .. */
1275
  if (alloclhs != NULL)
1276
    free (alloclhs);            /* .. */
1277
  if (status != 0)
1278
    decStatus (res, status, set);
1279
  return res;
1280
}
1281
 
1282
/* ------------------------------------------------------------------ */
1283
/* decNumberQuantize -- force exponent to requested value             */
1284
/*                                                                    */
1285
/*   This computes C = op(A, B), where op adjusts the coefficient     */
1286
/*   of C (by rounding or shifting) such that the exponent (-scale)   */
1287
/*   of C has exponent of B.  The numerical value of C will equal A,  */
1288
/*   except for the effects of any rounding that occurred.            */
1289
/*                                                                    */
1290
/*   res is C, the result.  C may be A or B                           */
1291
/*   lhs is A, the number to adjust                                   */
1292
/*   rhs is B, the number with exponent to match                      */
1293
/*   set is the context                                               */
1294
/*                                                                    */
1295
/* C must have space for set->digits digits.                          */
1296
/*                                                                    */
1297
/* Unless there is an error or the result is infinite, the exponent   */
1298
/* after the operation is guaranteed to be equal to that of B.        */
1299
/* ------------------------------------------------------------------ */
1300
decNumber *
1301
decNumberQuantize (decNumber * res, const decNumber * lhs,
1302
                   const decNumber * rhs, decContext * set)
1303
{
1304
  uInt status = 0;               /* accumulator */
1305
  decQuantizeOp (res, lhs, rhs, set, 1, &status);
1306
  if (status != 0)
1307
    decStatus (res, status, set);
1308
  return res;
1309
}
1310
 
1311
/* ------------------------------------------------------------------ */
1312
/* decNumberRescale -- force exponent to requested value              */
1313
/*                                                                    */
1314
/*   This computes C = op(A, B), where op adjusts the coefficient     */
1315
/*   of C (by rounding or shifting) such that the exponent (-scale)   */
1316
/*   of C has the value B.  The numerical value of C will equal A,    */
1317
/*   except for the effects of any rounding that occurred.            */
1318
/*                                                                    */
1319
/*   res is C, the result.  C may be A or B                           */
1320
/*   lhs is A, the number to adjust                                   */
1321
/*   rhs is B, the requested exponent                                 */
1322
/*   set is the context                                               */
1323
/*                                                                    */
1324
/* C must have space for set->digits digits.                          */
1325
/*                                                                    */
1326
/* Unless there is an error or the result is infinite, the exponent   */
1327
/* after the operation is guaranteed to be equal to B.                */
1328
/* ------------------------------------------------------------------ */
1329
decNumber *
1330
decNumberRescale (decNumber * res, const decNumber * lhs,
1331
                  const decNumber * rhs, decContext * set)
1332
{
1333
  uInt status = 0;               /* accumulator */
1334
  decQuantizeOp (res, lhs, rhs, set, 0, &status);
1335
  if (status != 0)
1336
    decStatus (res, status, set);
1337
  return res;
1338
}
1339
 
1340
/* ------------------------------------------------------------------ */
1341
/* decNumberRemainder -- divide and return remainder                  */
1342
/*                                                                    */
1343
/*   This computes C = A % B                                          */
1344
/*                                                                    */
1345
/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
1346
/*   lhs is A                                                         */
1347
/*   rhs is B                                                         */
1348
/*   set is the context                                               */
1349
/*                                                                    */
1350
/* C must have space for set->digits digits.                          */
1351
/* ------------------------------------------------------------------ */
1352
decNumber *
1353
decNumberRemainder (decNumber * res, const decNumber * lhs,
1354
                    const decNumber * rhs, decContext * set)
1355
{
1356
  uInt status = 0;               /* accumulator */
1357
  decDivideOp (res, lhs, rhs, set, REMAINDER, &status);
1358
  if (status != 0)
1359
    decStatus (res, status, set);
1360
  return res;
1361
}
1362
 
1363
/* ------------------------------------------------------------------ */
1364
/* decNumberRemainderNear -- divide and return remainder from nearest */
1365
/*                                                                    */
1366
/*   This computes C = A % B, where % is the IEEE remainder operator  */
1367
/*                                                                    */
1368
/*   res is C, the result.  C may be A and/or B (e.g., X=X%X)         */
1369
/*   lhs is A                                                         */
1370
/*   rhs is B                                                         */
1371
/*   set is the context                                               */
1372
/*                                                                    */
1373
/* C must have space for set->digits digits.                          */
1374
/* ------------------------------------------------------------------ */
1375
decNumber *
1376
decNumberRemainderNear (decNumber * res, const decNumber * lhs,
1377
                        const decNumber * rhs, decContext * set)
1378
{
1379
  uInt status = 0;               /* accumulator */
1380
  decDivideOp (res, lhs, rhs, set, REMNEAR, &status);
1381
  if (status != 0)
1382
    decStatus (res, status, set);
1383
  return res;
1384
}
1385
 
1386
/* ------------------------------------------------------------------ */
1387
/* decNumberSameQuantum -- test for equal exponents                   */
1388
/*                                                                    */
1389
/*   res is the result number, which will contain either 0 or 1       */
1390
/*   lhs is a number to test                                          */
1391
/*   rhs is the second (usually a pattern)                            */
1392
/*                                                                    */
1393
/* No errors are possible and no context is needed.                   */
1394
/* ------------------------------------------------------------------ */
1395
decNumber *
1396
decNumberSameQuantum (decNumber * res, const decNumber * lhs, const decNumber * rhs)
1397
{
1398
  uByte merged;                 /* merged flags */
1399
  Unit ret = 0;                  /* return value */
1400
 
1401
#if DECCHECK
1402
  if (decCheckOperands (res, lhs, rhs, DECUNUSED))
1403
    return res;
1404
#endif
1405
 
1406
  merged = (lhs->bits | rhs->bits) & DECSPECIAL;
1407
  if (merged)
1408
    {
1409
      if (decNumberIsNaN (lhs) && decNumberIsNaN (rhs))
1410
        ret = 1;
1411
      else if (decNumberIsInfinite (lhs) && decNumberIsInfinite (rhs))
1412
        ret = 1;
1413
      /* [anything else with a special gives 0] */
1414
    }
1415
  else if (lhs->exponent == rhs->exponent)
1416
    ret = 1;
1417
 
1418
  decNumberZero (res);          /* OK to overwrite an operand */
1419
  *res->lsu = ret;
1420
  return res;
1421
}
1422
 
1423
/* ------------------------------------------------------------------ */
1424
/* decNumberSquareRoot -- square root operator                        */
1425
/*                                                                    */
1426
/*   This computes C = squareroot(A)                                  */
1427
/*                                                                    */
1428
/*   res is C, the result.  C may be A                                */
1429
/*   rhs is A                                                         */
1430
/*   set is the context; note that rounding mode has no effect        */
1431
/*                                                                    */
1432
/* C must have space for set->digits digits.                          */
1433
/* ------------------------------------------------------------------ */
1434
/* This uses the following varying-precision algorithm in:            */
1435
/*                                                                    */
1436
/*   Properly Rounded Variable Precision Square Root, T. E. Hull and  */
1437
/*   A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
1438
/*   pp229-237, ACM, September 1985.                                  */
1439
/*                                                                    */
1440
/* % [Reformatted original Numerical Turing source code follows.]     */
1441
/* function sqrt(x : real) : real                                     */
1442
/* % sqrt(x) returns the properly rounded approximation to the square */
1443
/* % root of x, in the precision of the calling environment, or it    */
1444
/* % fails if x < 0.                                                  */
1445
/* % t e hull and a abrham, august, 1984                              */
1446
/* if x <= 0 then                                                     */
1447
/*   if x < 0 then                                                    */
1448
/*     assert false                                                   */
1449
/*   else                                                             */
1450
/*     result 0                                                       */
1451
/*   end if                                                           */
1452
/* end if                                                             */
1453
/* var f := setexp(x, 0)  % fraction part of x   [0.1 <= x < 1]       */
1454
/* var e := getexp(x)     % exponent part of x                        */
1455
/* var approx : real                                                  */
1456
/* if e mod 2 = 0  then                                               */
1457
/*   approx := .259 + .819 * f   % approx to root of f                */
1458
/* else                                                               */
1459
/*   f := f/l0                   % adjustments                        */
1460
/*   e := e + 1                  %   for odd                          */
1461
/*   approx := .0819 + 2.59 * f  %   exponent                         */
1462
/* end if                                                             */
1463
/*                                                                    */
1464
/* var p:= 3                                                          */
1465
/* const maxp := currentprecision + 2                                 */
1466
/* loop                                                               */
1467
/*   p := min(2*p - 2, maxp)     % p = 4,6,10, . . . , maxp           */
1468
/*   precision p                                                      */
1469
/*   approx := .5 * (approx + f/approx)                               */
1470
/*   exit when p = maxp                                               */
1471
/* end loop                                                           */
1472
/*                                                                    */
1473
/* % approx is now within 1 ulp of the properly rounded square root   */
1474
/* % of f; to ensure proper rounding, compare squares of (approx -    */
1475
/* % l/2 ulp) and (approx + l/2 ulp) with f.                          */
1476
/* p := currentprecision                                              */
1477
/* begin                                                              */
1478
/*   precision p + 2                                                  */
1479
/*   const approxsubhalf := approx - setexp(.5, -p)                   */
1480
/*   if mulru(approxsubhalf, approxsubhalf) > f then                  */
1481
/*     approx := approx - setexp(.l, -p + 1)                          */
1482
/*   else                                                             */
1483
/*     const approxaddhalf := approx + setexp(.5, -p)                 */
1484
/*     if mulrd(approxaddhalf, approxaddhalf) < f then                */
1485
/*       approx := approx + setexp(.l, -p + 1)                        */
1486
/*     end if                                                         */
1487
/*   end if                                                           */
1488
/* end                                                                */
1489
/* result setexp(approx, e div 2)  % fix exponent                     */
1490
/* end sqrt                                                           */
1491
/* ------------------------------------------------------------------ */
1492
decNumber *
1493
decNumberSquareRoot (decNumber * res, const decNumber * rhs, decContext * set)
1494
{
1495
  decContext workset, approxset;        /* work contexts */
1496
  decNumber dzero;              /* used for constant zero */
1497
  Int maxp = set->digits + 2;   /* largest working precision */
1498
  Int residue = 0;               /* rounding residue */
1499
  uInt status = 0, ignore = 0;    /* status accumulators */
1500
  Int exp;                      /* working exponent */
1501
  Int ideal;                    /* ideal (preferred) exponent */
1502
  uInt needbytes;               /* work */
1503
  Int dropped;                  /* .. */
1504
 
1505
  decNumber *allocrhs = NULL;   /* non-NULL if rounded rhs allocated */
1506
  /* buffer for f [needs +1 in case DECBUFFER 0] */
1507
  uByte buff[sizeof (decNumber) + (D2U (DECBUFFER + 1) - 1) * sizeof (Unit)];
1508
  /* buffer for a [needs +2 to match maxp] */
1509
  uByte bufa[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1510
  /* buffer for temporary, b [must be same size as a] */
1511
  uByte bufb[sizeof (decNumber) + (D2U (DECBUFFER + 2) - 1) * sizeof (Unit)];
1512
  decNumber *allocbuff = NULL;  /* -> allocated buff, iff allocated */
1513
  decNumber *allocbufa = NULL;  /* -> allocated bufa, iff allocated */
1514
  decNumber *allocbufb = NULL;  /* -> allocated bufb, iff allocated */
1515
  decNumber *f = (decNumber *) buff;    /* reduced fraction */
1516
  decNumber *a = (decNumber *) bufa;    /* approximation to result */
1517
  decNumber *b = (decNumber *) bufb;    /* intermediate result */
1518
  /* buffer for temporary variable, up to 3 digits */
1519
  uByte buft[sizeof (decNumber) + (D2U (3) - 1) * sizeof (Unit)];
1520
  decNumber *t = (decNumber *) buft;    /* up-to-3-digit constant or work */
1521
 
1522
#if DECCHECK
1523
  if (decCheckOperands (res, DECUNUSED, rhs, set))
1524
    return res;
1525
#endif
1526
 
1527
  do
1528
    {                           /* protect allocated storage */
1529
#if DECSUBSET
1530
      if (!set->extended)
1531
        {
1532
          /* reduce operand and set lostDigits status, as needed */
1533
          if (rhs->digits > set->digits)
1534
            {
1535
              allocrhs = decRoundOperand (rhs, set, &status);
1536
              if (allocrhs == NULL)
1537
                break;
1538
              /* [Note: 'f' allocation below could reuse this buffer if */
1539
              /* used, but as this is rare we keep them separate for clarity.] */
1540
              rhs = allocrhs;
1541
            }
1542
        }
1543
#endif
1544
      /* [following code does not require input rounding] */
1545
 
1546
      /* handle infinities and NaNs */
1547
      if (rhs->bits & DECSPECIAL)
1548
        {
1549
          if (decNumberIsInfinite (rhs))
1550
            {                   /* an infinity */
1551
              if (decNumberIsNegative (rhs))
1552
                status |= DEC_Invalid_operation;
1553
              else
1554
                decNumberCopy (res, rhs);       /* +Infinity */
1555
            }
1556
          else
1557
            decNaNs (res, rhs, NULL, &status);  /* a NaN */
1558
          break;
1559
        }
1560
 
1561
      /* calculate the ideal (preferred) exponent [floor(exp/2)] */
1562
      /* [We would like to write: ideal=rhs->exponent>>1, but this */
1563
      /* generates a compiler warning.  Generated code is the same.] */
1564
      ideal = (rhs->exponent & ~1) / 2; /* target */
1565
 
1566
      /* handle zeros */
1567
      if (ISZERO (rhs))
1568
        {
1569
          decNumberCopy (res, rhs);     /* could be 0 or -0 */
1570
          res->exponent = ideal;        /* use the ideal [safe] */
1571
          break;
1572
        }
1573
 
1574
      /* any other -x is an oops */
1575
      if (decNumberIsNegative (rhs))
1576
        {
1577
          status |= DEC_Invalid_operation;
1578
          break;
1579
        }
1580
 
1581
      /* we need space for three working variables */
1582
      /*   f -- the same precision as the RHS, reduced to 0.01->0.99... */
1583
      /*   a -- Hull's approx -- precision, when assigned, is */
1584
      /*        currentprecision (we allow +2 for use as temporary) */
1585
      /*   b -- intermediate temporary result */
1586
      /* if any is too long for local storage, then allocate */
1587
      needbytes =
1588
        sizeof (decNumber) + (D2U (rhs->digits) - 1) * sizeof (Unit);
1589
      if (needbytes > sizeof (buff))
1590
        {
1591
          allocbuff = (decNumber *) malloc (needbytes);
1592
          if (allocbuff == NULL)
1593
            {                   /* hopeless -- abandon */
1594
              status |= DEC_Insufficient_storage;
1595
              break;
1596
            }
1597
          f = allocbuff;        /* use the allocated space */
1598
        }
1599
      /* a and b both need to be able to hold a maxp-length number */
1600
      needbytes = sizeof (decNumber) + (D2U (maxp) - 1) * sizeof (Unit);
1601
      if (needbytes > sizeof (bufa))
1602
        {                       /* [same applies to b] */
1603
          allocbufa = (decNumber *) malloc (needbytes);
1604
          allocbufb = (decNumber *) malloc (needbytes);
1605
          if (allocbufa == NULL || allocbufb == NULL)
1606
            {                   /* hopeless */
1607
              status |= DEC_Insufficient_storage;
1608
              break;
1609
            }
1610
          a = allocbufa;        /* use the allocated space */
1611
          b = allocbufb;        /* .. */
1612
        }
1613
 
1614
      /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
1615
      decNumberCopy (f, rhs);
1616
      exp = f->exponent + f->digits;    /* adjusted to Hull rules */
1617
      f->exponent = -(f->digits);       /* to range */
1618
 
1619
      /* set up working contexts (the second is used for Numerical */
1620
      /* Turing assignment) */
1621
      decContextDefault (&workset, DEC_INIT_DECIMAL64);
1622
      decContextDefault (&approxset, DEC_INIT_DECIMAL64);
1623
      approxset.digits = set->digits;   /* approx's length */
1624
 
1625
      /* [Until further notice, no error is possible and status bits */
1626
      /* (Rounded, etc.) should be ignored, not accumulated.] */
1627
 
1628
      /* Calculate initial approximation, and allow for odd exponent */
1629
      workset.digits = set->digits;     /* p for initial calculation */
1630
      t->bits = 0;
1631
      t->digits = 3;
1632
      a->bits = 0;
1633
      a->digits = 3;
1634
      if ((exp & 1) == 0)
1635
        {                       /* even exponent */
1636
          /* Set t=0.259, a=0.819 */
1637
          t->exponent = -3;
1638
          a->exponent = -3;
1639
#if DECDPUN>=3
1640
          t->lsu[0] = 259;
1641
          a->lsu[0] = 819;
1642
#elif DECDPUN==2
1643
          t->lsu[0] = 59;
1644
          t->lsu[1] = 2;
1645
          a->lsu[0] = 19;
1646
          a->lsu[1] = 8;
1647
#else
1648
          t->lsu[0] = 9;
1649
          t->lsu[1] = 5;
1650
          t->lsu[2] = 2;
1651
          a->lsu[0] = 9;
1652
          a->lsu[1] = 1;
1653
          a->lsu[2] = 8;
1654
#endif
1655
        }
1656
      else
1657
        {                       /* odd exponent */
1658
          /* Set t=0.0819, a=2.59 */
1659
          f->exponent--;        /* f=f/10 */
1660
          exp++;                /* e=e+1 */
1661
          t->exponent = -4;
1662
          a->exponent = -2;
1663
#if DECDPUN>=3
1664
          t->lsu[0] = 819;
1665
          a->lsu[0] = 259;
1666
#elif DECDPUN==2
1667
          t->lsu[0] = 19;
1668
          t->lsu[1] = 8;
1669
          a->lsu[0] = 59;
1670
          a->lsu[1] = 2;
1671
#else
1672
          t->lsu[0] = 9;
1673
          t->lsu[1] = 1;
1674
          t->lsu[2] = 8;
1675
          a->lsu[0] = 9;
1676
          a->lsu[1] = 5;
1677
          a->lsu[2] = 2;
1678
#endif
1679
        }
1680
      decMultiplyOp (a, a, f, &workset, &ignore);       /* a=a*f */
1681
      decAddOp (a, a, t, &workset, 0, &ignore);  /* ..+t */
1682
      /* [a is now the initial approximation for sqrt(f), calculated with */
1683
      /* currentprecision, which is also a's precision.] */
1684
 
1685
      /* the main calculation loop */
1686
      decNumberZero (&dzero);   /* make 0 */
1687
      decNumberZero (t);        /* set t = 0.5 */
1688
      t->lsu[0] = 5;             /* .. */
1689
      t->exponent = -1;         /* .. */
1690
      workset.digits = 3;       /* initial p */
1691
      for (;;)
1692
        {
1693
          /* set p to min(2*p - 2, maxp)  [hence 3; or: 4, 6, 10, ... , maxp] */
1694
          workset.digits = workset.digits * 2 - 2;
1695
          if (workset.digits > maxp)
1696
            workset.digits = maxp;
1697
          /* a = 0.5 * (a + f/a) */
1698
          /* [calculated at p then rounded to currentprecision] */
1699
          decDivideOp (b, f, a, &workset, DIVIDE, &ignore);     /* b=f/a */
1700
          decAddOp (b, b, a, &workset, 0, &ignore);      /* b=b+a */
1701
          decMultiplyOp (a, b, t, &workset, &ignore);   /* a=b*0.5 */
1702
          /* assign to approx [round to length] */
1703
          decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1704
          if (workset.digits == maxp)
1705
            break;              /* just did final */
1706
        }                       /* loop */
1707
 
1708
      /* a is now at currentprecision and within 1 ulp of the properly */
1709
      /* rounded square root of f; to ensure proper rounding, compare */
1710
      /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
1711
      /* Here workset.digits=maxp and t=0.5 */
1712
      workset.digits--;         /* maxp-1 is OK now */
1713
      t->exponent = -set->digits - 1;   /* make 0.5 ulp */
1714
      decNumberCopy (b, a);
1715
      decAddOp (b, b, t, &workset, DECNEG, &ignore);    /* b = a - 0.5 ulp */
1716
      workset.round = DEC_ROUND_UP;
1717
      decMultiplyOp (b, b, b, &workset, &ignore);       /* b = mulru(b, b) */
1718
      decCompareOp (b, f, b, &workset, COMPARE, &ignore);       /* b ? f, reversed */
1719
      if (decNumberIsNegative (b))
1720
        {                       /* f < b [i.e., b > f] */
1721
          /* this is the more common adjustment, though both are rare */
1722
          t->exponent++;        /* make 1.0 ulp */
1723
          t->lsu[0] = 1; /* .. */
1724
          decAddOp (a, a, t, &workset, DECNEG, &ignore);        /* a = a - 1 ulp */
1725
          /* assign to approx [round to length] */
1726
          decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1727
        }
1728
      else
1729
        {
1730
          decNumberCopy (b, a);
1731
          decAddOp (b, b, t, &workset, 0, &ignore);      /* b = a + 0.5 ulp */
1732
          workset.round = DEC_ROUND_DOWN;
1733
          decMultiplyOp (b, b, b, &workset, &ignore);   /* b = mulrd(b, b) */
1734
          decCompareOp (b, b, f, &workset, COMPARE, &ignore);   /* b ? f */
1735
          if (decNumberIsNegative (b))
1736
            {                   /* b < f */
1737
              t->exponent++;    /* make 1.0 ulp */
1738
              t->lsu[0] = 1;     /* .. */
1739
              decAddOp (a, a, t, &workset, 0, &ignore);  /* a = a + 1 ulp */
1740
              /* assign to approx [round to length] */
1741
              decAddOp (a, &dzero, a, &approxset, 0, &ignore);
1742
            }
1743
        }
1744
      /* [no errors are possible in the above, and rounding/inexact during */
1745
      /* estimation are irrelevant, so status was not accumulated] */
1746
 
1747
      /* Here, 0.1 <= a < 1  [Hull] */
1748
      a->exponent += exp / 2;   /* set correct exponent */
1749
 
1750
      /* Process Subnormals */
1751
      decFinalize (a, set, &residue, &status);
1752
 
1753
      /* count dropable zeros [after any subnormal rounding] */
1754
      decNumberCopy (b, a);
1755
      decTrim (b, 1, &dropped); /* [drops trailing zeros] */
1756
 
1757
      /* Finally set Inexact and Rounded.  The answer can only be exact if */
1758
      /* it is short enough so that squaring it could fit in set->digits, */
1759
      /* so this is the only (relatively rare) time we have to check */
1760
      /* carefully */
1761
      if (b->digits * 2 - 1 > set->digits)
1762
        {                       /* cannot fit */
1763
          status |= DEC_Inexact | DEC_Rounded;
1764
        }
1765
      else
1766
        {                       /* could be exact/unrounded */
1767
          uInt mstatus = 0;      /* local status */
1768
          decMultiplyOp (b, b, b, &workset, &mstatus);  /* try the multiply */
1769
          if (mstatus != 0)
1770
            {                   /* result won't fit */
1771
              status |= DEC_Inexact | DEC_Rounded;
1772
            }
1773
          else
1774
            {                   /* plausible */
1775
              decCompareOp (t, b, rhs, &workset, COMPARE, &mstatus);    /* b ? rhs */
1776
              if (!ISZERO (t))
1777
                {
1778
                  status |= DEC_Inexact | DEC_Rounded;
1779
                }
1780
              else
1781
                {               /* is Exact */
1782
                  /* here, dropped is the count of trailing zeros in 'a' */
1783
                  /* use closest exponent to ideal... */
1784
                  Int todrop = ideal - a->exponent;     /* most we can drop */
1785
 
1786
                  if (todrop < 0)
1787
                    {           /* ideally would add 0s */
1788
                      status |= DEC_Rounded;
1789
                    }
1790
                  else
1791
                    {           /* unrounded */
1792
                      if (dropped < todrop)
1793
                        todrop = dropped;       /* clamp to those available */
1794
                      if (todrop > 0)
1795
                        {       /* OK, some to drop */
1796
                          decShiftToLeast (a->lsu, D2U (a->digits), todrop);
1797
                          a->exponent += todrop;        /* maintain numerical value */
1798
                          a->digits -= todrop;  /* new length */
1799
                        }
1800
                    }
1801
                }
1802
            }
1803
        }
1804
      decNumberCopy (res, a);   /* assume this is the result */
1805
    }
1806
  while (0);                     /* end protected */
1807
 
1808
  if (allocbuff != NULL)
1809
    free (allocbuff);           /* drop any storage we used */
1810
  if (allocbufa != NULL)
1811
    free (allocbufa);           /* .. */
1812
  if (allocbufb != NULL)
1813
    free (allocbufb);           /* .. */
1814
  if (allocrhs != NULL)
1815
    free (allocrhs);            /* .. */
1816
  if (status != 0)
1817
    decStatus (res, status, set);       /* then report status */
1818
  return res;
1819
}
1820
 
1821
/* ------------------------------------------------------------------ */
1822
/* decNumberSubtract -- subtract two Numbers                          */
1823
/*                                                                    */
1824
/*   This computes C = A - B                                          */
1825
/*                                                                    */
1826
/*   res is C, the result.  C may be A and/or B (e.g., X=X-X)         */
1827
/*   lhs is A                                                         */
1828
/*   rhs is B                                                         */
1829
/*   set is the context                                               */
1830
/*                                                                    */
1831
/* C must have space for set->digits digits.                          */
1832
/* ------------------------------------------------------------------ */
1833
decNumber *
1834
decNumberSubtract (decNumber * res, const decNumber * lhs,
1835
                   const decNumber * rhs, decContext * set)
1836
{
1837
  uInt status = 0;               /* accumulator */
1838
 
1839
  decAddOp (res, lhs, rhs, set, DECNEG, &status);
1840
  if (status != 0)
1841
    decStatus (res, status, set);
1842
  return res;
1843
}
1844
 
1845
/* ------------------------------------------------------------------ */
1846
/* decNumberToIntegralValue -- round-to-integral-value                */
1847
/*                                                                    */
1848
/*   res is the result                                                */
1849
/*   rhs is input number                                              */
1850
/*   set is the context                                               */
1851
/*                                                                    */
1852
/* res must have space for any value of rhs.                          */
1853
/*                                                                    */
1854
/* This implements the IEEE special operator and therefore treats     */
1855
/* special values as valid, and also never sets Inexact.  For finite  */
1856
/* numbers it returns rescale(rhs, 0) if rhs->exponent is <0.         */
1857
/* Otherwise the result is rhs (so no error is possible).             */
1858
/*                                                                    */
1859
/* The context is used for rounding mode and status after sNaN, but   */
1860
/* the digits setting is ignored.                                     */
1861
/* ------------------------------------------------------------------ */
1862
decNumber *
1863
decNumberToIntegralValue (decNumber * res, const decNumber * rhs, decContext * set)
1864
{
1865
  decNumber dn;
1866
  decContext workset;           /* working context */
1867
 
1868
#if DECCHECK
1869
  if (decCheckOperands (res, DECUNUSED, rhs, set))
1870
    return res;
1871
#endif
1872
 
1873
  /* handle infinities and NaNs */
1874
  if (rhs->bits & DECSPECIAL)
1875
    {
1876
      uInt status = 0;
1877
      if (decNumberIsInfinite (rhs))
1878
        decNumberCopy (res, rhs);       /* an Infinity */
1879
      else
1880
        decNaNs (res, rhs, NULL, &status);      /* a NaN */
1881
      if (status != 0)
1882
        decStatus (res, status, set);
1883
      return res;
1884
    }
1885
 
1886
  /* we have a finite number; no error possible */
1887
  if (rhs->exponent >= 0)
1888
    return decNumberCopy (res, rhs);
1889
  /* that was easy, but if negative exponent we have work to do... */
1890
  workset = *set;               /* clone rounding, etc. */
1891
  workset.digits = rhs->digits; /* no length rounding */
1892
  workset.traps = 0;             /* no traps */
1893
  decNumberZero (&dn);          /* make a number with exponent 0 */
1894
  return decNumberQuantize (res, rhs, &dn, &workset);
1895
}
1896
 
1897
/* ================================================================== */
1898
/* Utility routines                                                   */
1899
/* ================================================================== */
1900
 
1901
/* ------------------------------------------------------------------ */
1902
/* decNumberCopy -- copy a number                                     */
1903
/*                                                                    */
1904
/*   dest is the target decNumber                                     */
1905
/*   src  is the source decNumber                                     */
1906
/*   returns dest                                                     */
1907
/*                                                                    */
1908
/* (dest==src is allowed and is a no-op)                              */
1909
/* All fields are updated as required.  This is a utility operation,  */
1910
/* so special values are unchanged and no error is possible.          */
1911
/* ------------------------------------------------------------------ */
1912
decNumber *
1913
decNumberCopy (decNumber * dest, const decNumber * src)
1914
{
1915
 
1916
#if DECCHECK
1917
  if (src == NULL)
1918
    return decNumberZero (dest);
1919
#endif
1920
 
1921
  if (dest == src)
1922
    return dest;                /* no copy required */
1923
 
1924
  /* We use explicit assignments here as structure assignment can copy */
1925
  /* more than just the lsu (for small DECDPUN).  This would not affect */
1926
  /* the value of the results, but would disturb test harness spill */
1927
  /* checking. */
1928
  dest->bits = src->bits;
1929
  dest->exponent = src->exponent;
1930
  dest->digits = src->digits;
1931
  dest->lsu[0] = src->lsu[0];
1932
  if (src->digits > DECDPUN)
1933
    {                           /* more Units to come */
1934
      Unit *d;                  /* work */
1935
      const Unit *s, *smsup;    /* work */
1936
      /* memcpy for the remaining Units would be safe as they cannot */
1937
      /* overlap.  However, this explicit loop is faster in short cases. */
1938
      d = dest->lsu + 1;        /* -> first destination */
1939
      smsup = src->lsu + D2U (src->digits);     /* -> source msu+1 */
1940
      for (s = src->lsu + 1; s < smsup; s++, d++)
1941
        *d = *s;
1942
    }
1943
  return dest;
1944
}
1945
 
1946
/* ------------------------------------------------------------------ */
1947
/* decNumberTrim -- remove insignificant zeros                        */
1948
/*                                                                    */
1949
/*   dn is the number to trim                                         */
1950
/*   returns dn                                                       */
1951
/*                                                                    */
1952
/* All fields are updated as required.  This is a utility operation,  */
1953
/* so special values are unchanged and no error is possible.          */
1954
/* ------------------------------------------------------------------ */
1955
decNumber *
1956
decNumberTrim (decNumber * dn)
1957
{
1958
  Int dropped;                  /* work */
1959
  return decTrim (dn, 0, &dropped);
1960
}
1961
 
1962
/* ------------------------------------------------------------------ */
1963
/* decNumberVersion -- return the name and version of this module     */
1964
/*                                                                    */
1965
/* No error is possible.                                              */
1966
/* ------------------------------------------------------------------ */
1967
const char *
1968
decNumberVersion (void)
1969
{
1970
  return DECVERSION;
1971
}
1972
 
1973
/* ------------------------------------------------------------------ */
1974
/* decNumberZero -- set a number to 0                                 */
1975
/*                                                                    */
1976
/*   dn is the number to set, with space for one digit                */
1977
/*   returns dn                                                       */
1978
/*                                                                    */
1979
/* No error is possible.                                              */
1980
/* ------------------------------------------------------------------ */
1981
/* Memset is not used as it is much slower in some environments. */
1982
decNumber *
1983
decNumberZero (decNumber * dn)
1984
{
1985
 
1986
#if DECCHECK
1987
  if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
1988
    return dn;
1989
#endif
1990
 
1991
  dn->bits = 0;
1992
  dn->exponent = 0;
1993
  dn->digits = 1;
1994
  dn->lsu[0] = 0;
1995
  return dn;
1996
}
1997
 
1998
/* ================================================================== */
1999
/* Local routines                                                     */
2000
/* ================================================================== */
2001
 
2002
/* ------------------------------------------------------------------ */
2003
/* decToString -- lay out a number into a string                      */
2004
/*                                                                    */
2005
/*   dn     is the number to lay out                                  */
2006
/*   string is where to lay out the number                            */
2007
/*   eng    is 1 if Engineering, 0 if Scientific                      */
2008
/*                                                                    */
2009
/* str must be at least dn->digits+14 characters long                 */
2010
/* No error is possible.                                              */
2011
/*                                                                    */
2012
/* Note that this routine can generate a -0 or 0.000.  These are      */
2013
/* never generated in subset to-number or arithmetic, but can occur   */
2014
/* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234).              */
2015
/* ------------------------------------------------------------------ */
2016
/* If DECCHECK is enabled the string "?" is returned if a number is */
2017
/* invalid. */
2018
 
2019
/* TODIGIT -- macro to remove the leading digit from the unsigned */
2020
/* integer u at column cut (counting from the right, LSD=0) and place */
2021
/* it as an ASCII character into the character pointed to by c.  Note */
2022
/* that cut must be <= 9, and the maximum value for u is 2,000,000,000 */
2023
/* (as is needed for negative exponents of subnormals).  The unsigned */
2024
/* integer pow is used as a temporary variable. */
2025
#define TODIGIT(u, cut, c) {            \
2026
  *(c)='0';                             \
2027
  pow=powers[cut]*2;                    \
2028
  if ((u)>pow) {                        \
2029
    pow*=4;                             \
2030
    if ((u)>=pow) {(u)-=pow; *(c)+=8;}  \
2031
    pow/=2;                             \
2032
    if ((u)>=pow) {(u)-=pow; *(c)+=4;}  \
2033
    pow/=2;                             \
2034
    }                                   \
2035
  if ((u)>=pow) {(u)-=pow; *(c)+=2;}    \
2036
  pow/=2;                               \
2037
  if ((u)>=pow) {(u)-=pow; *(c)+=1;}    \
2038
  }
2039
 
2040
static void
2041
decToString (const decNumber * dn, char *string, Flag eng)
2042
{
2043
  Int exp = dn->exponent;       /* local copy */
2044
  Int e;                        /* E-part value */
2045
  Int pre;                      /* digits before the '.' */
2046
  Int cut;                      /* for counting digits in a Unit */
2047
  char *c = string;             /* work [output pointer] */
2048
  const Unit *up = dn->lsu + D2U (dn->digits) - 1;      /* -> msu [input pointer] */
2049
  uInt u, pow;                  /* work */
2050
 
2051
#if DECCHECK
2052
  if (decCheckOperands (DECUNUSED, dn, DECUNUSED, DECUNUSED))
2053
    {
2054
      strcpy (string, "?");
2055
      return;
2056
    }
2057
#endif
2058
 
2059
  if (decNumberIsNegative (dn))
2060
    {                           /* Negatives get a minus (except */
2061
      *c = '-';                 /* NaNs, which remove the '-' below) */
2062
      c++;
2063
    }
2064
  if (dn->bits & DECSPECIAL)
2065
    {                           /* Is a special value */
2066
      if (decNumberIsInfinite (dn))
2067
        {
2068
          strcpy (c, "Infinity");
2069
          return;
2070
        }
2071
      /* a NaN */
2072
      if (dn->bits & DECSNAN)
2073
        {                       /* signalling NaN */
2074
          *c = 's';
2075
          c++;
2076
        }
2077
      strcpy (c, "NaN");
2078
      c += 3;                   /* step past */
2079
      /* if not a clean non-zero coefficient, that's all we have in a */
2080
      /* NaN string */
2081
      if (exp != 0 || (*dn->lsu == 0 && dn->digits == 1))
2082
        return;
2083
      /* [drop through to add integer] */
2084
    }
2085
 
2086
  /* calculate how many digits in msu, and hence first cut */
2087
  cut = dn->digits % DECDPUN;
2088
  if (cut == 0)
2089
    cut = DECDPUN;              /* msu is full */
2090
  cut--;                        /* power of ten for digit */
2091
 
2092
  if (exp == 0)
2093
    {                           /* simple integer [common fastpath, */
2094
      /*   used for NaNs, too] */
2095
      for (; up >= dn->lsu; up--)
2096
        {                       /* each Unit from msu */
2097
          u = *up;              /* contains DECDPUN digits to lay out */
2098
          for (; cut >= 0; c++, cut--)
2099
            TODIGIT (u, cut, c);
2100
          cut = DECDPUN - 1;    /* next Unit has all digits */
2101
        }
2102
      *c = '\0';                /* terminate the string */
2103
      return;
2104
    }
2105
 
2106
  /* non-0 exponent -- assume plain form */
2107
  pre = dn->digits + exp;       /* digits before '.' */
2108
  e = 0;                 /* no E */
2109
  if ((exp > 0) || (pre < -5))
2110
    {                           /* need exponential form */
2111
      e = exp + dn->digits - 1; /* calculate E value */
2112
      pre = 1;                  /* assume one digit before '.' */
2113
      if (eng && (e != 0))
2114
        {                       /* may need to adjust */
2115
          Int adj;              /* adjustment */
2116
          /* The C remainder operator is undefined for negative numbers, so */
2117
          /* we must use positive remainder calculation here */
2118
          if (e < 0)
2119
            {
2120
              adj = (-e) % 3;
2121
              if (adj != 0)
2122
                adj = 3 - adj;
2123
            }
2124
          else
2125
            {                   /* e>0 */
2126
              adj = e % 3;
2127
            }
2128
          e = e - adj;
2129
          /* if we are dealing with zero we will use exponent which is a */
2130
          /* multiple of three, as expected, but there will only be the */
2131
          /* one zero before the E, still.  Otherwise note the padding. */
2132
          if (!ISZERO (dn))
2133
            pre += adj;
2134
          else
2135
            {                   /* is zero */
2136
              if (adj != 0)
2137
                {               /* 0.00Esnn needed */
2138
                  e = e + 3;
2139
                  pre = -(2 - adj);
2140
                }
2141
            }                   /* zero */
2142
        }                       /* eng */
2143
    }
2144
 
2145
  /* lay out the digits of the coefficient, adding 0s and . as needed */
2146
  u = *up;
2147
  if (pre > 0)
2148
    {                           /* xxx.xxx or xx00 (engineering) form */
2149
      for (; pre > 0; pre--, c++, cut--)
2150
        {
2151
          if (cut < 0)
2152
            {                   /* need new Unit */
2153
              if (up == dn->lsu)
2154
                break;          /* out of input digits (pre>digits) */
2155
              up--;
2156
              cut = DECDPUN - 1;
2157
              u = *up;
2158
            }
2159
          TODIGIT (u, cut, c);
2160
        }
2161
      if (up > dn->lsu || (up == dn->lsu && cut >= 0))
2162
        {                       /* more to come, after '.' */
2163
          *c = '.';
2164
          c++;
2165
          for (;; c++, cut--)
2166
            {
2167
              if (cut < 0)
2168
                {               /* need new Unit */
2169
                  if (up == dn->lsu)
2170
                    break;      /* out of input digits */
2171
                  up--;
2172
                  cut = DECDPUN - 1;
2173
                  u = *up;
2174
                }
2175
              TODIGIT (u, cut, c);
2176
            }
2177
        }
2178
      else
2179
        for (; pre > 0; pre--, c++)
2180
          *c = '0';             /* 0 padding (for engineering) needed */
2181
    }
2182
  else
2183
    {                           /* 0.xxx or 0.000xxx form */
2184
      *c = '0';
2185
      c++;
2186
      *c = '.';
2187
      c++;
2188
      for (; pre < 0; pre++, c++)
2189
        *c = '0';               /* add any 0's after '.' */
2190
      for (;; c++, cut--)
2191
        {
2192
          if (cut < 0)
2193
            {                   /* need new Unit */
2194
              if (up == dn->lsu)
2195
                break;          /* out of input digits */
2196
              up--;
2197
              cut = DECDPUN - 1;
2198
              u = *up;
2199
            }
2200
          TODIGIT (u, cut, c);
2201
        }
2202
    }
2203
 
2204
  /* Finally add the E-part, if needed.  It will never be 0, has a
2205
     base maximum and minimum of +999999999 through -999999999, but
2206
     could range down to -1999999998 for subnormal numbers */
2207
  if (e != 0)
2208
    {
2209
      Flag had = 0;              /* 1=had non-zero */
2210
      *c = 'E';
2211
      c++;
2212
      *c = '+';
2213
      c++;                      /* assume positive */
2214
      u = e;                    /* .. */
2215
      if (e < 0)
2216
        {
2217
          *(c - 1) = '-';       /* oops, need - */
2218
          u = -e;               /* uInt, please */
2219
        }
2220
      /* layout the exponent (_itoa is not ANSI C) */
2221
      for (cut = 9; cut >= 0; cut--)
2222
        {
2223
          TODIGIT (u, cut, c);
2224
          if (*c == '0' && !had)
2225
            continue;           /* skip leading zeros */
2226
          had = 1;              /* had non-0 */
2227
          c++;                  /* step for next */
2228
        }                       /* cut */
2229
    }
2230
  *c = '\0';                    /* terminate the string (all paths) */
2231
  return;
2232
}
2233
 
2234
/* ------------------------------------------------------------------ */
2235
/* decAddOp -- add/subtract operation                                 */
2236
/*                                                                    */
2237
/*   This computes C = A + B                                          */
2238
/*                                                                    */
2239
/*   res is C, the result.  C may be A and/or B (e.g., X=X+X)         */
2240
/*   lhs is A                                                         */
2241
/*   rhs is B                                                         */
2242
/*   set is the context                                               */
2243
/*   negate is DECNEG if rhs should be negated, or 0 otherwise        */
2244
/*   status accumulates status for the caller                         */
2245
/*                                                                    */
2246
/* C must have space for set->digits digits.                          */
2247
/* ------------------------------------------------------------------ */
2248
/* If possible, we calculate the coefficient directly into C.         */
2249
/* However, if:                                                       */
2250
/*   -- we need a digits+1 calculation because numbers are unaligned  */
2251
/*      and span more than set->digits digits                         */
2252
/*   -- a carry to digits+1 digits looks possible                     */
2253
/*   -- C is the same as A or B, and the result would destructively   */
2254
/*      overlap the A or B coefficient                                */
2255
/* then we must calculate into a temporary buffer.  In this latter    */
2256
/* case we use the local (stack) buffer if possible, and only if too  */
2257
/* long for that do we resort to malloc.                              */
2258
/*                                                                    */
2259
/* Misalignment is handled as follows:                                */
2260
/*   Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp.    */
2261
/*   BPad: Apply the padding by a combination of shifting (whole      */
2262
/*         units) and multiplication (part units).                    */
2263
/*                                                                    */
2264
/* Addition, especially x=x+1, is speed-critical, so we take pains    */
2265
/* to make returning as fast as possible, by flagging any allocation. */
2266
/* ------------------------------------------------------------------ */
2267
static decNumber *
2268
decAddOp (decNumber * res, const decNumber * lhs,
2269
          const decNumber * rhs, decContext * set, uByte negate, uInt * status)
2270
{
2271
  decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
2272
  decNumber *allocrhs = NULL;   /* .., rhs */
2273
  Int rhsshift;                 /* working shift (in Units) */
2274
  Int maxdigits;                /* longest logical length */
2275
  Int mult;                     /* multiplier */
2276
  Int residue;                  /* rounding accumulator */
2277
  uByte bits;                   /* result bits */
2278
  Flag diffsign;                /* non-0 if arguments have different sign */
2279
  Unit *acc;                    /* accumulator for result */
2280
  Unit accbuff[D2U (DECBUFFER + 1)];    /* local buffer [+1 is for possible */
2281
  /* final carry digit or DECBUFFER=0] */
2282
  Unit *allocacc = NULL;        /* -> allocated acc buffer, iff allocated */
2283
  Flag alloced = 0;              /* set non-0 if any allocations */
2284
  Int reqdigits = set->digits;  /* local copy; requested DIGITS */
2285
  uByte merged;                 /* merged flags */
2286
  Int padding;                  /* work */
2287
 
2288
#if DECCHECK
2289
  if (decCheckOperands (res, lhs, rhs, set))
2290
    return res;
2291
#endif
2292
 
2293
  do
2294
    {                           /* protect allocated storage */
2295
#if DECSUBSET
2296
      if (!set->extended)
2297
        {
2298
          /* reduce operands and set lostDigits status, as needed */
2299
          if (lhs->digits > reqdigits)
2300
            {
2301
              alloclhs = decRoundOperand (lhs, set, status);
2302
              if (alloclhs == NULL)
2303
                break;
2304
              lhs = alloclhs;
2305
              alloced = 1;
2306
            }
2307
          if (rhs->digits > reqdigits)
2308
            {
2309
              allocrhs = decRoundOperand (rhs, set, status);
2310
              if (allocrhs == NULL)
2311
                break;
2312
              rhs = allocrhs;
2313
              alloced = 1;
2314
            }
2315
        }
2316
#endif
2317
      /* [following code does not require input rounding] */
2318
 
2319
      /* note whether signs differ */
2320
      diffsign = (Flag) ((lhs->bits ^ rhs->bits ^ negate) & DECNEG);
2321
 
2322
      /* handle infinities and NaNs */
2323
      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2324
      if (merged)
2325
        {                       /* a special bit set */
2326
          if (merged & (DECSNAN | DECNAN))      /* a NaN */
2327
            decNaNs (res, lhs, rhs, status);
2328
          else
2329
            {                   /* one or two infinities */
2330
              if (decNumberIsInfinite (lhs))
2331
                {               /* LHS is infinity */
2332
                  /* two infinities with different signs is invalid */
2333
                  if (decNumberIsInfinite (rhs) && diffsign)
2334
                    {
2335
                      *status |= DEC_Invalid_operation;
2336
                      break;
2337
                    }
2338
                  bits = lhs->bits & DECNEG;    /* get sign from LHS */
2339
                }
2340
              else
2341
                bits = (rhs->bits ^ negate) & DECNEG;   /* RHS must be Infinity */
2342
              bits |= DECINF;
2343
              decNumberZero (res);
2344
              res->bits = bits; /* set +/- infinity */
2345
            }                   /* an infinity */
2346
          break;
2347
        }
2348
 
2349
      /* Quick exit for add 0s; return the non-0, modified as need be */
2350
      if (ISZERO (lhs))
2351
        {
2352
          Int adjust;           /* work */
2353
          Int lexp = lhs->exponent;     /* save in case LHS==RES */
2354
          bits = lhs->bits;     /* .. */
2355
          residue = 0;           /* clear accumulator */
2356
          decCopyFit (res, rhs, set, &residue, status); /* copy (as needed) */
2357
          res->bits ^= negate;  /* flip if rhs was negated */
2358
#if DECSUBSET
2359
          if (set->extended)
2360
            {                   /* exponents on zeros count */
2361
#endif
2362
              /* exponent will be the lower of the two */
2363
              adjust = lexp - res->exponent;    /* adjustment needed [if -ve] */
2364
              if (ISZERO (res))
2365
                {               /* both 0: special IEEE 854 rules */
2366
                  if (adjust < 0)
2367
                    res->exponent = lexp;       /* set exponent */
2368
                  /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */
2369
                  if (diffsign)
2370
                    {
2371
                      if (set->round != DEC_ROUND_FLOOR)
2372
                        res->bits = 0;
2373
                      else
2374
                        res->bits = DECNEG;     /* preserve 0 sign */
2375
                    }
2376
                }
2377
              else
2378
                {               /* non-0 res */
2379
                  if (adjust < 0)
2380
                    {           /* 0-padding needed */
2381
                      if ((res->digits - adjust) > set->digits)
2382
                        {
2383
                          adjust = res->digits - set->digits;   /* to fit exactly */
2384
                          *status |= DEC_Rounded;       /* [but exact] */
2385
                        }
2386
                      res->digits =
2387
                        decShiftToMost (res->lsu, res->digits, -adjust);
2388
                      res->exponent += adjust;  /* set the exponent. */
2389
                    }
2390
                }               /* non-0 res */
2391
#if DECSUBSET
2392
            }                   /* extended */
2393
#endif
2394
          decFinish (res, set, &residue, status);       /* clean and finalize */
2395
          break;
2396
        }
2397
 
2398
      if (ISZERO (rhs))
2399
        {                       /* [lhs is non-zero] */
2400
          Int adjust;           /* work */
2401
          Int rexp = rhs->exponent;     /* save in case RHS==RES */
2402
          bits = rhs->bits;     /* be clean */
2403
          residue = 0;           /* clear accumulator */
2404
          decCopyFit (res, lhs, set, &residue, status); /* copy (as needed) */
2405
#if DECSUBSET
2406
          if (set->extended)
2407
            {                   /* exponents on zeros count */
2408
#endif
2409
              /* exponent will be the lower of the two */
2410
              /* [0-0 case handled above] */
2411
              adjust = rexp - res->exponent;    /* adjustment needed [if -ve] */
2412
              if (adjust < 0)
2413
                {               /* 0-padding needed */
2414
                  if ((res->digits - adjust) > set->digits)
2415
                    {
2416
                      adjust = res->digits - set->digits;       /* to fit exactly */
2417
                      *status |= DEC_Rounded;   /* [but exact] */
2418
                    }
2419
                  res->digits =
2420
                    decShiftToMost (res->lsu, res->digits, -adjust);
2421
                  res->exponent += adjust;      /* set the exponent. */
2422
                }
2423
#if DECSUBSET
2424
            }                   /* extended */
2425
#endif
2426
          decFinish (res, set, &residue, status);       /* clean and finalize */
2427
          break;
2428
        }
2429
      /* [both fastpath and mainpath code below assume these cases */
2430
      /* (notably 0-0) have already been handled] */
2431
 
2432
      /* calculate the padding needed to align the operands */
2433
      padding = rhs->exponent - lhs->exponent;
2434
 
2435
      /* Fastpath cases where the numbers are aligned and normal, the RHS */
2436
      /* is all in one unit, no operand rounding is needed, and no carry, */
2437
      /* lengthening, or borrow is needed */
2438
      if (rhs->digits <= DECDPUN && padding == 0 && rhs->exponent >= set->emin   /* [some normals drop through] */
2439
          && rhs->digits <= reqdigits && lhs->digits <= reqdigits)
2440
        {
2441
          Int partial = *lhs->lsu;
2442
          if (!diffsign)
2443
            {                   /* adding */
2444
              Int maxv = DECDPUNMAX;    /* highest no-overflow */
2445
              if (lhs->digits < DECDPUN)
2446
                maxv = powers[lhs->digits] - 1;
2447
              partial += *rhs->lsu;
2448
              if (partial <= maxv)
2449
                {               /* no carry */
2450
                  if (res != lhs)
2451
                    decNumberCopy (res, lhs);   /* not in place */
2452
                  *res->lsu = (Unit) partial;   /* [copy could have overwritten RHS] */
2453
                  break;
2454
                }
2455
              /* else drop out for careful add */
2456
            }
2457
          else
2458
            {                   /* signs differ */
2459
              partial -= *rhs->lsu;
2460
              if (partial > 0)
2461
                {               /* no borrow needed, and non-0 result */
2462
                  if (res != lhs)
2463
                    decNumberCopy (res, lhs);   /* not in place */
2464
                  *res->lsu = (Unit) partial;
2465
                  /* this could have reduced digits [but result>0] */
2466
                  res->digits = decGetDigits (res->lsu, D2U (res->digits));
2467
                  break;
2468
                }
2469
              /* else drop out for careful subtract */
2470
            }
2471
        }
2472
 
2473
      /* Now align (pad) the lhs or rhs so we can add or subtract them, as
2474
         necessary.  If one number is much larger than the other (that is,
2475
         if in plain form there is a least one digit between the lowest
2476
         digit or one and the highest of the other) we need to pad with up
2477
         to DIGITS-1 trailing zeros, and then apply rounding (as exotic
2478
         rounding modes may be affected by the residue).
2479
       */
2480
      rhsshift = 0;              /* rhs shift to left (padding) in Units */
2481
      bits = lhs->bits;         /* assume sign is that of LHS */
2482
      mult = 1;                 /* likely multiplier */
2483
 
2484
      /* if padding==0 the operands are aligned; no padding needed */
2485
      if (padding != 0)
2486
        {
2487
          /* some padding needed */
2488
          /* We always pad the RHS, as we can then effect any required */
2489
          /* padding by a combination of shifts and a multiply */
2490
          Flag swapped = 0;
2491
          if (padding < 0)
2492
            {                   /* LHS needs the padding */
2493
              const decNumber *t;
2494
              padding = -padding;       /* will be +ve */
2495
              bits = (uByte) (rhs->bits ^ negate);      /* assumed sign is now that of RHS */
2496
              t = lhs;
2497
              lhs = rhs;
2498
              rhs = t;
2499
              swapped = 1;
2500
            }
2501
 
2502
          /* If, after pad, rhs would be longer than lhs by digits+1 or */
2503
          /* more then lhs cannot affect the answer, except as a residue, */
2504
          /* so we only need to pad up to a length of DIGITS+1. */
2505
          if (rhs->digits + padding > lhs->digits + reqdigits + 1)
2506
            {
2507
              /* The RHS is sufficient */
2508
              /* for residue we use the relative sign indication... */
2509
              Int shift = reqdigits - rhs->digits;      /* left shift needed */
2510
              residue = 1;      /* residue for rounding */
2511
              if (diffsign)
2512
                residue = -residue;     /* signs differ */
2513
              /* copy, shortening if necessary */
2514
              decCopyFit (res, rhs, set, &residue, status);
2515
              /* if it was already shorter, then need to pad with zeros */
2516
              if (shift > 0)
2517
                {
2518
                  res->digits = decShiftToMost (res->lsu, res->digits, shift);
2519
                  res->exponent -= shift;       /* adjust the exponent. */
2520
                }
2521
              /* flip the result sign if unswapped and rhs was negated */
2522
              if (!swapped)
2523
                res->bits ^= negate;
2524
              decFinish (res, set, &residue, status);   /* done */
2525
              break;
2526
            }
2527
 
2528
          /* LHS digits may affect result */
2529
          rhsshift = D2U (padding + 1) - 1;     /* this much by Unit shift .. */
2530
          mult = powers[padding - (rhsshift * DECDPUN)];        /* .. this by multiplication */
2531
        }                       /* padding needed */
2532
 
2533
      if (diffsign)
2534
        mult = -mult;           /* signs differ */
2535
 
2536
      /* determine the longer operand */
2537
      maxdigits = rhs->digits + padding;        /* virtual length of RHS */
2538
      if (lhs->digits > maxdigits)
2539
        maxdigits = lhs->digits;
2540
 
2541
      /* Decide on the result buffer to use; if possible place directly */
2542
      /* into result. */
2543
      acc = res->lsu;           /* assume build direct */
2544
      /* If destructive overlap, or the number is too long, or a carry or */
2545
      /* borrow to DIGITS+1 might be possible we must use a buffer. */
2546
      /* [Might be worth more sophisticated tests when maxdigits==reqdigits] */
2547
      if ((maxdigits >= reqdigits)      /* is, or could be, too large */
2548
          || (res == rhs && rhsshift > 0))
2549
        {                       /* destructive overlap */
2550
          /* buffer needed; choose it */
2551
          /* we'll need units for maxdigits digits, +1 Unit for carry or borrow */
2552
          Int need = D2U (maxdigits) + 1;
2553
          acc = accbuff;        /* assume use local buffer */
2554
          if (need * sizeof (Unit) > sizeof (accbuff))
2555
            {
2556
              allocacc = (Unit *) malloc (need * sizeof (Unit));
2557
              if (allocacc == NULL)
2558
                {               /* hopeless -- abandon */
2559
                  *status |= DEC_Insufficient_storage;
2560
                  break;
2561
                }
2562
              acc = allocacc;
2563
              alloced = 1;
2564
            }
2565
        }
2566
 
2567
      res->bits = (uByte) (bits & DECNEG);      /* it's now safe to overwrite.. */
2568
      res->exponent = lhs->exponent;    /* .. operands (even if aliased) */
2569
 
2570
#if DECTRACE
2571
      decDumpAr ('A', lhs->lsu, D2U (lhs->digits));
2572
      decDumpAr ('B', rhs->lsu, D2U (rhs->digits));
2573
      printf ("  :h: %d %d\n", rhsshift, mult);
2574
#endif
2575
 
2576
      /* add [A+B*m] or subtract [A+B*(-m)] */
2577
      res->digits = decUnitAddSub (lhs->lsu, D2U (lhs->digits), rhs->lsu, D2U (rhs->digits), rhsshift, acc, mult) * DECDPUN;    /* [units -> digits] */
2578
      if (res->digits < 0)
2579
        {                       /* we borrowed */
2580
          res->digits = -res->digits;
2581
          res->bits ^= DECNEG;  /* flip the sign */
2582
        }
2583
#if DECTRACE
2584
      decDumpAr ('+', acc, D2U (res->digits));
2585
#endif
2586
 
2587
      /* If we used a buffer we need to copy back, possibly shortening */
2588
      /* (If we didn't use buffer it must have fit, so can't need rounding */
2589
      /* and residue must be 0.) */
2590
      residue = 0;               /* clear accumulator */
2591
      if (acc != res->lsu)
2592
        {
2593
#if DECSUBSET
2594
          if (set->extended)
2595
            {                   /* round from first significant digit */
2596
#endif
2597
              /* remove leading zeros that we added due to rounding up to */
2598
              /* integral Units -- before the test for rounding. */
2599
              if (res->digits > reqdigits)
2600
                res->digits = decGetDigits (acc, D2U (res->digits));
2601
              decSetCoeff (res, set, acc, res->digits, &residue, status);
2602
#if DECSUBSET
2603
            }
2604
          else
2605
            {                   /* subset arithmetic rounds from original significant digit */
2606
              /* We may have an underestimate.  This only occurs when both */
2607
              /* numbers fit in DECDPUN digits and we are padding with a */
2608
              /* negative multiple (-10, -100...) and the top digit(s) become */
2609
              /* 0.  (This only matters if we are using X3.274 rules where the */
2610
              /* leading zero could be included in the rounding.) */
2611
              if (res->digits < maxdigits)
2612
                {
2613
                  *(acc + D2U (res->digits)) = 0;        /* ensure leading 0 is there */
2614
                  res->digits = maxdigits;
2615
                }
2616
              else
2617
                {
2618
                  /* remove leading zeros that we added due to rounding up to */
2619
                  /* integral Units (but only those in excess of the original */
2620
                  /* maxdigits length, unless extended) before test for rounding. */
2621
                  if (res->digits > reqdigits)
2622
                    {
2623
                      res->digits = decGetDigits (acc, D2U (res->digits));
2624
                      if (res->digits < maxdigits)
2625
                        res->digits = maxdigits;
2626
                    }
2627
                }
2628
              decSetCoeff (res, set, acc, res->digits, &residue, status);
2629
              /* Now apply rounding if needed before removing leading zeros. */
2630
              /* This is safe because subnormals are not a possibility */
2631
              if (residue != 0)
2632
                {
2633
                  decApplyRound (res, set, residue, status);
2634
                  residue = 0;   /* we did what we had to do */
2635
                }
2636
            }                   /* subset */
2637
#endif
2638
        }                       /* used buffer */
2639
 
2640
      /* strip leading zeros [these were left on in case of subset subtract] */
2641
      res->digits = decGetDigits (res->lsu, D2U (res->digits));
2642
 
2643
      /* apply checks and rounding */
2644
      decFinish (res, set, &residue, status);
2645
 
2646
      /* "When the sum of two operands with opposite signs is exactly */
2647
      /* zero, the sign of that sum shall be '+' in all rounding modes */
2648
      /* except round toward -Infinity, in which mode that sign shall be */
2649
      /* '-'."  [Subset zeros also never have '-', set by decFinish.] */
2650
      if (ISZERO (res) && diffsign
2651
#if DECSUBSET
2652
          && set->extended
2653
#endif
2654
          && (*status & DEC_Inexact) == 0)
2655
        {
2656
          if (set->round == DEC_ROUND_FLOOR)
2657
            res->bits |= DECNEG;        /* sign - */
2658
          else
2659
            res->bits &= ~DECNEG;       /* sign + */
2660
        }
2661
    }
2662
  while (0);                     /* end protected */
2663
 
2664
  if (alloced)
2665
    {
2666
      if (allocacc != NULL)
2667
        free (allocacc);        /* drop any storage we used */
2668
      if (allocrhs != NULL)
2669
        free (allocrhs);        /* .. */
2670
      if (alloclhs != NULL)
2671
        free (alloclhs);        /* .. */
2672
    }
2673
  return res;
2674
}
2675
 
2676
/* ------------------------------------------------------------------ */
2677
/* decDivideOp -- division operation                                  */
2678
/*                                                                    */
2679
/*  This routine performs the calculations for all four division      */
2680
/*  operators (divide, divideInteger, remainder, remainderNear).      */
2681
/*                                                                    */
2682
/*  C=A op B                                                          */
2683
/*                                                                    */
2684
/*   res is C, the result.  C may be A and/or B (e.g., X=X/X)         */
2685
/*   lhs is A                                                         */
2686
/*   rhs is B                                                         */
2687
/*   set is the context                                               */
2688
/*   op  is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively.    */
2689
/*   status is the usual accumulator                                  */
2690
/*                                                                    */
2691
/* C must have space for set->digits digits.                          */
2692
/*                                                                    */
2693
/* ------------------------------------------------------------------ */
2694
/*   The underlying algorithm of this routine is the same as in the   */
2695
/*   1981 S/370 implementation, that is, non-restoring long division  */
2696
/*   with bi-unit (rather than bi-digit) estimation for each unit     */
2697
/*   multiplier.  In this pseudocode overview, complications for the  */
2698
/*   Remainder operators and division residues for exact rounding are */
2699
/*   omitted for clarity.                                             */
2700
/*                                                                    */
2701
/*     Prepare operands and handle special values                     */
2702
/*     Test for x/0 and then 0/x                                      */
2703
/*     Exp =Exp1 - Exp2                                               */
2704
/*     Exp =Exp +len(var1) -len(var2)                                 */
2705
/*     Sign=Sign1 * Sign2                                             */
2706
/*     Pad accumulator (Var1) to double-length with 0's (pad1)        */
2707
/*     Pad Var2 to same length as Var1                                */
2708
/*     msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round  */
2709
/*     have=0                                                         */
2710
/*     Do until (have=digits+1 OR residue=0)                          */
2711
/*       if exp<0 then if integer divide/residue then leave           */
2712
/*       this_unit=0                                                  */
2713
/*       Do forever                                                   */
2714
/*          compare numbers                                           */
2715
/*          if <0 then leave inner_loop                               */
2716
/*          if =0 then (* quick exit without subtract *) do           */
2717
/*             this_unit=this_unit+1; output this_unit                */
2718
/*             leave outer_loop; end                                  */
2719
/*          Compare lengths of numbers (mantissae):                   */
2720
/*          If same then tops2=msu2pair -- {units 1&2 of var2}        */
2721
/*                  else tops2=msu2plus -- {0, unit 1 of var2}        */
2722
/*          tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
2723
/*          mult=tops1/tops2  -- Good and safe guess at divisor       */
2724
/*          if mult=0 then mult=1                                     */
2725
/*          this_unit=this_unit+mult                                  */
2726
/*          subtract                                                  */
2727
/*          end inner_loop                                            */
2728
/*        if have\=0 | this_unit\=0 then do                           */
2729
/*          output this_unit                                          */
2730
/*          have=have+1; end                                          */
2731
/*        var2=var2/10                                                */
2732
/*        exp=exp-1                                                   */
2733
/*        end outer_loop                                              */
2734
/*     exp=exp+1   -- set the proper exponent                         */
2735
/*     if have=0 then generate answer=0                               */
2736
/*     Return (Result is defined by Var1)                             */
2737
/*                                                                    */
2738
/* ------------------------------------------------------------------ */
2739
/* We need two working buffers during the long division; one (digits+ */
2740
/* 1) to accumulate the result, and the other (up to 2*digits+1) for  */
2741
/* long subtractions.  These are acc and var1 respectively.           */
2742
/* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
2743
/* ------------------------------------------------------------------ */
2744
static decNumber *
2745
decDivideOp (decNumber * res,
2746
             const decNumber * lhs, const decNumber * rhs,
2747
             decContext * set, Flag op, uInt * status)
2748
{
2749
  decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
2750
  decNumber *allocrhs = NULL;   /* .., rhs */
2751
  Unit accbuff[D2U (DECBUFFER + DECDPUN)];      /* local buffer */
2752
  Unit *acc = accbuff;          /* -> accumulator array for result */
2753
  Unit *allocacc = NULL;        /* -> allocated buffer, iff allocated */
2754
  Unit *accnext;                /* -> where next digit will go */
2755
  Int acclength;                /* length of acc needed [Units] */
2756
  Int accunits;                 /* count of units accumulated */
2757
  Int accdigits;                /* count of digits accumulated */
2758
 
2759
  Unit varbuff[D2U (DECBUFFER * 2 + DECDPUN) * sizeof (Unit)];  /* buffer for var1 */
2760
  Unit *var1 = varbuff;         /* -> var1 array for long subtraction */
2761
  Unit *varalloc = NULL;        /* -> allocated buffer, iff used */
2762
 
2763
  const Unit *var2;             /* -> var2 array */
2764
 
2765
  Int var1units, var2units;     /* actual lengths */
2766
  Int var2ulen;                 /* logical length (units) */
2767
  Int var1initpad = 0;           /* var1 initial padding (digits) */
2768
  Unit *msu1;                   /* -> msu of each var */
2769
  const Unit *msu2;             /* -> msu of each var */
2770
  Int msu2plus;                 /* msu2 plus one [does not vary] */
2771
  eInt msu2pair;                /* msu2 pair plus one [does not vary] */
2772
  Int maxdigits;                /* longest LHS or required acc length */
2773
  Int mult;                     /* multiplier for subtraction */
2774
  Unit thisunit;                /* current unit being accumulated */
2775
  Int residue;                  /* for rounding */
2776
  Int reqdigits = set->digits;  /* requested DIGITS */
2777
  Int exponent;                 /* working exponent */
2778
  Int maxexponent = 0;           /* DIVIDE maximum exponent if unrounded */
2779
  uByte bits;                   /* working sign */
2780
  uByte merged;                 /* merged flags */
2781
  Unit *target;                 /* work */
2782
  const Unit *source;           /* work */
2783
  uInt const *pow;              /* .. */
2784
  Int shift, cut;               /* .. */
2785
#if DECSUBSET
2786
  Int dropped;                  /* work */
2787
#endif
2788
 
2789
#if DECCHECK
2790
  if (decCheckOperands (res, lhs, rhs, set))
2791
    return res;
2792
#endif
2793
 
2794
  do
2795
    {                           /* protect allocated storage */
2796
#if DECSUBSET
2797
      if (!set->extended)
2798
        {
2799
          /* reduce operands and set lostDigits status, as needed */
2800
          if (lhs->digits > reqdigits)
2801
            {
2802
              alloclhs = decRoundOperand (lhs, set, status);
2803
              if (alloclhs == NULL)
2804
                break;
2805
              lhs = alloclhs;
2806
            }
2807
          if (rhs->digits > reqdigits)
2808
            {
2809
              allocrhs = decRoundOperand (rhs, set, status);
2810
              if (allocrhs == NULL)
2811
                break;
2812
              rhs = allocrhs;
2813
            }
2814
        }
2815
#endif
2816
      /* [following code does not require input rounding] */
2817
 
2818
      bits = (lhs->bits ^ rhs->bits) & DECNEG;  /* assumed sign for divisions */
2819
 
2820
      /* handle infinities and NaNs */
2821
      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
2822
      if (merged)
2823
        {                       /* a special bit set */
2824
          if (merged & (DECSNAN | DECNAN))
2825
            {                   /* one or two NaNs */
2826
              decNaNs (res, lhs, rhs, status);
2827
              break;
2828
            }
2829
          /* one or two infinities */
2830
          if (decNumberIsInfinite (lhs))
2831
            {                   /* LHS (dividend) is infinite */
2832
              if (decNumberIsInfinite (rhs) ||  /* two infinities are invalid .. */
2833
                  op & (REMAINDER | REMNEAR))
2834
                {               /* as is remainder of infinity */
2835
                  *status |= DEC_Invalid_operation;
2836
                  break;
2837
                }
2838
              /* [Note that infinity/0 raises no exceptions] */
2839
              decNumberZero (res);
2840
              res->bits = bits | DECINF;        /* set +/- infinity */
2841
              break;
2842
            }
2843
          else
2844
            {                   /* RHS (divisor) is infinite */
2845
              residue = 0;
2846
              if (op & (REMAINDER | REMNEAR))
2847
                {
2848
                  /* result is [finished clone of] lhs */
2849
                  decCopyFit (res, lhs, set, &residue, status);
2850
                }
2851
              else
2852
                {               /* a division */
2853
                  decNumberZero (res);
2854
                  res->bits = bits;     /* set +/- zero */
2855
                  /* for DIVIDEINT the exponent is always 0.  For DIVIDE, result */
2856
                  /* is a 0 with infinitely negative exponent, clamped to minimum */
2857
                  if (op & DIVIDE)
2858
                    {
2859
                      res->exponent = set->emin - set->digits + 1;
2860
                      *status |= DEC_Clamped;
2861
                    }
2862
                }
2863
              decFinish (res, set, &residue, status);
2864
              break;
2865
            }
2866
        }
2867
 
2868
      /* handle 0 rhs (x/0) */
2869
      if (ISZERO (rhs))
2870
        {                       /* x/0 is always exceptional */
2871
          if (ISZERO (lhs))
2872
            {
2873
              decNumberZero (res);      /* [after lhs test] */
2874
              *status |= DEC_Division_undefined;        /* 0/0 will become NaN */
2875
            }
2876
          else
2877
            {
2878
              decNumberZero (res);
2879
              if (op & (REMAINDER | REMNEAR))
2880
                *status |= DEC_Invalid_operation;
2881
              else
2882
                {
2883
                  *status |= DEC_Division_by_zero;      /* x/0 */
2884
                  res->bits = bits | DECINF;    /* .. is +/- Infinity */
2885
                }
2886
            }
2887
          break;
2888
        }
2889
 
2890
      /* handle 0 lhs (0/x) */
2891
      if (ISZERO (lhs))
2892
        {                       /* 0/x [x!=0] */
2893
#if DECSUBSET
2894
          if (!set->extended)
2895
            decNumberZero (res);
2896
          else
2897
            {
2898
#endif
2899
              if (op & DIVIDE)
2900
                {
2901
                  residue = 0;
2902
                  exponent = lhs->exponent - rhs->exponent;     /* ideal exponent */
2903
                  decNumberCopy (res, lhs);     /* [zeros always fit] */
2904
                  res->bits = bits;     /* sign as computed */
2905
                  res->exponent = exponent;     /* exponent, too */
2906
                  decFinalize (res, set, &residue, status);     /* check exponent */
2907
                }
2908
              else if (op & DIVIDEINT)
2909
                {
2910
                  decNumberZero (res);  /* integer 0 */
2911
                  res->bits = bits;     /* sign as computed */
2912
                }
2913
              else
2914
                {               /* a remainder */
2915
                  exponent = rhs->exponent;     /* [save in case overwrite] */
2916
                  decNumberCopy (res, lhs);     /* [zeros always fit] */
2917
                  if (exponent < res->exponent)
2918
                    res->exponent = exponent;   /* use lower */
2919
                }
2920
#if DECSUBSET
2921
            }
2922
#endif
2923
          break;
2924
        }
2925
 
2926
      /* Precalculate exponent.  This starts off adjusted (and hence fits */
2927
      /* in 31 bits) and becomes the usual unadjusted exponent as the */
2928
      /* division proceeds.  The order of evaluation is important, here, */
2929
      /* to avoid wrap. */
2930
      exponent =
2931
        (lhs->exponent + lhs->digits) - (rhs->exponent + rhs->digits);
2932
 
2933
      /* If the working exponent is -ve, then some quick exits are */
2934
      /* possible because the quotient is known to be <1 */
2935
      /* [for REMNEAR, it needs to be < -1, as -0.5 could need work] */
2936
      if (exponent < 0 && !(op == DIVIDE))
2937
        {
2938
          if (op & DIVIDEINT)
2939
            {
2940
              decNumberZero (res);      /* integer part is 0 */
2941
#if DECSUBSET
2942
              if (set->extended)
2943
#endif
2944
                res->bits = bits;       /* set +/- zero */
2945
              break;
2946
            }
2947
          /* we can fastpath remainders so long as the lhs has the */
2948
          /* smaller (or equal) exponent */
2949
          if (lhs->exponent <= rhs->exponent)
2950
            {
2951
              if (op & REMAINDER || exponent < -1)
2952
                {
2953
                  /* It is REMAINDER or safe REMNEAR; result is [finished */
2954
                  /* clone of] lhs  (r = x - 0*y) */
2955
                  residue = 0;
2956
                  decCopyFit (res, lhs, set, &residue, status);
2957
                  decFinish (res, set, &residue, status);
2958
                  break;
2959
                }
2960
              /* [unsafe REMNEAR drops through] */
2961
            }
2962
        }                       /* fastpaths */
2963
 
2964
      /* We need long (slow) division; roll up the sleeves... */
2965
 
2966
      /* The accumulator will hold the quotient of the division. */
2967
      /* If it needs to be too long for stack storage, then allocate. */
2968
      acclength = D2U (reqdigits + DECDPUN);    /* in Units */
2969
      if (acclength * sizeof (Unit) > sizeof (accbuff))
2970
        {
2971
          allocacc = (Unit *) malloc (acclength * sizeof (Unit));
2972
          if (allocacc == NULL)
2973
            {                   /* hopeless -- abandon */
2974
              *status |= DEC_Insufficient_storage;
2975
              break;
2976
            }
2977
          acc = allocacc;       /* use the allocated space */
2978
        }
2979
 
2980
      /* var1 is the padded LHS ready for subtractions. */
2981
      /* If it needs to be too long for stack storage, then allocate. */
2982
      /* The maximum units we need for var1 (long subtraction) is: */
2983
      /* Enough for */
2984
      /*     (rhs->digits+reqdigits-1) -- to allow full slide to right */
2985
      /* or  (lhs->digits)             -- to allow for long lhs */
2986
      /* whichever is larger */
2987
      /*   +1                -- for rounding of slide to right */
2988
      /*   +1                -- for leading 0s */
2989
      /*   +1                -- for pre-adjust if a remainder or DIVIDEINT */
2990
      /* [Note: unused units do not participate in decUnitAddSub data] */
2991
      maxdigits = rhs->digits + reqdigits - 1;
2992
      if (lhs->digits > maxdigits)
2993
        maxdigits = lhs->digits;
2994
      var1units = D2U (maxdigits) + 2;
2995
      /* allocate a guard unit above msu1 for REMAINDERNEAR */
2996
      if (!(op & DIVIDE))
2997
        var1units++;
2998
      if ((var1units + 1) * sizeof (Unit) > sizeof (varbuff))
2999
        {
3000
          varalloc = (Unit *) malloc ((var1units + 1) * sizeof (Unit));
3001
          if (varalloc == NULL)
3002
            {                   /* hopeless -- abandon */
3003
              *status |= DEC_Insufficient_storage;
3004
              break;
3005
            }
3006
          var1 = varalloc;      /* use the allocated space */
3007
        }
3008
 
3009
      /* Extend the lhs and rhs to full long subtraction length.  The lhs */
3010
      /* is truly extended into the var1 buffer, with 0 padding, so we can */
3011
      /* subtract in place.  The rhs (var2) has virtual padding */
3012
      /* (implemented by decUnitAddSub). */
3013
      /* We allocated one guard unit above msu1 for rem=rem+rem in REMAINDERNEAR */
3014
      msu1 = var1 + var1units - 1;      /* msu of var1 */
3015
      source = lhs->lsu + D2U (lhs->digits) - 1;        /* msu of input array */
3016
      for (target = msu1; source >= lhs->lsu; source--, target--)
3017
        *target = *source;
3018
      for (; target >= var1; target--)
3019
        *target = 0;
3020
 
3021
      /* rhs (var2) is left-aligned with var1 at the start */
3022
      var2ulen = var1units;     /* rhs logical length (units) */
3023
      var2units = D2U (rhs->digits);    /* rhs actual length (units) */
3024
      var2 = rhs->lsu;          /* -> rhs array */
3025
      msu2 = var2 + var2units - 1;      /* -> msu of var2 [never changes] */
3026
      /* now set up the variables which we'll use for estimating the */
3027
      /* multiplication factor.  If these variables are not exact, we add */
3028
      /* 1 to make sure that we never overestimate the multiplier. */
3029
      msu2plus = *msu2;         /* it's value .. */
3030
      if (var2units > 1)
3031
        msu2plus++;             /* .. +1 if any more */
3032
      msu2pair = (eInt) * msu2 * (DECDPUNMAX + 1);      /* top two pair .. */
3033
      if (var2units > 1)
3034
        {                       /* .. [else treat 2nd as 0] */
3035
          msu2pair += *(msu2 - 1);      /* .. */
3036
          if (var2units > 2)
3037
            msu2pair++;         /* .. +1 if any more */
3038
        }
3039
 
3040
      /* Since we are working in units, the units may have leading zeros, */
3041
      /* but we calculated the exponent on the assumption that they are */
3042
      /* both left-aligned.  Adjust the exponent to compensate: add the */
3043
      /* number of leading zeros in var1 msu and subtract those in var2 msu. */
3044
      /* [We actually do this by counting the digits and negating, as */
3045
      /* lead1=DECDPUN-digits1, and similarly for lead2.] */
3046
      for (pow = &powers[1]; *msu1 >= *pow; pow++)
3047
        exponent--;
3048
      for (pow = &powers[1]; *msu2 >= *pow; pow++)
3049
        exponent++;
3050
 
3051
      /* Now, if doing an integer divide or remainder, we want to ensure */
3052
      /* that the result will be Unit-aligned.  To do this, we shift the */
3053
      /* var1 accumulator towards least if need be.  (It's much easier to */
3054
      /* do this now than to reassemble the residue afterwards, if we are */
3055
      /* doing a remainder.)  Also ensure the exponent is not negative. */
3056
      if (!(op & DIVIDE))
3057
        {
3058
          Unit *u;
3059
          /* save the initial 'false' padding of var1, in digits */
3060
          var1initpad = (var1units - D2U (lhs->digits)) * DECDPUN;
3061
          /* Determine the shift to do. */
3062
          if (exponent < 0)
3063
            cut = -exponent;
3064
          else
3065
            cut = DECDPUN - exponent % DECDPUN;
3066
          decShiftToLeast (var1, var1units, cut);
3067
          exponent += cut;      /* maintain numerical value */
3068
          var1initpad -= cut;   /* .. and reduce padding */
3069
          /* clean any most-significant units we just emptied */
3070
          for (u = msu1; cut >= DECDPUN; cut -= DECDPUN, u--)
3071
            *u = 0;
3072
        }                       /* align */
3073
      else
3074
        {                       /* is DIVIDE */
3075
          maxexponent = lhs->exponent - rhs->exponent;  /* save */
3076
          /* optimization: if the first iteration will just produce 0, */
3077
          /* preadjust to skip it [valid for DIVIDE only] */
3078
          if (*msu1 < *msu2)
3079
            {
3080
              var2ulen--;       /* shift down */
3081
              exponent -= DECDPUN;      /* update the exponent */
3082
            }
3083
        }
3084
 
3085
      /* ---- start the long-division loops ------------------------------ */
3086
      accunits = 0;              /* no units accumulated yet */
3087
      accdigits = 0;             /* .. or digits */
3088
      accnext = acc + acclength - 1;    /* -> msu of acc [NB: allows digits+1] */
3089
      for (;;)
3090
        {                       /* outer forever loop */
3091
          thisunit = 0;          /* current unit assumed 0 */
3092
          /* find the next unit */
3093
          for (;;)
3094
            {                   /* inner forever loop */
3095
              /* strip leading zero units [from either pre-adjust or from */
3096
              /* subtract last time around].  Leave at least one unit. */
3097
              for (; *msu1 == 0 && msu1 > var1; msu1--)
3098
                var1units--;
3099
 
3100
              if (var1units < var2ulen)
3101
                break;          /* var1 too low for subtract */
3102
              if (var1units == var2ulen)
3103
                {               /* unit-by-unit compare needed */
3104
                  /* compare the two numbers, from msu */
3105
                  Unit *pv1, v2;        /* units to compare */
3106
                  const Unit *pv2;      /* units to compare */
3107
                  pv2 = msu2;   /* -> msu */
3108
                  for (pv1 = msu1;; pv1--, pv2--)
3109
                    {
3110
                      /* v1=*pv1 -- always OK */
3111
                      v2 = 0;    /* assume in padding */
3112
                      if (pv2 >= var2)
3113
                        v2 = *pv2;      /* in range */
3114
                      if (*pv1 != v2)
3115
                        break;  /* no longer the same */
3116
                      if (pv1 == var1)
3117
                        break;  /* done; leave pv1 as is */
3118
                    }
3119
                  /* here when all inspected or a difference seen */
3120
                  if (*pv1 < v2)
3121
                    break;      /* var1 too low to subtract */
3122
                  if (*pv1 == v2)
3123
                    {           /* var1 == var2 */
3124
                      /* reach here if var1 and var2 are identical; subtraction */
3125
                      /* would increase digit by one, and the residue will be 0 so */
3126
                      /* we are done; leave the loop with residue set to 0. */
3127
                      thisunit++;       /* as though subtracted */
3128
                      *var1 = 0; /* set var1 to 0 */
3129
                      var1units = 1;    /* .. */
3130
                      break;    /* from inner */
3131
                    }           /* var1 == var2 */
3132
                  /* *pv1>v2.  Prepare for real subtraction; the lengths are equal */
3133
                  /* Estimate the multiplier (there's always a msu1-1)... */
3134
                  /* Bring in two units of var2 to provide a good estimate. */
3135
                  mult =
3136
                    (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
3137
                            *(msu1 - 1)) / msu2pair);
3138
                }               /* lengths the same */
3139
              else
3140
                {               /* var1units > var2ulen, so subtraction is safe */
3141
                  /* The var2 msu is one unit towards the lsu of the var1 msu, */
3142
                  /* so we can only use one unit for var2. */
3143
                  mult =
3144
                    (Int) (((eInt) * msu1 * (DECDPUNMAX + 1) +
3145
                            *(msu1 - 1)) / msu2plus);
3146
                }
3147
              if (mult == 0)
3148
                mult = 1;       /* must always be at least 1 */
3149
              /* subtraction needed; var1 is > var2 */
3150
              thisunit = (Unit) (thisunit + mult);      /* accumulate */
3151
              /* subtract var1-var2, into var1; only the overlap needs */
3152
              /* processing, as we are in place */
3153
              shift = var2ulen - var2units;
3154
#if DECTRACE
3155
              decDumpAr ('1', &var1[shift], var1units - shift);
3156
              decDumpAr ('2', var2, var2units);
3157
              printf ("m=%d\n", -mult);
3158
#endif
3159
              decUnitAddSub (&var1[shift], var1units - shift,
3160
                             var2, var2units, 0, &var1[shift], -mult);
3161
#if DECTRACE
3162
              decDumpAr ('#', &var1[shift], var1units - shift);
3163
#endif
3164
              /* var1 now probably has leading zeros; these are removed at the */
3165
              /* top of the inner loop. */
3166
            }                   /* inner loop */
3167
 
3168
          /* We have the next unit; unless it's a leading zero, add to acc */
3169
          if (accunits != 0 || thisunit != 0)
3170
            {                   /* put the unit we got */
3171
              *accnext = thisunit;      /* store in accumulator */
3172
              /* account exactly for the digits we got */
3173
              if (accunits == 0)
3174
                {
3175
                  accdigits++;  /* at least one */
3176
                  for (pow = &powers[1]; thisunit >= *pow; pow++)
3177
                    accdigits++;
3178
                }
3179
              else
3180
                accdigits += DECDPUN;
3181
              accunits++;       /* update count */
3182
              accnext--;        /* ready for next */
3183
              if (accdigits > reqdigits)
3184
                break;          /* we have all we need */
3185
            }
3186
 
3187
          /* if the residue is zero, we're done (unless divide or */
3188
          /* divideInteger and we haven't got enough digits yet) */
3189
          if (*var1 == 0 && var1units == 1)
3190
            {                   /* residue is 0 */
3191
              if (op & (REMAINDER | REMNEAR))
3192
                break;
3193
              if ((op & DIVIDE) && (exponent <= maxexponent))
3194
                break;
3195
              /* [drop through if divideInteger] */
3196
            }
3197
          /* we've also done enough if calculating remainder or integer */
3198
          /* divide and we just did the last ('units') unit */
3199
          if (exponent == 0 && !(op & DIVIDE))
3200
            break;
3201
 
3202
          /* to get here, var1 is less than var2, so divide var2 by the per- */
3203
          /* Unit power of ten and go for the next digit */
3204
          var2ulen--;           /* shift down */
3205
          exponent -= DECDPUN;  /* update the exponent */
3206
        }                       /* outer loop */
3207
 
3208
      /* ---- division is complete --------------------------------------- */
3209
      /* here: acc      has at least reqdigits+1 of good results (or fewer */
3210
      /*                if early stop), starting at accnext+1 (its lsu) */
3211
      /*       var1     has any residue at the stopping point */
3212
      /*       accunits is the number of digits we collected in acc */
3213
      if (accunits == 0)
3214
        {                       /* acc is 0 */
3215
          accunits = 1;         /* show we have one .. */
3216
          accdigits = 1;        /* .. */
3217
          *accnext = 0;          /* .. whose value is 0 */
3218
        }
3219
      else
3220
        accnext++;              /* back to last placed */
3221
      /* accnext now -> lowest unit of result */
3222
 
3223
      residue = 0;               /* assume no residue */
3224
      if (op & DIVIDE)
3225
        {
3226
          /* record the presence of any residue, for rounding */
3227
          if (*var1 != 0 || var1units > 1)
3228
            residue = 1;
3229
          else
3230
            {                   /* no residue */
3231
              /* We had an exact division; clean up spurious trailing 0s. */
3232
              /* There will be at most DECDPUN-1, from the final multiply, */
3233
              /* and then only if the result is non-0 (and even) and the */
3234
              /* exponent is 'loose'. */
3235
#if DECDPUN>1
3236
              Unit lsu = *accnext;
3237
              if (!(lsu & 0x01) && (lsu != 0))
3238
                {
3239
                  /* count the trailing zeros */
3240
                  Int drop = 0;
3241
                  for (;; drop++)
3242
                    {           /* [will terminate because lsu!=0] */
3243
                      if (exponent >= maxexponent)
3244
                        break;  /* don't chop real 0s */
3245
#if DECDPUN<=4
3246
                      if ((lsu - QUOT10 (lsu, drop + 1)
3247
                           * powers[drop + 1]) != 0)
3248
                        break;  /* found non-0 digit */
3249
#else
3250
                      if (lsu % powers[drop + 1] != 0)
3251
                        break;  /* found non-0 digit */
3252
#endif
3253
                      exponent++;
3254
                    }
3255
                  if (drop > 0)
3256
                    {
3257
                      accunits = decShiftToLeast (accnext, accunits, drop);
3258
                      accdigits = decGetDigits (accnext, accunits);
3259
                      accunits = D2U (accdigits);
3260
                      /* [exponent was adjusted in the loop] */
3261
                    }
3262
                }               /* neither odd nor 0 */
3263
#endif
3264
            }                   /* exact divide */
3265
        }                       /* divide */
3266
      else                      /* op!=DIVIDE */
3267
        {
3268
          /* check for coefficient overflow */
3269
          if (accdigits + exponent > reqdigits)
3270
            {
3271
              *status |= DEC_Division_impossible;
3272
              break;
3273
            }
3274
          if (op & (REMAINDER | REMNEAR))
3275
            {
3276
              /* [Here, the exponent will be 0, because we adjusted var1 */
3277
              /* appropriately.] */
3278
              Int postshift;    /* work */
3279
              Flag wasodd = 0;   /* integer was odd */
3280
              Unit *quotlsu;    /* for save */
3281
              Int quotdigits;   /* .. */
3282
 
3283
              /* Fastpath when residue is truly 0 is worthwhile [and */
3284
              /* simplifies the code below] */
3285
              if (*var1 == 0 && var1units == 1)
3286
                {               /* residue is 0 */
3287
                  Int exp = lhs->exponent;      /* save min(exponents) */
3288
                  if (rhs->exponent < exp)
3289
                    exp = rhs->exponent;
3290
                  decNumberZero (res);  /* 0 coefficient */
3291
#if DECSUBSET
3292
                  if (set->extended)
3293
#endif
3294
                    res->exponent = exp;        /* .. with proper exponent */
3295
                  break;
3296
                }
3297
              /* note if the quotient was odd */
3298
              if (*accnext & 0x01)
3299
                wasodd = 1;     /* acc is odd */
3300
              quotlsu = accnext;        /* save in case need to reinspect */
3301
              quotdigits = accdigits;   /* .. */
3302
 
3303
              /* treat the residue, in var1, as the value to return, via acc */
3304
              /* calculate the unused zero digits.  This is the smaller of: */
3305
              /*   var1 initial padding (saved above) */
3306
              /*   var2 residual padding, which happens to be given by: */
3307
              postshift =
3308
                var1initpad + exponent - lhs->exponent + rhs->exponent;
3309
              /* [the 'exponent' term accounts for the shifts during divide] */
3310
              if (var1initpad < postshift)
3311
                postshift = var1initpad;
3312
 
3313
              /* shift var1 the requested amount, and adjust its digits */
3314
              var1units = decShiftToLeast (var1, var1units, postshift);
3315
              accnext = var1;
3316
              accdigits = decGetDigits (var1, var1units);
3317
              accunits = D2U (accdigits);
3318
 
3319
              exponent = lhs->exponent; /* exponent is smaller of lhs & rhs */
3320
              if (rhs->exponent < exponent)
3321
                exponent = rhs->exponent;
3322
              bits = lhs->bits; /* remainder sign is always as lhs */
3323
 
3324
              /* Now correct the result if we are doing remainderNear; if it */
3325
              /* (looking just at coefficients) is > rhs/2, or == rhs/2 and */
3326
              /* the integer was odd then the result should be rem-rhs. */
3327
              if (op & REMNEAR)
3328
                {
3329
                  Int compare, tarunits;        /* work */
3330
                  Unit *up;     /* .. */
3331
 
3332
 
3333
                  /* calculate remainder*2 into the var1 buffer (which has */
3334
                  /* 'headroom' of an extra unit and hence enough space) */
3335
                  /* [a dedicated 'double' loop would be faster, here] */
3336
                  tarunits =
3337
                    decUnitAddSub (accnext, accunits, accnext, accunits, 0,
3338
                                   accnext, 1);
3339
                  /* decDumpAr('r', accnext, tarunits); */
3340
 
3341
                  /* Here, accnext (var1) holds tarunits Units with twice the */
3342
                  /* remainder's coefficient, which we must now compare to the */
3343
                  /* RHS.  The remainder's exponent may be smaller than the RHS's. */
3344
                  compare =
3345
                    decUnitCompare (accnext, tarunits, rhs->lsu,
3346
                                    D2U (rhs->digits),
3347
                                    rhs->exponent - exponent);
3348
                  if (compare == BADINT)
3349
                    {           /* deep trouble */
3350
                      *status |= DEC_Insufficient_storage;
3351
                      break;
3352
                    }
3353
 
3354
                  /* now restore the remainder by dividing by two; we know the */
3355
                  /* lsu is even. */
3356
                  for (up = accnext; up < accnext + tarunits; up++)
3357
                    {
3358
                      Int half; /* half to add to lower unit */
3359
                      half = *up & 0x01;
3360
                      *up /= 2; /* [shift] */
3361
                      if (!half)
3362
                        continue;
3363
                      *(up - 1) += (DECDPUNMAX + 1) / 2;
3364
                    }
3365
                  /* [accunits still describes the original remainder length] */
3366
 
3367
                  if (compare > 0 || (compare == 0 && wasodd))
3368
                    {           /* adjustment needed */
3369
                      Int exp, expunits, exprem;        /* work */
3370
                      /* This is effectively causing round-up of the quotient, */
3371
                      /* so if it was the rare case where it was full and all */
3372
                      /* nines, it would overflow and hence division-impossible */
3373
                      /* should be raised */
3374
                      Flag allnines = 0; /* 1 if quotient all nines */
3375
                      if (quotdigits == reqdigits)
3376
                        {       /* could be borderline */
3377
                          for (up = quotlsu;; up++)
3378
                            {
3379
                              if (quotdigits > DECDPUN)
3380
                                {
3381
                                  if (*up != DECDPUNMAX)
3382
                                    break;      /* non-nines */
3383
                                }
3384
                              else
3385
                                {       /* this is the last Unit */
3386
                                  if (*up == powers[quotdigits] - 1)
3387
                                    allnines = 1;
3388
                                  break;
3389
                                }
3390
                              quotdigits -= DECDPUN;    /* checked those digits */
3391
                            }   /* up */
3392
                        }       /* borderline check */
3393
                      if (allnines)
3394
                        {
3395
                          *status |= DEC_Division_impossible;
3396
                          break;
3397
                        }
3398
 
3399
                      /* we need rem-rhs; the sign will invert.  Again we can */
3400
                      /* safely use var1 for the working Units array. */
3401
                      exp = rhs->exponent - exponent;   /* RHS padding needed */
3402
                      /* Calculate units and remainder from exponent. */
3403
                      expunits = exp / DECDPUN;
3404
                      exprem = exp % DECDPUN;
3405
                      /* subtract [A+B*(-m)]; the result will always be negative */
3406
                      accunits = -decUnitAddSub (accnext, accunits,
3407
                                                 rhs->lsu, D2U (rhs->digits),
3408
                                                 expunits, accnext,
3409
                                                 -(Int) powers[exprem]);
3410
                      accdigits = decGetDigits (accnext, accunits);     /* count digits exactly */
3411
                      accunits = D2U (accdigits);       /* and recalculate the units for copy */
3412
                      /* [exponent is as for original remainder] */
3413
                      bits ^= DECNEG;   /* flip the sign */
3414
                    }
3415
                }               /* REMNEAR */
3416
            }                   /* REMAINDER or REMNEAR */
3417
        }                       /* not DIVIDE */
3418
 
3419
      /* Set exponent and bits */
3420
      res->exponent = exponent;
3421
      res->bits = (uByte) (bits & DECNEG);      /* [cleaned] */
3422
 
3423
      /* Now the coefficient. */
3424
      decSetCoeff (res, set, accnext, accdigits, &residue, status);
3425
 
3426
      decFinish (res, set, &residue, status);   /* final cleanup */
3427
 
3428
#if DECSUBSET
3429
      /* If a divide then strip trailing zeros if subset [after round] */
3430
      if (!set->extended && (op == DIVIDE))
3431
        decTrim (res, 0, &dropped);
3432
#endif
3433
    }
3434
  while (0);                     /* end protected */
3435
 
3436
  if (varalloc != NULL)
3437
    free (varalloc);            /* drop any storage we used */
3438
  if (allocacc != NULL)
3439
    free (allocacc);            /* .. */
3440
  if (allocrhs != NULL)
3441
    free (allocrhs);            /* .. */
3442
  if (alloclhs != NULL)
3443
    free (alloclhs);            /* .. */
3444
  return res;
3445
}
3446
 
3447
/* ------------------------------------------------------------------ */
3448
/* decMultiplyOp -- multiplication operation                          */
3449
/*                                                                    */
3450
/*  This routine performs the multiplication C=A x B.                 */
3451
/*                                                                    */
3452
/*   res is C, the result.  C may be A and/or B (e.g., X=X*X)         */
3453
/*   lhs is A                                                         */
3454
/*   rhs is B                                                         */
3455
/*   set is the context                                               */
3456
/*   status is the usual accumulator                                  */
3457
/*                                                                    */
3458
/* C must have space for set->digits digits.                          */
3459
/*                                                                    */
3460
/* ------------------------------------------------------------------ */
3461
/* Note: We use 'long' multiplication rather than Karatsuba, as the   */
3462
/* latter would give only a minor improvement for the short numbers   */
3463
/* we expect to handle most (and uses much more memory).              */
3464
/*                                                                    */
3465
/* We always have to use a buffer for the accumulator.                */
3466
/* ------------------------------------------------------------------ */
3467
static decNumber *
3468
decMultiplyOp (decNumber * res, const decNumber * lhs,
3469
               const decNumber * rhs, decContext * set, uInt * status)
3470
{
3471
  decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
3472
  decNumber *allocrhs = NULL;   /* .., rhs */
3473
  Unit accbuff[D2U (DECBUFFER * 2 + 1)];        /* local buffer (+1 in case DECBUFFER==0) */
3474
  Unit *acc = accbuff;          /* -> accumulator array for exact result */
3475
  Unit *allocacc = NULL;        /* -> allocated buffer, iff allocated */
3476
  const Unit *mer, *mermsup;    /* work */
3477
  Int accunits;                 /* Units of accumulator in use */
3478
  Int madlength;                /* Units in multiplicand */
3479
  Int shift;                    /* Units to shift multiplicand by */
3480
  Int need;                     /* Accumulator units needed */
3481
  Int exponent;                 /* work */
3482
  Int residue = 0;               /* rounding residue */
3483
  uByte bits;                   /* result sign */
3484
  uByte merged;                 /* merged flags */
3485
 
3486
#if DECCHECK
3487
  if (decCheckOperands (res, lhs, rhs, set))
3488
    return res;
3489
#endif
3490
 
3491
  do
3492
    {                           /* protect allocated storage */
3493
#if DECSUBSET
3494
      if (!set->extended)
3495
        {
3496
          /* reduce operands and set lostDigits status, as needed */
3497
          if (lhs->digits > set->digits)
3498
            {
3499
              alloclhs = decRoundOperand (lhs, set, status);
3500
              if (alloclhs == NULL)
3501
                break;
3502
              lhs = alloclhs;
3503
            }
3504
          if (rhs->digits > set->digits)
3505
            {
3506
              allocrhs = decRoundOperand (rhs, set, status);
3507
              if (allocrhs == NULL)
3508
                break;
3509
              rhs = allocrhs;
3510
            }
3511
        }
3512
#endif
3513
      /* [following code does not require input rounding] */
3514
 
3515
      /* precalculate result sign */
3516
      bits = (uByte) ((lhs->bits ^ rhs->bits) & DECNEG);
3517
 
3518
      /* handle infinities and NaNs */
3519
      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
3520
      if (merged)
3521
        {                       /* a special bit set */
3522
          if (merged & (DECSNAN | DECNAN))
3523
            {                   /* one or two NaNs */
3524
              decNaNs (res, lhs, rhs, status);
3525
              break;
3526
            }
3527
          /* one or two infinities. Infinity * 0 is invalid */
3528
          if (((lhs->bits & DECSPECIAL) == 0 && ISZERO (lhs))
3529
              || ((rhs->bits & DECSPECIAL) == 0 && ISZERO (rhs)))
3530
            {
3531
              *status |= DEC_Invalid_operation;
3532
              break;
3533
            }
3534
          decNumberZero (res);
3535
          res->bits = bits | DECINF;    /* infinity */
3536
          break;
3537
        }
3538
 
3539
      /* For best speed, as in DMSRCN, we use the shorter number as the */
3540
      /* multiplier (rhs) and the longer as the multiplicand (lhs) */
3541
      if (lhs->digits < rhs->digits)
3542
        {                       /* swap... */
3543
          const decNumber *hold = lhs;
3544
          lhs = rhs;
3545
          rhs = hold;
3546
        }
3547
 
3548
      /* if accumulator is too long for local storage, then allocate */
3549
      need = D2U (lhs->digits) + D2U (rhs->digits);     /* maximum units in result */
3550
      if (need * sizeof (Unit) > sizeof (accbuff))
3551
        {
3552
          allocacc = (Unit *) malloc (need * sizeof (Unit));
3553
          if (allocacc == NULL)
3554
            {
3555
              *status |= DEC_Insufficient_storage;
3556
              break;
3557
            }
3558
          acc = allocacc;       /* use the allocated space */
3559
        }
3560
 
3561
      /* Now the main long multiplication loop */
3562
      /* Unlike the equivalent in the IBM Java implementation, there */
3563
      /* is no advantage in calculating from msu to lsu.  So we do it */
3564
      /* by the book, as it were. */
3565
      /* Each iteration calculates ACC=ACC+MULTAND*MULT */
3566
      accunits = 1;             /* accumulator starts at '0' */
3567
      *acc = 0;                  /* .. (lsu=0) */
3568
      shift = 0;         /* no multiplicand shift at first */
3569
      madlength = D2U (lhs->digits);    /* we know this won't change */
3570
      mermsup = rhs->lsu + D2U (rhs->digits);   /* -> msu+1 of multiplier */
3571
 
3572
      for (mer = rhs->lsu; mer < mermsup; mer++)
3573
        {
3574
          /* Here, *mer is the next Unit in the multiplier to use */
3575
          /* If non-zero [optimization] add it... */
3576
          if (*mer != 0)
3577
            {
3578
              accunits =
3579
                decUnitAddSub (&acc[shift], accunits - shift, lhs->lsu,
3580
                               madlength, 0, &acc[shift], *mer) + shift;
3581
            }
3582
          else
3583
            {                   /* extend acc with a 0; we'll use it shortly */
3584
              /* [this avoids length of <=0 later] */
3585
              *(acc + accunits) = 0;
3586
              accunits++;
3587
            }
3588
          /* multiply multiplicand by 10**DECDPUN for next Unit to left */
3589
          shift++;              /* add this for 'logical length' */
3590
        }                       /* n */
3591
#if DECTRACE
3592
      /* Show exact result */
3593
      decDumpAr ('*', acc, accunits);
3594
#endif
3595
 
3596
      /* acc now contains the exact result of the multiplication */
3597
      /* Build a decNumber from it, noting if any residue */
3598
      res->bits = bits;         /* set sign */
3599
      res->digits = decGetDigits (acc, accunits);       /* count digits exactly */
3600
 
3601
      /* We might have a 31-bit wrap in calculating the exponent. */
3602
      /* This can only happen if both input exponents are negative and */
3603
      /* both their magnitudes are large.  If we did wrap, we set a safe */
3604
      /* very negative exponent, from which decFinalize() will raise a */
3605
      /* hard underflow. */
3606
      exponent = lhs->exponent + rhs->exponent; /* calculate exponent */
3607
      if (lhs->exponent < 0 && rhs->exponent < 0 && exponent > 0)
3608
        exponent = -2 * DECNUMMAXE;     /* force underflow */
3609
      res->exponent = exponent; /* OK to overwrite now */
3610
 
3611
      /* Set the coefficient.  If any rounding, residue records */
3612
      decSetCoeff (res, set, acc, res->digits, &residue, status);
3613
 
3614
      decFinish (res, set, &residue, status);   /* final cleanup */
3615
    }
3616
  while (0);                     /* end protected */
3617
 
3618
  if (allocacc != NULL)
3619
    free (allocacc);            /* drop any storage we used */
3620
  if (allocrhs != NULL)
3621
    free (allocrhs);            /* .. */
3622
  if (alloclhs != NULL)
3623
    free (alloclhs);            /* .. */
3624
  return res;
3625
}
3626
 
3627
/* ------------------------------------------------------------------ */
3628
/* decQuantizeOp  -- force exponent to requested value                */
3629
/*                                                                    */
3630
/*   This computes C = op(A, B), where op adjusts the coefficient     */
3631
/*   of C (by rounding or shifting) such that the exponent (-scale)   */
3632
/*   of C has the value B or matches the exponent of B.               */
3633
/*   The numerical value of C will equal A, except for the effects of */
3634
/*   any rounding that occurred.                                      */
3635
/*                                                                    */
3636
/*   res is C, the result.  C may be A or B                           */
3637
/*   lhs is A, the number to adjust                                   */
3638
/*   rhs is B, the requested exponent                                 */
3639
/*   set is the context                                               */
3640
/*   quant is 1 for quantize or 0 for rescale                         */
3641
/*   status is the status accumulator (this can be called without     */
3642
/*          risk of control loss)                                     */
3643
/*                                                                    */
3644
/* C must have space for set->digits digits.                          */
3645
/*                                                                    */
3646
/* Unless there is an error or the result is infinite, the exponent   */
3647
/* after the operation is guaranteed to be that requested.            */
3648
/* ------------------------------------------------------------------ */
3649
static decNumber *
3650
decQuantizeOp (decNumber * res, const decNumber * lhs,
3651
               const decNumber * rhs, decContext * set, Flag quant, uInt * status)
3652
{
3653
  decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
3654
  decNumber *allocrhs = NULL;   /* .., rhs */
3655
  const decNumber *inrhs = rhs; /* save original rhs */
3656
  Int reqdigits = set->digits;  /* requested DIGITS */
3657
  Int reqexp;                   /* requested exponent [-scale] */
3658
  Int residue = 0;               /* rounding residue */
3659
  uByte merged;                 /* merged flags */
3660
  Int etiny = set->emin - (set->digits - 1);
3661
 
3662
#if DECCHECK
3663
  if (decCheckOperands (res, lhs, rhs, set))
3664
    return res;
3665
#endif
3666
 
3667
  do
3668
    {                           /* protect allocated storage */
3669
#if DECSUBSET
3670
      if (!set->extended)
3671
        {
3672
          /* reduce operands and set lostDigits status, as needed */
3673
          if (lhs->digits > reqdigits)
3674
            {
3675
              alloclhs = decRoundOperand (lhs, set, status);
3676
              if (alloclhs == NULL)
3677
                break;
3678
              lhs = alloclhs;
3679
            }
3680
          if (rhs->digits > reqdigits)
3681
            {                   /* [this only checks lostDigits] */
3682
              allocrhs = decRoundOperand (rhs, set, status);
3683
              if (allocrhs == NULL)
3684
                break;
3685
              rhs = allocrhs;
3686
            }
3687
        }
3688
#endif
3689
      /* [following code does not require input rounding] */
3690
 
3691
      /* Handle special values */
3692
      merged = (lhs->bits | rhs->bits) & DECSPECIAL;
3693
      if ((lhs->bits | rhs->bits) & DECSPECIAL)
3694
        {
3695
          /* NaNs get usual processing */
3696
          if (merged & (DECSNAN | DECNAN))
3697
            decNaNs (res, lhs, rhs, status);
3698
          /* one infinity but not both is bad */
3699
          else if ((lhs->bits ^ rhs->bits) & DECINF)
3700
            *status |= DEC_Invalid_operation;
3701
          /* both infinity: return lhs */
3702
          else
3703
            decNumberCopy (res, lhs);   /* [nop if in place] */
3704
          break;
3705
        }
3706
 
3707
      /* set requested exponent */
3708
      if (quant)
3709
        reqexp = inrhs->exponent;       /* quantize -- match exponents */
3710
      else
3711
        {                       /* rescale -- use value of rhs */
3712
          /* Original rhs must be an integer that fits and is in range */
3713
#if DECSUBSET
3714
          reqexp = decGetInt (inrhs, set);
3715
#else
3716
          reqexp = decGetInt (inrhs);
3717
#endif
3718
        }
3719
 
3720
#if DECSUBSET
3721
      if (!set->extended)
3722
        etiny = set->emin;      /* no subnormals */
3723
#endif
3724
 
3725
      if (reqexp == BADINT      /* bad (rescale only) or .. */
3726
          || (reqexp < etiny)   /* < lowest */
3727
          || (reqexp > set->emax))
3728
        {                       /* > Emax */
3729
          *status |= DEC_Invalid_operation;
3730
          break;
3731
        }
3732
 
3733
      /* we've processed the RHS, so we can overwrite it now if necessary */
3734
      if (ISZERO (lhs))
3735
        {                       /* zero coefficient unchanged */
3736
          decNumberCopy (res, lhs);     /* [nop if in place] */
3737
          res->exponent = reqexp;       /* .. just set exponent */
3738
#if DECSUBSET
3739
          if (!set->extended)
3740
            res->bits = 0;       /* subset specification; no -0 */
3741
#endif
3742
        }
3743
      else
3744
        {                       /* non-zero lhs */
3745
          Int adjust = reqexp - lhs->exponent;  /* digit adjustment needed */
3746
          /* if adjusted coefficient will not fit, give up now */
3747
          if ((lhs->digits - adjust) > reqdigits)
3748
            {
3749
              *status |= DEC_Invalid_operation;
3750
              break;
3751
            }
3752
 
3753
          if (adjust > 0)
3754
            {                   /* increasing exponent */
3755
              /* this will decrease the length of the coefficient by adjust */
3756
              /* digits, and must round as it does so */
3757
              decContext workset;       /* work */
3758
              workset = *set;   /* clone rounding, etc. */
3759
              workset.digits = lhs->digits - adjust;    /* set requested length */
3760
              /* [note that the latter can be <1, here] */
3761
              decCopyFit (res, lhs, &workset, &residue, status);        /* fit to result */
3762
              decApplyRound (res, &workset, residue, status);   /* .. and round */
3763
              residue = 0;       /* [used] */
3764
              /* If we rounded a 999s case, exponent will be off by one; */
3765
              /* adjust back if so. */
3766
              if (res->exponent > reqexp)
3767
                {
3768
                  res->digits = decShiftToMost (res->lsu, res->digits, 1);      /* shift */
3769
                  res->exponent--;      /* (re)adjust the exponent. */
3770
                }
3771
#if DECSUBSET
3772
              if (ISZERO (res) && !set->extended)
3773
                res->bits = 0;   /* subset; no -0 */
3774
#endif
3775
            }                   /* increase */
3776
          else                  /* adjust<=0 */
3777
            {                   /* decreasing or = exponent */
3778
              /* this will increase the length of the coefficient by -adjust */
3779
              /* digits, by adding trailing zeros. */
3780
              decNumberCopy (res, lhs); /* [it will fit] */
3781
              /* if padding needed (adjust<0), add it now... */
3782
              if (adjust < 0)
3783
                {
3784
                  res->digits =
3785
                    decShiftToMost (res->lsu, res->digits, -adjust);
3786
                  res->exponent += adjust;      /* adjust the exponent */
3787
                }
3788
            }                   /* decrease */
3789
        }                       /* non-zero */
3790
 
3791
      /* Check for overflow [do not use Finalize in this case, as an */
3792
      /* overflow here is a "don't fit" situation] */
3793
      if (res->exponent > set->emax - res->digits + 1)
3794
        {                       /* too big */
3795
          *status |= DEC_Invalid_operation;
3796
          break;
3797
        }
3798
      else
3799
        {
3800
          decFinalize (res, set, &residue, status);     /* set subnormal flags */
3801
          *status &= ~DEC_Underflow;    /* suppress Underflow [754r] */
3802
        }
3803
    }
3804
  while (0);                     /* end protected */
3805
 
3806
  if (allocrhs != NULL)
3807
    free (allocrhs);            /* drop any storage we used */
3808
  if (alloclhs != NULL)
3809
    free (alloclhs);            /* .. */
3810
  return res;
3811
}
3812
 
3813
/* ------------------------------------------------------------------ */
3814
/* decCompareOp -- compare, min, or max two Numbers                   */
3815
/*                                                                    */
3816
/*   This computes C = A ? B and returns the signum (as a Number)     */
3817
/*   for COMPARE or the maximum or minimum (for COMPMAX and COMPMIN). */
3818
/*                                                                    */
3819
/*   res is C, the result.  C may be A and/or B (e.g., X=X?X)         */
3820
/*   lhs is A                                                         */
3821
/*   rhs is B                                                         */
3822
/*   set is the context                                               */
3823
/*   op  is the operation flag                                        */
3824
/*   status is the usual accumulator                                  */
3825
/*                                                                    */
3826
/* C must have space for one digit for COMPARE or set->digits for     */
3827
/* COMPMAX and COMPMIN.                                               */
3828
/* ------------------------------------------------------------------ */
3829
/* The emphasis here is on speed for common cases, and avoiding       */
3830
/* coefficient comparison if possible.                                */
3831
/* ------------------------------------------------------------------ */
3832
decNumber *
3833
decCompareOp (decNumber * res, const decNumber * lhs, const decNumber * rhs,
3834
              decContext * set, Flag op, uInt * status)
3835
{
3836
  decNumber *alloclhs = NULL;   /* non-NULL if rounded lhs allocated */
3837
  decNumber *allocrhs = NULL;   /* .., rhs */
3838
  Int result = 0;                /* default result value */
3839
  uByte merged;                 /* merged flags */
3840
  uByte bits = 0;                /* non-0 for NaN */
3841
 
3842
#if DECCHECK
3843
  if (decCheckOperands (res, lhs, rhs, set))
3844
    return res;
3845
#endif
3846
 
3847
  do
3848
    {                           /* protect allocated storage */
3849
#if DECSUBSET
3850
      if (!set->extended)
3851
        {
3852
          /* reduce operands and set lostDigits status, as needed */
3853
          if (lhs->digits > set->digits)
3854
            {
3855
              alloclhs = decRoundOperand (lhs, set, status);
3856
              if (alloclhs == NULL)
3857
                {
3858
                  result = BADINT;
3859
                  break;
3860
                }
3861
              lhs = alloclhs;
3862
            }
3863
          if (rhs->digits > set->digits)
3864
            {
3865
              allocrhs = decRoundOperand (rhs, set, status);
3866
              if (allocrhs == NULL)
3867
                {
3868
                  result = BADINT;
3869
                  break;
3870
                }
3871
              rhs = allocrhs;
3872
            }
3873
        }
3874
#endif
3875
      /* [following code does not require input rounding] */
3876
 
3877
      /* handle NaNs now; let infinities drop through */
3878
      /* +++ review sNaN handling with 754r, for now assumes sNaN */
3879
      /* (even just one) leads to NaN. */
3880
      merged = (lhs->bits | rhs->bits) & (DECSNAN | DECNAN);
3881
      if (merged)
3882
        {                       /* a NaN bit set */
3883
          if (op == COMPARE);
3884
          else if (merged & DECSNAN);
3885
          else
3886
            {                   /* 754r rules for MIN and MAX ignore single NaN */
3887
              /* here if MIN or MAX, and one or two quiet NaNs */
3888
              if (lhs->bits & rhs->bits & DECNAN);
3889
              else
3890
                {               /* just one quiet NaN */
3891
                  /* force choice to be the non-NaN operand */
3892
                  op = COMPMAX;
3893
                  if (lhs->bits & DECNAN)
3894
                    result = -1;        /* pick rhs */
3895
                  else
3896
                    result = +1;        /* pick lhs */
3897
                  break;
3898
                }
3899
            }
3900
          op = COMPNAN;         /* use special path */
3901
          decNaNs (res, lhs, rhs, status);
3902
          break;
3903
        }
3904
 
3905
      result = decCompare (lhs, rhs);   /* we have numbers */
3906
    }
3907
  while (0);                     /* end protected */
3908
 
3909
  if (result == BADINT)
3910
    *status |= DEC_Insufficient_storage;        /* rare */
3911
  else
3912
    {
3913
      if (op == COMPARE)
3914
        {                       /* return signum */
3915
          decNumberZero (res);  /* [always a valid result] */
3916
          if (result == 0)
3917
            res->bits = bits;   /* (maybe qNaN) */
3918
          else
3919
            {
3920
              *res->lsu = 1;
3921
              if (result < 0)
3922
                res->bits = DECNEG;
3923
            }
3924
        }
3925
      else if (op == COMPNAN);  /* special, drop through */
3926
      else
3927
        {                       /* MAX or MIN, non-NaN result */
3928
          Int residue = 0;       /* rounding accumulator */
3929
          /* choose the operand for the result */
3930
          const decNumber *choice;
3931
          if (result == 0)
3932
            {                   /* operands are numerically equal */
3933
              /* choose according to sign then exponent (see 754r) */
3934
              uByte slhs = (lhs->bits & DECNEG);
3935
              uByte srhs = (rhs->bits & DECNEG);
3936
#if DECSUBSET
3937
              if (!set->extended)
3938
                {               /* subset: force left-hand */
3939
                  op = COMPMAX;
3940
                  result = +1;
3941
                }
3942
              else
3943
#endif
3944
              if (slhs != srhs)
3945
                {               /* signs differ */
3946
                  if (slhs)
3947
                    result = -1;        /* rhs is max */
3948
                  else
3949
                    result = +1;        /* lhs is max */
3950
                }
3951
              else if (slhs && srhs)
3952
                {               /* both negative */
3953
                  if (lhs->exponent < rhs->exponent)
3954
                    result = +1;
3955
                  else
3956
                    result = -1;
3957
                  /* [if equal, we use lhs, technically identical] */
3958
                }
3959
              else
3960
                {               /* both positive */
3961
                  if (lhs->exponent > rhs->exponent)
3962
                    result = +1;
3963
                  else
3964
                    result = -1;
3965
                  /* [ditto] */
3966
                }
3967
            }                   /* numerically equal */
3968
          /* here result will be non-0 */
3969
          if (op == COMPMIN)
3970
            result = -result;   /* reverse if looking for MIN */
3971
          choice = (result > 0 ? lhs : rhs);     /* choose */
3972
          /* copy chosen to result, rounding if need be */
3973
          decCopyFit (res, choice, set, &residue, status);
3974
          decFinish (res, set, &residue, status);
3975
        }
3976
    }
3977
  if (allocrhs != NULL)
3978
    free (allocrhs);            /* free any storage we used */
3979
  if (alloclhs != NULL)
3980
    free (alloclhs);            /* .. */
3981
  return res;
3982
}
3983
 
3984
/* ------------------------------------------------------------------ */
3985
/* decCompare -- compare two decNumbers by numerical value            */
3986
/*                                                                    */
3987
/*  This routine compares A ? B without altering them.                */
3988
/*                                                                    */
3989
/*  Arg1 is A, a decNumber which is not a NaN                         */
3990
/*  Arg2 is B, a decNumber which is not a NaN                         */
3991
/*                                                                    */
3992
/*  returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure   */
3993
/*  (the only possible failure is an allocation error)                */
3994
/* ------------------------------------------------------------------ */
3995
/* This could be merged into decCompareOp */
3996
static Int
3997
decCompare (const decNumber * lhs, const decNumber * rhs)
3998
{
3999
  Int result;                   /* result value */
4000
  Int sigr;                     /* rhs signum */
4001
  Int compare;                  /* work */
4002
  result = 1;                   /* assume signum(lhs) */
4003
  if (ISZERO (lhs))
4004
    result = 0;
4005
  else if (decNumberIsNegative (lhs))
4006
    result = -1;
4007
  sigr = 1;                     /* compute signum(rhs) */
4008
  if (ISZERO (rhs))
4009
    sigr = 0;
4010
  else if (decNumberIsNegative (rhs))
4011
    sigr = -1;
4012
  if (result > sigr)
4013
    return +1;                  /* L > R, return 1 */
4014
  if (result < sigr)
4015
    return -1;                  /* R < L, return -1 */
4016
 
4017
  /* signums are the same */
4018
  if (result == 0)
4019
    return 0;                    /* both 0 */
4020
  /* Both non-zero */
4021
  if ((lhs->bits | rhs->bits) & DECINF)
4022
    {                           /* one or more infinities */
4023
      if (lhs->bits == rhs->bits)
4024
        result = 0;              /* both the same */
4025
      else if (decNumberIsInfinite (rhs))
4026
        result = -result;
4027
      return result;
4028
    }
4029
 
4030
  /* we must compare the coefficients, allowing for exponents */
4031
  if (lhs->exponent > rhs->exponent)
4032
    {                           /* LHS exponent larger */
4033
      /* swap sides, and sign */
4034
      const decNumber *temp = lhs;
4035
      lhs = rhs;
4036
      rhs = temp;
4037
      result = -result;
4038
    }
4039
 
4040
  compare = decUnitCompare (lhs->lsu, D2U (lhs->digits),
4041
                            rhs->lsu, D2U (rhs->digits),
4042
                            rhs->exponent - lhs->exponent);
4043
 
4044
  if (compare != BADINT)
4045
    compare *= result;          /* comparison succeeded */
4046
  return compare;               /* what we got */
4047
}
4048
 
4049
/* ------------------------------------------------------------------ */
4050
/* decUnitCompare -- compare two >=0 integers in Unit arrays          */
4051
/*                                                                    */
4052
/*  This routine compares A ? B*10**E where A and B are unit arrays   */
4053
/*  A is a plain integer                                              */
4054
/*  B has an exponent of E (which must be non-negative)               */
4055
/*                                                                    */
4056
/*  Arg1 is A first Unit (lsu)                                        */
4057
/*  Arg2 is A length in Units                                         */
4058
/*  Arg3 is B first Unit (lsu)                                        */
4059
/*  Arg4 is B length in Units                                         */
4060
/*  Arg5 is E                                                         */
4061
/*                                                                    */
4062
/*  returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure   */
4063
/*  (the only possible failure is an allocation error)                */
4064
/* ------------------------------------------------------------------ */
4065
static Int
4066
decUnitCompare (const Unit * a, Int alength, const Unit * b, Int blength, Int exp)
4067
{
4068
  Unit *acc;                    /* accumulator for result */
4069
  Unit accbuff[D2U (DECBUFFER + 1)];    /* local buffer */
4070
  Unit *allocacc = NULL;        /* -> allocated acc buffer, iff allocated */
4071
  Int accunits, need;           /* units in use or needed for acc */
4072
  const Unit *l, *r, *u;        /* work */
4073
  Int expunits, exprem, result; /* .. */
4074
 
4075
  if (exp == 0)
4076
    {                           /* aligned; fastpath */
4077
      if (alength > blength)
4078
        return 1;
4079
      if (alength < blength)
4080
        return -1;
4081
      /* same number of units in both -- need unit-by-unit compare */
4082
      l = a + alength - 1;
4083
      r = b + alength - 1;
4084
      for (; l >= a; l--, r--)
4085
        {
4086
          if (*l > *r)
4087
            return 1;
4088
          if (*l < *r)
4089
            return -1;
4090
        }
4091
      return 0;                  /* all units match */
4092
    }                           /* aligned */
4093
 
4094
  /* Unaligned.  If one is >1 unit longer than the other, padded */
4095
  /* approximately, then we can return easily */
4096
  if (alength > blength + (Int) D2U (exp))
4097
    return 1;
4098
  if (alength + 1 < blength + (Int) D2U (exp))
4099
    return -1;
4100
 
4101
  /* We need to do a real subtract.  For this, we need a result buffer */
4102
  /* even though we only are interested in the sign.  Its length needs */
4103
  /* to be the larger of alength and padded blength, +2 */
4104
  need = blength + D2U (exp);   /* maximum real length of B */
4105
  if (need < alength)
4106
    need = alength;
4107
  need += 2;
4108
  acc = accbuff;                /* assume use local buffer */
4109
  if (need * sizeof (Unit) > sizeof (accbuff))
4110
    {
4111
      allocacc = (Unit *) malloc (need * sizeof (Unit));
4112
      if (allocacc == NULL)
4113
        return BADINT;          /* hopeless -- abandon */
4114
      acc = allocacc;
4115
    }
4116
  /* Calculate units and remainder from exponent. */
4117
  expunits = exp / DECDPUN;
4118
  exprem = exp % DECDPUN;
4119
  /* subtract [A+B*(-m)] */
4120
  accunits = decUnitAddSub (a, alength, b, blength, expunits, acc,
4121
                            -(Int) powers[exprem]);
4122
  /* [UnitAddSub result may have leading zeros, even on zero] */
4123
  if (accunits < 0)
4124
    result = -1;                /* negative result */
4125
  else
4126
    {                           /* non-negative result */
4127
      /* check units of the result before freeing any storage */
4128
      for (u = acc; u < acc + accunits - 1 && *u == 0;)
4129
        u++;
4130
      result = (*u == 0 ? 0 : +1);
4131
    }
4132
  /* clean up and return the result */
4133
  if (allocacc != NULL)
4134
    free (allocacc);            /* drop any storage we used */
4135
  return result;
4136
}
4137
 
4138
/* ------------------------------------------------------------------ */
4139
/* decUnitAddSub -- add or subtract two >=0 integers in Unit arrays   */
4140
/*                                                                    */
4141
/*  This routine performs the calculation:                            */
4142
/*                                                                    */
4143
/*  C=A+(B*M)                                                         */
4144
/*                                                                    */
4145
/*  Where M is in the range -DECDPUNMAX through +DECDPUNMAX.          */
4146
/*                                                                    */
4147
/*  A may be shorter or longer than B.                                */
4148
/*                                                                    */
4149
/*  Leading zeros are not removed after a calculation.  The result is */
4150
/*  either the same length as the longer of A and B (adding any       */
4151
/*  shift), or one Unit longer than that (if a Unit carry occurred).  */
4152
/*                                                                    */
4153
/*  A and B content are not altered unless C is also A or B.          */
4154
/*  C may be the same array as A or B, but only if no zero padding is */
4155
/*  requested (that is, C may be B only if bshift==0).                */
4156
/*  C is filled from the lsu; only those units necessary to complete  */
4157
/*  the calculation are referenced.                                   */
4158
/*                                                                    */
4159
/*  Arg1 is A first Unit (lsu)                                        */
4160
/*  Arg2 is A length in Units                                         */
4161
/*  Arg3 is B first Unit (lsu)                                        */
4162
/*  Arg4 is B length in Units                                         */
4163
/*  Arg5 is B shift in Units  (>=0; pads with 0 units if positive)    */
4164
/*  Arg6 is C first Unit (lsu)                                        */
4165
/*  Arg7 is M, the multiplier                                         */
4166
/*                                                                    */
4167
/*  returns the count of Units written to C, which will be non-zero   */
4168
/*  and negated if the result is negative.  That is, the sign of the  */
4169
/*  returned Int is the sign of the result (positive for zero) and    */
4170
/*  the absolute value of the Int is the count of Units.              */
4171
/*                                                                    */
4172
/*  It is the caller's responsibility to make sure that C size is     */
4173
/*  safe, allowing space if necessary for a one-Unit carry.           */
4174
/*                                                                    */
4175
/*  This routine is severely performance-critical; *any* change here  */
4176
/*  must be measured (timed) to assure no performance degradation.    */
4177
/*  In particular, trickery here tends to be counter-productive, as   */
4178
/*  increased complexity of code hurts register optimizations on      */
4179
/*  register-poor architectures.  Avoiding divisions is nearly        */
4180
/*  always a Good Idea, however.                                      */
4181
/*                                                                    */
4182
/* Special thanks to Rick McGuire (IBM Cambridge, MA) and Dave Clark  */
4183
/* (IBM Warwick, UK) for some of the ideas used in this routine.      */
4184
/* ------------------------------------------------------------------ */
4185
static Int
4186
decUnitAddSub (const Unit * a, Int alength,
4187
               const Unit * b, Int blength, Int bshift, Unit * c, Int m)
4188
{
4189
  const Unit *alsu = a;         /* A lsu [need to remember it] */
4190
  Unit *clsu = c;               /* C ditto */
4191
  Unit *minC;                   /* low water mark for C */
4192
  Unit *maxC;                   /* high water mark for C */
4193
  eInt carry = 0;                /* carry integer (could be Long) */
4194
  Int add;                      /* work */
4195
#if DECDPUN==4                  /* myriadal */
4196
  Int est;                      /* estimated quotient */
4197
#endif
4198
 
4199
#if DECTRACE
4200
  if (alength < 1 || blength < 1)
4201
    printf ("decUnitAddSub: alen blen m %d %d [%d]\n", alength, blength, m);
4202
#endif
4203
 
4204
  maxC = c + alength;           /* A is usually the longer */
4205
  minC = c + blength;           /* .. and B the shorter */
4206
  if (bshift != 0)
4207
    {                           /* B is shifted; low As copy across */
4208
      minC += bshift;
4209
      /* if in place [common], skip copy unless there's a gap [rare] */
4210
      if (a == c && bshift <= alength)
4211
        {
4212
          c += bshift;
4213
          a += bshift;
4214
        }
4215
      else
4216
        for (; c < clsu + bshift; a++, c++)
4217
          {                     /* copy needed */
4218
            if (a < alsu + alength)
4219
              *c = *a;
4220
            else
4221
              *c = 0;
4222
          }
4223
    }
4224
  if (minC > maxC)
4225
    {                           /* swap */
4226
      Unit *hold = minC;
4227
      minC = maxC;
4228
      maxC = hold;
4229
    }
4230
 
4231
  /* For speed, we do the addition as two loops; the first where both A */
4232
  /* and B contribute, and the second (if necessary) where only one or */
4233
  /* other of the numbers contribute. */
4234
  /* Carry handling is the same (i.e., duplicated) in each case. */
4235
  for (; c < minC; c++)
4236
    {
4237
      carry += *a;
4238
      a++;
4239
      carry += ((eInt) * b) * m;        /* [special-casing m=1/-1 */
4240
      b++;                      /* here is not a win] */
4241
      /* here carry is new Unit of digits; it could be +ve or -ve */
4242
      if ((ueInt) carry <= DECDPUNMAX)
4243
        {                       /* fastpath 0-DECDPUNMAX */
4244
          *c = (Unit) carry;
4245
          carry = 0;
4246
          continue;
4247
        }
4248
      /* remainder operator is undefined if negative, so we must test */
4249
#if DECDPUN==4                  /* use divide-by-multiply */
4250
      if (carry >= 0)
4251
        {
4252
          est = (((ueInt) carry >> 11) * 53687) >> 18;
4253
          *c = (Unit) (carry - est * (DECDPUNMAX + 1)); /* remainder */
4254
          carry = est;          /* likely quotient [89%] */
4255
          if (*c < DECDPUNMAX + 1)
4256
            continue;           /* estimate was correct */
4257
          carry++;
4258
          *c -= DECDPUNMAX + 1;
4259
          continue;
4260
        }
4261
      /* negative case */
4262
      carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);       /* make positive */
4263
      est = (((ueInt) carry >> 11) * 53687) >> 18;
4264
      *c = (Unit) (carry - est * (DECDPUNMAX + 1));
4265
      carry = est - (DECDPUNMAX + 1);   /* correctly negative */
4266
      if (*c < DECDPUNMAX + 1)
4267
        continue;               /* was OK */
4268
      carry++;
4269
      *c -= DECDPUNMAX + 1;
4270
#else
4271
      if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
4272
        {                       /* fastpath carry +1 */
4273
          *c = (Unit) (carry - (DECDPUNMAX + 1));       /* [helps additions] */
4274
          carry = 1;
4275
          continue;
4276
        }
4277
      if (carry >= 0)
4278
        {
4279
          *c = (Unit) (carry % (DECDPUNMAX + 1));
4280
          carry = carry / (DECDPUNMAX + 1);
4281
          continue;
4282
        }
4283
      /* negative case */
4284
      carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);       /* make positive */
4285
      *c = (Unit) (carry % (DECDPUNMAX + 1));
4286
      carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
4287
#endif
4288
    }                           /* c */
4289
 
4290
  /* we now may have one or other to complete */
4291
  /* [pretest to avoid loop setup/shutdown] */
4292
  if (c < maxC)
4293
    for (; c < maxC; c++)
4294
      {
4295
        if (a < alsu + alength)
4296
          {                     /* still in A */
4297
            carry += *a;
4298
            a++;
4299
          }
4300
        else
4301
          {                     /* inside B */
4302
            carry += ((eInt) * b) * m;
4303
            b++;
4304
          }
4305
        /* here carry is new Unit of digits; it could be +ve or -ve and */
4306
        /* magnitude up to DECDPUNMAX squared */
4307
        if ((ueInt) carry <= DECDPUNMAX)
4308
          {                     /* fastpath 0-DECDPUNMAX */
4309
            *c = (Unit) carry;
4310
            carry = 0;
4311
            continue;
4312
          }
4313
        /* result for this unit is negative or >DECDPUNMAX */
4314
#if DECDPUN==4                  /* use divide-by-multiply */
4315
        /* remainder is undefined if negative, so we must test */
4316
        if (carry >= 0)
4317
          {
4318
            est = (((ueInt) carry >> 11) * 53687) >> 18;
4319
            *c = (Unit) (carry - est * (DECDPUNMAX + 1));       /* remainder */
4320
            carry = est;        /* likely quotient [79.7%] */
4321
            if (*c < DECDPUNMAX + 1)
4322
              continue;         /* estimate was correct */
4323
            carry++;
4324
            *c -= DECDPUNMAX + 1;
4325
            continue;
4326
          }
4327
        /* negative case */
4328
        carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);     /* make positive */
4329
        est = (((ueInt) carry >> 11) * 53687) >> 18;
4330
        *c = (Unit) (carry - est * (DECDPUNMAX + 1));
4331
        carry = est - (DECDPUNMAX + 1); /* correctly negative */
4332
        if (*c < DECDPUNMAX + 1)
4333
          continue;             /* was OK */
4334
        carry++;
4335
        *c -= DECDPUNMAX + 1;
4336
#else
4337
        if ((ueInt) carry < (DECDPUNMAX + 1) * 2)
4338
          {                     /* fastpath carry 1 */
4339
            *c = (Unit) (carry - (DECDPUNMAX + 1));
4340
            carry = 1;
4341
            continue;
4342
          }
4343
        /* remainder is undefined if negative, so we must test */
4344
        if (carry >= 0)
4345
          {
4346
            *c = (Unit) (carry % (DECDPUNMAX + 1));
4347
            carry = carry / (DECDPUNMAX + 1);
4348
            continue;
4349
          }
4350
        /* negative case */
4351
        carry = carry + (eInt) (DECDPUNMAX + 1) * (DECDPUNMAX + 1);     /* make positive */
4352
        *c = (Unit) (carry % (DECDPUNMAX + 1));
4353
        carry = carry / (DECDPUNMAX + 1) - (DECDPUNMAX + 1);
4354
#endif
4355
      }                         /* c */
4356
 
4357
  /* OK, all A and B processed; might still have carry or borrow */
4358
  /* return number of Units in the result, negated if a borrow */
4359
  if (carry == 0)
4360
    return c - clsu;            /* no carry, we're done */
4361
  if (carry > 0)
4362
    {                           /* positive carry */
4363
      *c = (Unit) carry;        /* place as new unit */
4364
      c++;                      /* .. */
4365
      return c - clsu;
4366
    }
4367
  /* -ve carry: it's a borrow; complement needed */
4368
  add = 1;                      /* temporary carry... */
4369
  for (c = clsu; c < maxC; c++)
4370
    {
4371
      add = DECDPUNMAX + add - *c;
4372
      if (add <= DECDPUNMAX)
4373
        {
4374
          *c = (Unit) add;
4375
          add = 0;
4376
        }
4377
      else
4378
        {
4379
          *c = 0;
4380
          add = 1;
4381
        }
4382
    }
4383
  /* add an extra unit iff it would be non-zero */
4384
#if DECTRACE
4385
  printf ("UAS borrow: add %d, carry %d\n", add, carry);
4386
#endif
4387
  if ((add - carry - 1) != 0)
4388
    {
4389
      *c = (Unit) (add - carry - 1);
4390
      c++;                      /* interesting, include it */
4391
    }
4392
  return clsu - c;              /* -ve result indicates borrowed */
4393
}
4394
 
4395
/* ------------------------------------------------------------------ */
4396
/* decTrim -- trim trailing zeros or normalize                        */
4397
/*                                                                    */
4398
/*   dn is the number to trim or normalize                            */
4399
/*   all is 1 to remove all trailing zeros, 0 for just fraction ones  */
4400
/*   dropped returns the number of discarded trailing zeros           */
4401
/*   returns dn                                                       */
4402
/*                                                                    */
4403
/* All fields are updated as required.  This is a utility operation,  */
4404
/* so special values are unchanged and no error is possible.          */
4405
/* ------------------------------------------------------------------ */
4406
static decNumber *
4407
decTrim (decNumber * dn, Flag all, Int * dropped)
4408
{
4409
  Int d, exp;                   /* work */
4410
  uInt cut;                     /* .. */
4411
  Unit *up;                     /* -> current Unit */
4412
 
4413
#if DECCHECK
4414
  if (decCheckOperands (dn, DECUNUSED, DECUNUSED, DECUNUSED))
4415
    return dn;
4416
#endif
4417
 
4418
  *dropped = 0;                  /* assume no zeros dropped */
4419
  if ((dn->bits & DECSPECIAL)   /* fast exit if special .. */
4420
      || (*dn->lsu & 0x01))
4421
    return dn;                  /* .. or odd */
4422
  if (ISZERO (dn))
4423
    {                           /* .. or 0 */
4424
      dn->exponent = 0;          /* (sign is preserved) */
4425
      return dn;
4426
    }
4427
 
4428
  /* we have a finite number which is even */
4429
  exp = dn->exponent;
4430
  cut = 1;                      /* digit (1-DECDPUN) in Unit */
4431
  up = dn->lsu;                 /* -> current Unit */
4432
  for (d = 0; d < dn->digits - 1; d++)
4433
    {                           /* [don't strip the final digit] */
4434
      /* slice by powers */
4435
#if DECDPUN<=4
4436
      uInt quot = QUOT10 (*up, cut);
4437
      if ((*up - quot * powers[cut]) != 0)
4438
        break;                  /* found non-0 digit */
4439
#else
4440
      if (*up % powers[cut] != 0)
4441
        break;                  /* found non-0 digit */
4442
#endif
4443
      /* have a trailing 0 */
4444
      if (!all)
4445
        {                       /* trimming */
4446
          /* [if exp>0 then all trailing 0s are significant for trim] */
4447
          if (exp <= 0)
4448
            {                   /* if digit might be significant */
4449
              if (exp == 0)
4450
                break;          /* then quit */
4451
              exp++;            /* next digit might be significant */
4452
            }
4453
        }
4454
      cut++;                    /* next power */
4455
      if (cut > DECDPUN)
4456
        {                       /* need new Unit */
4457
          up++;
4458
          cut = 1;
4459
        }
4460
    }                           /* d */
4461
  if (d == 0)
4462
    return dn;                  /* none dropped */
4463
 
4464
  /* effect the drop */
4465
  decShiftToLeast (dn->lsu, D2U (dn->digits), d);
4466
  dn->exponent += d;            /* maintain numerical value */
4467
  dn->digits -= d;              /* new length */
4468
  *dropped = d;                 /* report the count */
4469
  return dn;
4470
}
4471
 
4472
/* ------------------------------------------------------------------ */
4473
/* decShiftToMost -- shift digits in array towards most significant   */
4474
/*                                                                    */
4475
/*   uar    is the array                                              */
4476
/*   digits is the count of digits in use in the array                */
4477
/*   shift  is the number of zeros to pad with (least significant);   */
4478
/*     it must be zero or positive                                    */
4479
/*                                                                    */
4480
/*   returns the new length of the integer in the array, in digits    */
4481
/*                                                                    */
4482
/* No overflow is permitted (that is, the uar array must be known to  */
4483
/* be large enough to hold the result, after shifting).               */
4484
/* ------------------------------------------------------------------ */
4485
static Int
4486
decShiftToMost (Unit * uar, Int digits, Int shift)
4487
{
4488
  Unit *target, *source, *first;        /* work */
4489
  uInt rem;                     /* for division */
4490
  Int cut;                      /* odd 0's to add */
4491
  uInt next;                    /* work */
4492
 
4493
  if (shift == 0)
4494
    return digits;              /* [fastpath] nothing to do */
4495
  if ((digits + shift) <= DECDPUN)
4496
    {                           /* [fastpath] single-unit case */
4497
      *uar = (Unit) (*uar * powers[shift]);
4498
      return digits + shift;
4499
    }
4500
 
4501
  cut = (DECDPUN - shift % DECDPUN) % DECDPUN;
4502
  source = uar + D2U (digits) - 1;      /* where msu comes from */
4503
  first = uar + D2U (digits + shift) - 1;       /* where msu of source will end up */
4504
  target = source + D2U (shift);        /* where upper part of first cut goes */
4505
  next = 0;
4506
 
4507
  for (; source >= uar; source--, target--)
4508
    {
4509
      /* split the source Unit and accumulate remainder for next */
4510
#if DECDPUN<=4
4511
      uInt quot = QUOT10 (*source, cut);
4512
      rem = *source - quot * powers[cut];
4513
      next += quot;
4514
#else
4515
      rem = *source % powers[cut];
4516
      next += *source / powers[cut];
4517
#endif
4518
      if (target <= first)
4519
        *target = (Unit) next;  /* write to target iff valid */
4520
      next = rem * powers[DECDPUN - cut];       /* save remainder for next Unit */
4521
    }
4522
  /* propagate to one below and clear the rest */
4523
  for (; target >= uar; target--)
4524
    {
4525
      *target = (Unit) next;
4526
      next = 0;
4527
    }
4528
  return digits + shift;
4529
}
4530
 
4531
/* ------------------------------------------------------------------ */
4532
/* decShiftToLeast -- shift digits in array towards least significant */
4533
/*                                                                    */
4534
/*   uar   is the array                                               */
4535
/*   units is length of the array, in units                           */
4536
/*   shift is the number of digits to remove from the lsu end; it     */
4537
/*     must be zero or positive and less than units*DECDPUN.          */
4538
/*                                                                    */
4539
/*   returns the new length of the integer in the array, in units     */
4540
/*                                                                    */
4541
/* Removed digits are discarded (lost).  Units not required to hold   */
4542
/* the final result are unchanged.                                    */
4543
/* ------------------------------------------------------------------ */
4544
static Int
4545
decShiftToLeast (Unit * uar, Int units, Int shift)
4546
{
4547
  Unit *target, *up;            /* work */
4548
  Int cut, count;               /* work */
4549
  Int quot, rem;                /* for division */
4550
 
4551
  if (shift == 0)
4552
    return units;               /* [fastpath] nothing to do */
4553
 
4554
  up = uar + shift / DECDPUN;   /* source; allow for whole Units */
4555
  cut = shift % DECDPUN;        /* odd 0's to drop */
4556
  target = uar;                 /* both paths */
4557
  if (cut == 0)
4558
    {                           /* whole units shift */
4559
      for (; up < uar + units; target++, up++)
4560
        *target = *up;
4561
      return target - uar;
4562
    }
4563
  /* messier */
4564
  count = units * DECDPUN - shift;      /* the maximum new length */
4565
#if DECDPUN<=4
4566
  quot = QUOT10 (*up, cut);
4567
#else
4568
  quot = *up / powers[cut];
4569
#endif
4570
  for (;; target++)
4571
    {
4572
      *target = (Unit) quot;
4573
      count -= (DECDPUN - cut);
4574
      if (count <= 0)
4575
        break;
4576
      up++;
4577
      quot = *up;
4578
#if DECDPUN<=4
4579
      quot = QUOT10 (quot, cut);
4580
      rem = *up - quot * powers[cut];
4581
#else
4582
      rem = quot % powers[cut];
4583
      quot = quot / powers[cut];
4584
#endif
4585
      *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
4586
      count -= cut;
4587
      if (count <= 0)
4588
        break;
4589
    }
4590
  return target - uar + 1;
4591
}
4592
 
4593
#if DECSUBSET
4594
/* ------------------------------------------------------------------ */
4595
/* decRoundOperand -- round an operand  [used for subset only]        */
4596
/*                                                                    */
4597
/*   dn is the number to round (dn->digits is > set->digits)          */
4598
/*   set is the relevant context                                      */
4599
/*   status is the status accumulator                                 */
4600
/*                                                                    */
4601
/*   returns an allocated decNumber with the rounded result.          */
4602
/*                                                                    */
4603
/* lostDigits and other status may be set by this.                    */
4604
/*                                                                    */
4605
/* Since the input is an operand, we are not permitted to modify it.  */
4606
/* We therefore return an allocated decNumber, rounded as required.   */
4607
/* It is the caller's responsibility to free the allocated storage.   */
4608
/*                                                                    */
4609
/* If no storage is available then the result cannot be used, so NULL */
4610
/* is returned.                                                       */
4611
/* ------------------------------------------------------------------ */
4612
static decNumber *
4613
decRoundOperand (const decNumber * dn, decContext * set, uInt * status)
4614
{
4615
  decNumber *res;               /* result structure */
4616
  uInt newstatus = 0;            /* status from round */
4617
  Int residue = 0;               /* rounding accumulator */
4618
 
4619
  /* Allocate storage for the returned decNumber, big enough for the */
4620
  /* length specified by the context */
4621
  res = (decNumber *) malloc (sizeof (decNumber)
4622
                              + (D2U (set->digits) - 1) * sizeof (Unit));
4623
  if (res == NULL)
4624
    {
4625
      *status |= DEC_Insufficient_storage;
4626
      return NULL;
4627
    }
4628
  decCopyFit (res, dn, set, &residue, &newstatus);
4629
  decApplyRound (res, set, residue, &newstatus);
4630
 
4631
  /* If that set Inexact then we "lost digits" */
4632
  if (newstatus & DEC_Inexact)
4633
    newstatus |= DEC_Lost_digits;
4634
  *status |= newstatus;
4635
  return res;
4636
}
4637
#endif
4638
 
4639
/* ------------------------------------------------------------------ */
4640
/* decCopyFit -- copy a number, shortening the coefficient if needed  */
4641
/*                                                                    */
4642
/*   dest is the target decNumber                                     */
4643
/*   src  is the source decNumber                                     */
4644
/*   set is the context [used for length (digits) and rounding mode]  */
4645
/*   residue is the residue accumulator                               */
4646
/*   status contains the current status to be updated                 */
4647
/*                                                                    */
4648
/* (dest==src is allowed and will be a no-op if fits)                 */
4649
/* All fields are updated as required.                                */
4650
/* ------------------------------------------------------------------ */
4651
static void
4652
decCopyFit (decNumber * dest, const decNumber * src, decContext * set,
4653
            Int * residue, uInt * status)
4654
{
4655
  dest->bits = src->bits;
4656
  dest->exponent = src->exponent;
4657
  decSetCoeff (dest, set, src->lsu, src->digits, residue, status);
4658
}
4659
 
4660
/* ------------------------------------------------------------------ */
4661
/* decSetCoeff -- set the coefficient of a number                     */
4662
/*                                                                    */
4663
/*   dn    is the number whose coefficient array is to be set.        */
4664
/*         It must have space for set->digits digits                  */
4665
/*   set   is the context [for size]                                  */
4666
/*   lsu   -> lsu of the source coefficient [may be dn->lsu]          */
4667
/*   len   is digits in the source coefficient [may be dn->digits]    */
4668
/*   residue is the residue accumulator.  This has values as in       */
4669
/*         decApplyRound, and will be unchanged unless the            */
4670
/*         target size is less than len.  In this case, the           */
4671
/*         coefficient is truncated and the residue is updated to     */
4672
/*         reflect the previous residue and the dropped digits.       */
4673
/*   status is the status accumulator, as usual                       */
4674
/*                                                                    */
4675
/* The coefficient may already be in the number, or it can be an      */
4676
/* external intermediate array.  If it is in the number, lsu must ==  */
4677
/* dn->lsu and len must == dn->digits.                                */
4678
/*                                                                    */
4679
/* Note that the coefficient length (len) may be < set->digits, and   */
4680
/* in this case this merely copies the coefficient (or is a no-op     */
4681
/* if dn->lsu==lsu).                                                  */
4682
/*                                                                    */
4683
/* Note also that (only internally, from decNumberRescale and         */
4684
/* decSetSubnormal) the value of set->digits may be less than one,    */
4685
/* indicating a round to left.                                        */
4686
/* This routine handles that case correctly; caller ensures space.    */
4687
/*                                                                    */
4688
/* dn->digits, dn->lsu (and as required), and dn->exponent are        */
4689
/* updated as necessary.   dn->bits (sign) is unchanged.              */
4690
/*                                                                    */
4691
/* DEC_Rounded status is set if any digits are discarded.             */
4692
/* DEC_Inexact status is set if any non-zero digits are discarded, or */
4693
/*                       incoming residue was non-0 (implies rounded) */
4694
/* ------------------------------------------------------------------ */
4695
/* mapping array: maps 0-9 to canonical residues, so that we can */
4696
/* adjust by a residue in range [-1, +1] and achieve correct rounding */
4697
/*                             0  1  2  3  4  5  6  7  8  9 */
4698
static const uByte resmap[10] = { 0, 3, 3, 3, 3, 5, 7, 7, 7, 7 };
4699
static void
4700
decSetCoeff (decNumber * dn, decContext * set, const Unit * lsu,
4701
             Int len, Int * residue, uInt * status)
4702
{
4703
  Int discard;                  /* number of digits to discard */
4704
  uInt discard1;                /* first discarded digit */
4705
  uInt cut;                     /* cut point in Unit */
4706
  uInt quot, rem;               /* for divisions */
4707
  Unit *target;                 /* work */
4708
  const Unit *up;               /* work */
4709
  Int count;                    /* .. */
4710
#if DECDPUN<=4
4711
  uInt temp;                    /* .. */
4712
#endif
4713
 
4714
  discard = len - set->digits;  /* digits to discard */
4715
  if (discard <= 0)
4716
    {                           /* no digits are being discarded */
4717
      if (dn->lsu != lsu)
4718
        {                       /* copy needed */
4719
          /* copy the coefficient array to the result number; no shift needed */
4720
          up = lsu;
4721
          for (target = dn->lsu; target < dn->lsu + D2U (len); target++, up++)
4722
            {
4723
              *target = *up;
4724
            }
4725
          dn->digits = len;     /* set the new length */
4726
        }
4727
      /* dn->exponent and residue are unchanged */
4728
      if (*residue != 0)
4729
        *status |= (DEC_Inexact | DEC_Rounded); /* record inexactitude */
4730
      return;
4731
    }
4732
 
4733
  /* we have to discard some digits */
4734
  *status |= DEC_Rounded;       /* accumulate Rounded status */
4735
  if (*residue > 1)
4736
    *residue = 1;               /* previous residue now to right, so -1 to +1 */
4737
 
4738
  if (discard > len)
4739
    {                           /* everything, +1, is being discarded */
4740
      /* guard digit is 0 */
4741
      /* residue is all the number [NB could be all 0s] */
4742
      if (*residue <= 0)
4743
        for (up = lsu + D2U (len) - 1; up >= lsu; up--)
4744
          {
4745
            if (*up != 0)
4746
              {                 /* found a non-0 */
4747
                *residue = 1;
4748
                break;          /* no need to check any others */
4749
              }
4750
          }
4751
      if (*residue != 0)
4752
        *status |= DEC_Inexact; /* record inexactitude */
4753
      *dn->lsu = 0;              /* coefficient will now be 0 */
4754
      dn->digits = 1;           /* .. */
4755
      dn->exponent += discard;  /* maintain numerical value */
4756
      return;
4757
    }                           /* total discard */
4758
 
4759
  /* partial discard [most common case] */
4760
  /* here, at least the first (most significant) discarded digit exists */
4761
 
4762
  /* spin up the number, noting residue as we pass, until we get to */
4763
  /* the Unit with the first discarded digit.  When we get there, */
4764
  /* extract it and remember where we're at */
4765
  count = 0;
4766
  for (up = lsu;; up++)
4767
    {
4768
      count += DECDPUN;
4769
      if (count >= discard)
4770
        break;                  /* full ones all checked */
4771
      if (*up != 0)
4772
        *residue = 1;
4773
    }                           /* up */
4774
 
4775
  /* here up -> Unit with discarded digit */
4776
  cut = discard - (count - DECDPUN) - 1;
4777
  if (cut == DECDPUN - 1)
4778
    {                           /* discard digit is at top */
4779
#if DECDPUN<=4
4780
      discard1 = QUOT10 (*up, DECDPUN - 1);
4781
      rem = *up - discard1 * powers[DECDPUN - 1];
4782
#else
4783
      rem = *up % powers[DECDPUN - 1];
4784
      discard1 = *up / powers[DECDPUN - 1];
4785
#endif
4786
      if (rem != 0)
4787
        *residue = 1;
4788
      up++;                     /* move to next */
4789
      cut = 0;                   /* bottom digit of result */
4790
      quot = 0;                  /* keep a certain compiler happy */
4791
    }
4792
  else
4793
    {
4794
      /* discard digit is in low digit(s), not top digit */
4795
      if (cut == 0)
4796
        quot = *up;
4797
      else                      /* cut>0 */
4798
        {                       /* it's not at bottom of Unit */
4799
#if DECDPUN<=4
4800
          quot = QUOT10 (*up, cut);
4801
          rem = *up - quot * powers[cut];
4802
#else
4803
          rem = *up % powers[cut];
4804
          quot = *up / powers[cut];
4805
#endif
4806
          if (rem != 0)
4807
            *residue = 1;
4808
        }
4809
      /* discard digit is now at bottom of quot */
4810
#if DECDPUN<=4
4811
      temp = (quot * 6554) >> 16;       /* fast /10 */
4812
      /* Vowels algorithm here not a win (9 instructions) */
4813
      discard1 = quot - X10 (temp);
4814
      quot = temp;
4815
#else
4816
      discard1 = quot % 10;
4817
      quot = quot / 10;
4818
#endif
4819
      cut++;                    /* update cut */
4820
    }
4821
 
4822
  /* here: up -> Unit of the array with discarded digit */
4823
  /*       cut is the division point for each Unit */
4824
  /*       quot holds the uncut high-order digits for the current */
4825
  /*            Unit, unless cut==0 in which case it's still in *up */
4826
  /* copy the coefficient array to the result number, shifting as we go */
4827
  count = set->digits;          /* digits to end up with */
4828
  if (count <= 0)
4829
    {                           /* special for Rescale/Subnormal :-( */
4830
      *dn->lsu = 0;              /* .. result is 0 */
4831
      dn->digits = 1;           /* .. */
4832
    }
4833
  else
4834
    {                           /* shift to least */
4835
      /* [this is similar to decShiftToLeast code, with copy] */
4836
      dn->digits = count;       /* set the new length */
4837
      if (cut == 0)
4838
        {
4839
          /* on unit boundary, so simple shift down copy loop suffices */
4840
          for (target = dn->lsu; target < dn->lsu + D2U (count);
4841
               target++, up++)
4842
            {
4843
              *target = *up;
4844
            }
4845
        }
4846
      else
4847
        for (target = dn->lsu;; target++)
4848
          {
4849
            *target = (Unit) quot;
4850
            count -= (DECDPUN - cut);
4851
            if (count <= 0)
4852
              break;
4853
            up++;
4854
            quot = *up;
4855
#if DECDPUN<=4
4856
            quot = QUOT10 (quot, cut);
4857
            rem = *up - quot * powers[cut];
4858
#else
4859
            rem = quot % powers[cut];
4860
            quot = quot / powers[cut];
4861
#endif
4862
            *target = (Unit) (*target + rem * powers[DECDPUN - cut]);
4863
            count -= cut;
4864
            if (count <= 0)
4865
              break;
4866
          }
4867
    }                           /* shift to least needed */
4868
  dn->exponent += discard;      /* maintain numerical value */
4869
 
4870
  /* here, discard1 is the guard digit, and residue is everything else */
4871
  /* [use mapping to accumulate residue safely] */
4872
  *residue += resmap[discard1];
4873
 
4874
  if (*residue != 0)
4875
    *status |= DEC_Inexact;     /* record inexactitude */
4876
  return;
4877
}
4878
 
4879
/* ------------------------------------------------------------------ */
4880
/* decApplyRound -- apply pending rounding to a number                */
4881
/*                                                                    */
4882
/*   dn    is the number, with space for set->digits digits           */
4883
/*   set   is the context [for size and rounding mode]                */
4884
/*   residue indicates pending rounding, being any accumulated        */
4885
/*         guard and sticky information.  It may be:                  */
4886
/*         6-9: rounding digit is >5                                  */
4887
/*         5:   rounding digit is exactly half-way                    */
4888
/*         1-4: rounding digit is <5 and >0                           */
4889
/*         0:   the coefficient is exact                              */
4890
/*        -1:   as 1, but the hidden digits are subtractive, that     */
4891
/*              is, of the opposite sign to dn.  In this case the     */
4892
/*              coefficient must be non-0.                            */
4893
/*   status is the status accumulator, as usual                       */
4894
/*                                                                    */
4895
/* This routine applies rounding while keeping the length of the      */
4896
/* coefficient constant.  The exponent and status are unchanged       */
4897
/* except if:                                                         */
4898
/*                                                                    */
4899
/*   -- the coefficient was increased and is all nines (in which      */
4900
/*      case Overflow could occur, and is handled directly here so    */
4901
/*      the caller does not need to re-test for overflow)             */
4902
/*                                                                    */
4903
/*   -- the coefficient was decreased and becomes all nines (in which */
4904
/*      case Underflow could occur, and is also handled directly).    */
4905
/*                                                                    */
4906
/* All fields in dn are updated as required.                          */
4907
/*                                                                    */
4908
/* ------------------------------------------------------------------ */
4909
static void
4910
decApplyRound (decNumber * dn, decContext * set, Int residue, uInt * status)
4911
{
4912
  Int bump;                     /* 1 if coefficient needs to be incremented */
4913
  /* -1 if coefficient needs to be decremented */
4914
 
4915
  if (residue == 0)
4916
    return;                     /* nothing to apply */
4917
 
4918
  bump = 0;                      /* assume a smooth ride */
4919
 
4920
  /* now decide whether, and how, to round, depending on mode */
4921
  switch (set->round)
4922
    {
4923
    case DEC_ROUND_DOWN:
4924
      {
4925
        /* no change, except if negative residue */
4926
        if (residue < 0)
4927
          bump = -1;
4928
        break;
4929
      }                         /* r-d */
4930
 
4931
    case DEC_ROUND_HALF_DOWN:
4932
      {
4933
        if (residue > 5)
4934
          bump = 1;
4935
        break;
4936
      }                         /* r-h-d */
4937
 
4938
    case DEC_ROUND_HALF_EVEN:
4939
      {
4940
        if (residue > 5)
4941
          bump = 1;             /* >0.5 goes up */
4942
        else if (residue == 5)
4943
          {                     /* exactly 0.5000... */
4944
            /* 0.5 goes up iff [new] lsd is odd */
4945
            if (*dn->lsu & 0x01)
4946
              bump = 1;
4947
          }
4948
        break;
4949
      }                         /* r-h-e */
4950
 
4951
    case DEC_ROUND_HALF_UP:
4952
      {
4953
        if (residue >= 5)
4954
          bump = 1;
4955
        break;
4956
      }                         /* r-h-u */
4957
 
4958
    case DEC_ROUND_UP:
4959
      {
4960
        if (residue > 0)
4961
          bump = 1;
4962
        break;
4963
      }                         /* r-u */
4964
 
4965
    case DEC_ROUND_CEILING:
4966
      {
4967
        /* same as _UP for positive numbers, and as _DOWN for negatives */
4968
        /* [negative residue cannot occur on 0] */
4969
        if (decNumberIsNegative (dn))
4970
          {
4971
            if (residue < 0)
4972
              bump = -1;
4973
          }
4974
        else
4975
          {
4976
            if (residue > 0)
4977
              bump = 1;
4978
          }
4979
        break;
4980
      }                         /* r-c */
4981
 
4982
    case DEC_ROUND_FLOOR:
4983
      {
4984
        /* same as _UP for negative numbers, and as _DOWN for positive */
4985
        /* [negative residue cannot occur on 0] */
4986
        if (!decNumberIsNegative (dn))
4987
          {
4988
            if (residue < 0)
4989
              bump = -1;
4990
          }
4991
        else
4992
          {
4993
            if (residue > 0)
4994
              bump = 1;
4995
          }
4996
        break;
4997
      }                         /* r-f */
4998
 
4999
    default:
5000
      {                         /* e.g., DEC_ROUND_MAX */
5001
        *status |= DEC_Invalid_context;
5002
#if DECTRACE
5003
        printf ("Unknown rounding mode: %d\n", set->round);
5004
#endif
5005
        break;
5006
      }
5007
    }                           /* switch */
5008
 
5009
  /* now bump the number, up or down, if need be */
5010
  if (bump == 0)
5011
    return;                     /* no action required */
5012
 
5013
  /* Simply use decUnitAddSub unless we are bumping up and the number */
5014
  /* is all nines.  In this special case we set to 1000... and adjust */
5015
  /* the exponent by one (as otherwise we could overflow the array) */
5016
  /* Similarly handle all-nines result if bumping down. */
5017
  if (bump > 0)
5018
    {
5019
      Unit *up;                 /* work */
5020
      uInt count = dn->digits;  /* digits to be checked */
5021
      for (up = dn->lsu;; up++)
5022
        {
5023
          if (count <= DECDPUN)
5024
            {
5025
              /* this is the last Unit (the msu) */
5026
              if (*up != powers[count] - 1)
5027
                break;          /* not still 9s */
5028
              /* here if it, too, is all nines */
5029
              *up = (Unit) powers[count - 1];   /* here 999 -> 100 etc. */
5030
              for (up = up - 1; up >= dn->lsu; up--)
5031
                *up = 0; /* others all to 0 */
5032
              dn->exponent++;   /* and bump exponent */
5033
              /* [which, very rarely, could cause Overflow...] */
5034
              if ((dn->exponent + dn->digits) > set->emax + 1)
5035
                {
5036
                  decSetOverflow (dn, set, status);
5037
                }
5038
              return;           /* done */
5039
            }
5040
          /* a full unit to check, with more to come */
5041
          if (*up != DECDPUNMAX)
5042
            break;              /* not still 9s */
5043
          count -= DECDPUN;
5044
        }                       /* up */
5045
    }                           /* bump>0 */
5046
  else
5047
    {                           /* -1 */
5048
      /* here we are lookng for a pre-bump of 1000... (leading 1, */
5049
      /* all other digits zero) */
5050
      Unit *up, *sup;           /* work */
5051
      uInt count = dn->digits;  /* digits to be checked */
5052
      for (up = dn->lsu;; up++)
5053
        {
5054
          if (count <= DECDPUN)
5055
            {
5056
              /* this is the last Unit (the msu) */
5057
              if (*up != powers[count - 1])
5058
                break;          /* not 100.. */
5059
              /* here if we have the 1000... case */
5060
              sup = up;         /* save msu pointer */
5061
              *up = (Unit) powers[count] - 1;   /* here 100 in msu -> 999 */
5062
              /* others all to all-nines, too */
5063
              for (up = up - 1; up >= dn->lsu; up--)
5064
                *up = (Unit) powers[DECDPUN] - 1;
5065
              dn->exponent--;   /* and bump exponent */
5066
 
5067
              /* iff the number was at the subnormal boundary (exponent=etiny) */
5068
              /* then the exponent is now out of range, so it will in fact get */
5069
              /* clamped to etiny and the final 9 dropped. */
5070
              /* printf(">> emin=%d exp=%d sdig=%d\n", set->emin, */
5071
              /*        dn->exponent, set->digits); */
5072
              if (dn->exponent + 1 == set->emin - set->digits + 1)
5073
                {
5074
                  if (count == 1 && dn->digits == 1)
5075
                    *sup = 0;    /* here 9 -> 0[.9] */
5076
                  else
5077
                    {
5078
                      *sup = (Unit) powers[count - 1] - 1;      /* here 999.. in msu -> 99.. */
5079
                      dn->digits--;
5080
                    }
5081
                  dn->exponent++;
5082
                  *status |=
5083
                    DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5084
                }
5085
              return;           /* done */
5086
            }
5087
 
5088
          /* a full unit to check, with more to come */
5089
          if (*up != 0)
5090
            break;              /* not still 0s */
5091
          count -= DECDPUN;
5092
        }                       /* up */
5093
 
5094
    }                           /* bump<0 */
5095
 
5096
  /* Actual bump needed.  Do it. */
5097
  decUnitAddSub (dn->lsu, D2U (dn->digits), one, 1, 0, dn->lsu, bump);
5098
}
5099
 
5100
#if DECSUBSET
5101
/* ------------------------------------------------------------------ */
5102
/* decFinish -- finish processing a number                            */
5103
/*                                                                    */
5104
/*   dn is the number                                                 */
5105
/*   set is the context                                               */
5106
/*   residue is the rounding accumulator (as in decApplyRound)        */
5107
/*   status is the accumulator                                        */
5108
/*                                                                    */
5109
/* This finishes off the current number by:                           */
5110
/*    1. If not extended:                                             */
5111
/*       a. Converting a zero result to clean '0'                     */
5112
/*       b. Reducing positive exponents to 0, if would fit in digits  */
5113
/*    2. Checking for overflow and subnormals (always)                */
5114
/* Note this is just Finalize when no subset arithmetic.              */
5115
/* All fields are updated as required.                                */
5116
/* ------------------------------------------------------------------ */
5117
static void
5118
decFinish (decNumber * dn, decContext * set, Int * residue, uInt * status)
5119
{
5120
  if (!set->extended)
5121
    {
5122
      if ISZERO
5123
        (dn)
5124
        {                       /* value is zero */
5125
          dn->exponent = 0;      /* clean exponent .. */
5126
          dn->bits = 0;          /* .. and sign */
5127
          return;               /* no error possible */
5128
        }
5129
      if (dn->exponent >= 0)
5130
        {                       /* non-negative exponent */
5131
          /* >0; reduce to integer if possible */
5132
          if (set->digits >= (dn->exponent + dn->digits))
5133
            {
5134
              dn->digits = decShiftToMost (dn->lsu, dn->digits, dn->exponent);
5135
              dn->exponent = 0;
5136
            }
5137
        }
5138
    }                           /* !extended */
5139
 
5140
  decFinalize (dn, set, residue, status);
5141
}
5142
#endif
5143
 
5144
/* ------------------------------------------------------------------ */
5145
/* decFinalize -- final check, clamp, and round of a number           */
5146
/*                                                                    */
5147
/*   dn is the number                                                 */
5148
/*   set is the context                                               */
5149
/*   residue is the rounding accumulator (as in decApplyRound)        */
5150
/*   status is the status accumulator                                 */
5151
/*                                                                    */
5152
/* This finishes off the current number by checking for subnormal     */
5153
/* results, applying any pending rounding, checking for overflow,     */
5154
/* and applying any clamping.                                         */
5155
/* Underflow and overflow conditions are raised as appropriate.       */
5156
/* All fields are updated as required.                                */
5157
/* ------------------------------------------------------------------ */
5158
static void
5159
decFinalize (decNumber * dn, decContext * set, Int * residue, uInt * status)
5160
{
5161
  Int shift;                    /* shift needed if clamping */
5162
 
5163
  /* We have to be careful when checking the exponent as the adjusted */
5164
  /* exponent could overflow 31 bits [because it may already be up */
5165
  /* to twice the expected]. */
5166
 
5167
  /* First test for subnormal.  This must be done before any final */
5168
  /* round as the result could be rounded to Nmin or 0. */
5169
  if (dn->exponent < 0           /* negative exponent */
5170
      && (dn->exponent < set->emin - dn->digits + 1))
5171
    {
5172
      /* Go handle subnormals; this will apply round if needed. */
5173
      decSetSubnormal (dn, set, residue, status);
5174
      return;
5175
    }
5176
 
5177
  /* now apply any pending round (this could raise overflow). */
5178
  if (*residue != 0)
5179
    decApplyRound (dn, set, *residue, status);
5180
 
5181
  /* Check for overflow [redundant in the 'rare' case] or clamp */
5182
  if (dn->exponent <= set->emax - set->digits + 1)
5183
    return;                     /* neither needed */
5184
 
5185
  /* here when we might have an overflow or clamp to do */
5186
  if (dn->exponent > set->emax - dn->digits + 1)
5187
    {                           /* too big */
5188
      decSetOverflow (dn, set, status);
5189
      return;
5190
    }
5191
  /* here when the result is normal but in clamp range */
5192
  if (!set->clamp)
5193
    return;
5194
 
5195
  /* here when we need to apply the IEEE exponent clamp (fold-down) */
5196
  shift = dn->exponent - (set->emax - set->digits + 1);
5197
 
5198
  /* shift coefficient (if non-zero) */
5199
  if (!ISZERO (dn))
5200
    {
5201
      dn->digits = decShiftToMost (dn->lsu, dn->digits, shift);
5202
    }
5203
  dn->exponent -= shift;        /* adjust the exponent to match */
5204
  *status |= DEC_Clamped;       /* and record the dirty deed */
5205
  return;
5206
}
5207
 
5208
/* ------------------------------------------------------------------ */
5209
/* decSetOverflow -- set number to proper overflow value              */
5210
/*                                                                    */
5211
/*   dn is the number (used for sign [only] and result)               */
5212
/*   set is the context [used for the rounding mode]                  */
5213
/*   status contains the current status to be updated                 */
5214
/*                                                                    */
5215
/* This sets the sign of a number and sets its value to either        */
5216
/* Infinity or the maximum finite value, depending on the sign of     */
5217
/* dn and therounding mode, following IEEE 854 rules.                 */
5218
/* ------------------------------------------------------------------ */
5219
static void
5220
decSetOverflow (decNumber * dn, decContext * set, uInt * status)
5221
{
5222
  Flag needmax = 0;              /* result is maximum finite value */
5223
  uByte sign = dn->bits & DECNEG;       /* clean and save sign bit */
5224
 
5225
  if (ISZERO (dn))
5226
    {                           /* zero does not overflow magnitude */
5227
      Int emax = set->emax;     /* limit value */
5228
      if (set->clamp)
5229
        emax -= set->digits - 1;        /* lower if clamping */
5230
      if (dn->exponent > emax)
5231
        {                       /* clamp required */
5232
          dn->exponent = emax;
5233
          *status |= DEC_Clamped;
5234
        }
5235
      return;
5236
    }
5237
 
5238
  decNumberZero (dn);
5239
  switch (set->round)
5240
    {
5241
    case DEC_ROUND_DOWN:
5242
      {
5243
        needmax = 1;            /* never Infinity */
5244
        break;
5245
      }                         /* r-d */
5246
    case DEC_ROUND_CEILING:
5247
      {
5248
        if (sign)
5249
          needmax = 1;          /* Infinity if non-negative */
5250
        break;
5251
      }                         /* r-c */
5252
    case DEC_ROUND_FLOOR:
5253
      {
5254
        if (!sign)
5255
          needmax = 1;          /* Infinity if negative */
5256
        break;
5257
      }                         /* r-f */
5258
    default:
5259
      break;                    /* Infinity in all other cases */
5260
    }
5261
  if (needmax)
5262
    {
5263
      Unit *up;                 /* work */
5264
      Int count = set->digits;  /* nines to add */
5265
      dn->digits = count;
5266
      /* fill in all nines to set maximum value */
5267
      for (up = dn->lsu;; up++)
5268
        {
5269
          if (count > DECDPUN)
5270
            *up = DECDPUNMAX;   /* unit full o'nines */
5271
          else
5272
            {                   /* this is the msu */
5273
              *up = (Unit) (powers[count] - 1);
5274
              break;
5275
            }
5276
          count -= DECDPUN;     /* we filled those digits */
5277
        }                       /* up */
5278
      dn->bits = sign;          /* sign */
5279
      dn->exponent = set->emax - set->digits + 1;
5280
    }
5281
  else
5282
    dn->bits = sign | DECINF;   /* Value is +/-Infinity */
5283
  *status |= DEC_Overflow | DEC_Inexact | DEC_Rounded;
5284
}
5285
 
5286
/* ------------------------------------------------------------------ */
5287
/* decSetSubnormal -- process value whose exponent is <Emin           */
5288
/*                                                                    */
5289
/*   dn is the number (used as input as well as output; it may have   */
5290
/*         an allowed subnormal value, which may need to be rounded)  */
5291
/*   set is the context [used for the rounding mode]                  */
5292
/*   residue is any pending residue                                   */
5293
/*   status contains the current status to be updated                 */
5294
/*                                                                    */
5295
/* If subset mode, set result to zero and set Underflow flags.        */
5296
/*                                                                    */
5297
/* Value may be zero with a low exponent; this does not set Subnormal */
5298
/* but the exponent will be clamped to Etiny.                         */
5299
/*                                                                    */
5300
/* Otherwise ensure exponent is not out of range, and round as        */
5301
/* necessary.  Underflow is set if the result is Inexact.             */
5302
/* ------------------------------------------------------------------ */
5303
static void
5304
decSetSubnormal (decNumber * dn, decContext * set,
5305
                 Int * residue, uInt * status)
5306
{
5307
  decContext workset;           /* work */
5308
  Int etiny, adjust;            /* .. */
5309
 
5310
#if DECSUBSET
5311
  /* simple set to zero and 'hard underflow' for subset */
5312
  if (!set->extended)
5313
    {
5314
      decNumberZero (dn);
5315
      /* always full overflow */
5316
      *status |= DEC_Underflow | DEC_Subnormal | DEC_Inexact | DEC_Rounded;
5317
      return;
5318
    }
5319
#endif
5320
 
5321
  /* Full arithmetic -- allow subnormals, rounded to minimum exponent */
5322
  /* (Etiny) if needed */
5323
  etiny = set->emin - (set->digits - 1);        /* smallest allowed exponent */
5324
 
5325
  if ISZERO
5326
    (dn)
5327
    {                           /* value is zero */
5328
      /* residue can never be non-zero here */
5329
#if DECCHECK
5330
      if (*residue != 0)
5331
        {
5332
          printf ("++ Subnormal 0 residue %d\n", *residue);
5333
          *status |= DEC_Invalid_operation;
5334
        }
5335
#endif
5336
      if (dn->exponent < etiny)
5337
        {                       /* clamp required */
5338
          dn->exponent = etiny;
5339
          *status |= DEC_Clamped;
5340
        }
5341
      return;
5342
    }
5343
 
5344
  *status |= DEC_Subnormal;     /* we have a non-zero subnormal */
5345
 
5346
  adjust = etiny - dn->exponent;        /* calculate digits to remove */
5347
  if (adjust <= 0)
5348
    {                           /* not out of range; unrounded */
5349
      /* residue can never be non-zero here, so fast-path out */
5350
#if DECCHECK
5351
      if (*residue != 0)
5352
        {
5353
          printf ("++ Subnormal no-adjust residue %d\n", *residue);
5354
          *status |= DEC_Invalid_operation;
5355
        }
5356
#endif
5357
      /* it may already be inexact (from setting the coefficient) */
5358
      if (*status & DEC_Inexact)
5359
        *status |= DEC_Underflow;
5360
      return;
5361
    }
5362
 
5363
  /* adjust>0.  we need to rescale the result so exponent becomes Etiny */
5364
  /* [this code is similar to that in rescale] */
5365
  workset = *set;               /* clone rounding, etc. */
5366
  workset.digits = dn->digits - adjust; /* set requested length */
5367
  workset.emin -= adjust;       /* and adjust emin to match */
5368
  /* [note that the latter can be <1, here, similar to Rescale case] */
5369
  decSetCoeff (dn, &workset, dn->lsu, dn->digits, residue, status);
5370
  decApplyRound (dn, &workset, *residue, status);
5371
 
5372
  /* Use 754R/854 default rule: Underflow is set iff Inexact */
5373
  /* [independent of whether trapped] */
5374
  if (*status & DEC_Inexact)
5375
    *status |= DEC_Underflow;
5376
 
5377
  /* if we rounded up a 999s case, exponent will be off by one; adjust */
5378
  /* back if so [it will fit, because we shortened] */
5379
  if (dn->exponent > etiny)
5380
    {
5381
      dn->digits = decShiftToMost (dn->lsu, dn->digits, 1);
5382
      dn->exponent--;           /* (re)adjust the exponent. */
5383
    }
5384
}
5385
 
5386
/* ------------------------------------------------------------------ */
5387
/* decGetInt -- get integer from a number                             */
5388
/*                                                                    */
5389
/*   dn is the number [which will not be altered]                     */
5390
/*   set is the context [requested digits], subset only               */
5391
/*   returns the converted integer, or BADINT if error                */
5392
/*                                                                    */
5393
/* This checks and gets a whole number from the input decNumber.      */
5394
/* The magnitude of the integer must be <2^31.                        */
5395
/* Any discarded fractional part must be 0.                           */
5396
/* If subset it must also fit in set->digits                          */
5397
/* ------------------------------------------------------------------ */
5398
#if DECSUBSET
5399
static Int
5400
decGetInt (const decNumber * dn, decContext * set)
5401
{
5402
#else
5403
static Int
5404
decGetInt (const decNumber * dn)
5405
{
5406
#endif
5407
  Int theInt;                   /* result accumulator */
5408
  const Unit *up;               /* work */
5409
  Int got;                      /* digits (real or not) processed */
5410
  Int ilength = dn->digits + dn->exponent;      /* integral length */
5411
 
5412
  /* The number must be an integer that fits in 10 digits */
5413
  /* Assert, here, that 10 is enough for any rescale Etiny */
5414
#if DEC_MAX_EMAX > 999999999
5415
#error GetInt may need updating [for Emax]
5416
#endif
5417
#if DEC_MIN_EMIN < -999999999
5418
#error GetInt may need updating [for Emin]
5419
#endif
5420
  if (ISZERO (dn))
5421
    return 0;                    /* zeros are OK, with any exponent */
5422
  if (ilength > 10)
5423
    return BADINT;              /* always too big */
5424
#if DECSUBSET
5425
  if (!set->extended && ilength > set->digits)
5426
    return BADINT;
5427
#endif
5428
 
5429
  up = dn->lsu;                 /* ready for lsu */
5430
  theInt = 0;                    /* ready to accumulate */
5431
  if (dn->exponent >= 0)
5432
    {                           /* relatively easy */
5433
      /* no fractional part [usual]; allow for positive exponent */
5434
      got = dn->exponent;
5435
    }
5436
  else
5437
    {                           /* -ve exponent; some fractional part to check and discard */
5438
      Int count = -dn->exponent;        /* digits to discard */
5439
      /* spin up whole units until we get to the Unit with the unit digit */
5440
      for (; count >= DECDPUN; up++)
5441
        {
5442
          if (*up != 0)
5443
            return BADINT;      /* non-zero Unit to discard */
5444
          count -= DECDPUN;
5445
        }
5446
      if (count == 0)
5447
        got = 0;         /* [a multiple of DECDPUN] */
5448
      else
5449
        {                       /* [not multiple of DECDPUN] */
5450
          Int rem;              /* work */
5451
          /* slice off fraction digits and check for non-zero */
5452
#if DECDPUN<=4
5453
          theInt = QUOT10 (*up, count);
5454
          rem = *up - theInt * powers[count];
5455
#else
5456
          rem = *up % powers[count];    /* slice off discards */
5457
          theInt = *up / powers[count];
5458
#endif
5459
          if (rem != 0)
5460
            return BADINT;      /* non-zero fraction */
5461
          /* OK, we're good */
5462
          got = DECDPUN - count;        /* number of digits so far */
5463
          up++;                 /* ready for next */
5464
        }
5465
    }
5466
  /* collect the rest */
5467
  for (; got < ilength; up++)
5468
    {
5469
      theInt += *up * powers[got];
5470
      got += DECDPUN;
5471
    }
5472
  if ((ilength == 10)           /* check no wrap */
5473
      && (theInt / (Int) powers[got - DECDPUN] != *(up - 1)))
5474
    return BADINT;
5475
  /* [that test also disallows the BADINT result case] */
5476
 
5477
  /* apply any sign and return */
5478
  if (decNumberIsNegative (dn))
5479
    theInt = -theInt;
5480
  return theInt;
5481
}
5482
 
5483
/* ------------------------------------------------------------------ */
5484
/* decStrEq -- caseless comparison of strings                         */
5485
/*                                                                    */
5486
/*   str1 is one of the strings to compare                            */
5487
/*   str2 is the other                                                */
5488
/*                                                                    */
5489
/*   returns 1 if strings caseless-compare equal, 0 otherwise         */
5490
/*                                                                    */
5491
/* Note that the strings must be the same length if they are to       */
5492
/* compare equal; there is no padding.                                */
5493
/* ------------------------------------------------------------------ */
5494
/* [strcmpi is not in ANSI C] */
5495
static Flag
5496
decStrEq (const char *str1, const char *str2)
5497
{
5498
  for (;; str1++, str2++)
5499
    {
5500
      unsigned char u1 = (unsigned char) *str1;
5501
      unsigned char u2 = (unsigned char) *str2;
5502
      if (u1 == u2)
5503
        {
5504
          if (u1 == '\0')
5505
            break;
5506
        }
5507
      else
5508
        {
5509
          if (tolower (u1) != tolower (u2))
5510
            return 0;
5511
        }
5512
    }                           /* stepping */
5513
  return 1;
5514
}
5515
 
5516
/* ------------------------------------------------------------------ */
5517
/* decNaNs -- handle NaN operand or operands                          */
5518
/*                                                                    */
5519
/*   res    is the result number                                      */
5520
/*   lhs    is the first operand                                      */
5521
/*   rhs    is the second operand, or NULL if none                    */
5522
/*   status contains the current status                               */
5523
/*   returns res in case convenient                                   */
5524
/*                                                                    */
5525
/* Called when one or both operands is a NaN, and propagates the      */
5526
/* appropriate result to res.  When an sNaN is found, it is changed   */
5527
/* to a qNaN and Invalid operation is set.                            */
5528
/* ------------------------------------------------------------------ */
5529
static decNumber *
5530
decNaNs (decNumber * res, const decNumber * lhs, const decNumber * rhs, uInt * status)
5531
{
5532
  /* This decision tree ends up with LHS being the source pointer, */
5533
  /* and status updated if need be */
5534
  if (lhs->bits & DECSNAN)
5535
    *status |= DEC_Invalid_operation | DEC_sNaN;
5536
  else if (rhs == NULL);
5537
  else if (rhs->bits & DECSNAN)
5538
    {
5539
      lhs = rhs;
5540
      *status |= DEC_Invalid_operation | DEC_sNaN;
5541
    }
5542
  else if (lhs->bits & DECNAN);
5543
  else
5544
    lhs = rhs;
5545
 
5546
  decNumberCopy (res, lhs);
5547
  res->bits &= ~DECSNAN;        /* convert any sNaN to NaN, while */
5548
  res->bits |= DECNAN;          /* .. preserving sign */
5549
  res->exponent = 0;             /* clean exponent */
5550
  /* [coefficient was copied] */
5551
  return res;
5552
}
5553
 
5554
/* ------------------------------------------------------------------ */
5555
/* decStatus -- apply non-zero status                                 */
5556
/*                                                                    */
5557
/*   dn     is the number to set if error                             */
5558
/*   status contains the current status (not yet in context)          */
5559
/*   set    is the context                                            */
5560
/*                                                                    */
5561
/* If the status is an error status, the number is set to a NaN,      */
5562
/* unless the error was an overflow, divide-by-zero, or underflow,    */
5563
/* in which case the number will have already been set.               */
5564
/*                                                                    */
5565
/* The context status is then updated with the new status.  Note that */
5566
/* this may raise a signal, so control may never return from this     */
5567
/* routine (hence resources must be recovered before it is called).   */
5568
/* ------------------------------------------------------------------ */
5569
static void
5570
decStatus (decNumber * dn, uInt status, decContext * set)
5571
{
5572
  if (status & DEC_NaNs)
5573
    {                           /* error status -> NaN */
5574
      /* if cause was an sNaN, clear and propagate [NaN is already set up] */
5575
      if (status & DEC_sNaN)
5576
        status &= ~DEC_sNaN;
5577
      else
5578
        {
5579
          decNumberZero (dn);   /* other error: clean throughout */
5580
          dn->bits = DECNAN;    /* and make a quiet NaN */
5581
        }
5582
    }
5583
  decContextSetStatus (set, status);
5584
  return;
5585
}
5586
 
5587
/* ------------------------------------------------------------------ */
5588
/* decGetDigits -- count digits in a Units array                      */
5589
/*                                                                    */
5590
/*   uar is the Unit array holding the number [this is often an       */
5591
/*          accumulator of some sort]                                 */
5592
/*   len is the length of the array in units                          */
5593
/*                                                                    */
5594
/*   returns the number of (significant) digits in the array          */
5595
/*                                                                    */
5596
/* All leading zeros are excluded, except the last if the array has   */
5597
/* only zero Units.                                                   */
5598
/* ------------------------------------------------------------------ */
5599
/* This may be called twice during some operations. */
5600
static Int
5601
decGetDigits (const Unit * uar, Int len)
5602
{
5603
  const Unit *up = uar + len - 1;       /* -> msu */
5604
  Int digits = len * DECDPUN;   /* maximum possible digits */
5605
  uInt const *pow;              /* work */
5606
 
5607
  for (; up >= uar; up--)
5608
    {
5609
      digits -= DECDPUN;
5610
      if (*up == 0)
5611
        {                       /* unit is 0 */
5612
          if (digits != 0)
5613
            continue;           /* more to check */
5614
          /* all units were 0 */
5615
          digits++;             /* .. so bump digits to 1 */
5616
          break;
5617
        }
5618
      /* found the first non-zero Unit */
5619
      digits++;
5620
      if (*up < 10)
5621
        break;                  /* fastpath 1-9 */
5622
      digits++;
5623
      for (pow = &powers[2]; *up >= *pow; pow++)
5624
        digits++;
5625
      break;
5626
    }                           /* up */
5627
 
5628
  return digits;
5629
}
5630
 
5631
 
5632
#if DECTRACE | DECCHECK
5633
/* ------------------------------------------------------------------ */
5634
/* decNumberShow -- display a number [debug aid]                      */
5635
/*   dn is the number to show                                         */
5636
/*                                                                    */
5637
/* Shows: sign, exponent, coefficient (msu first), digits             */
5638
/*    or: sign, special-value                                         */
5639
/* ------------------------------------------------------------------ */
5640
/* this is public so other modules can use it */
5641
void
5642
decNumberShow (const decNumber * dn)
5643
{
5644
  const Unit *up;               /* work */
5645
  uInt u, d;                    /* .. */
5646
  Int cut;                      /* .. */
5647
  char isign = '+';             /* main sign */
5648
  if (dn == NULL)
5649
    {
5650
      printf ("NULL\n");
5651
      return;
5652
    }
5653
  if (decNumberIsNegative (dn))
5654
    isign = '-';
5655
  printf (" >> %c ", isign);
5656
  if (dn->bits & DECSPECIAL)
5657
    {                           /* Is a special value */
5658
      if (decNumberIsInfinite (dn))
5659
        printf ("Infinity");
5660
      else
5661
        {                       /* a NaN */
5662
          if (dn->bits & DECSNAN)
5663
            printf ("sNaN");    /* signalling NaN */
5664
          else
5665
            printf ("NaN");
5666
        }
5667
      /* if coefficient and exponent are 0, we're done */
5668
      if (dn->exponent == 0 && dn->digits == 1 && *dn->lsu == 0)
5669
        {
5670
          printf ("\n");
5671
          return;
5672
        }
5673
      /* drop through to report other information */
5674
      printf (" ");
5675
    }
5676
 
5677
  /* now carefully display the coefficient */
5678
  up = dn->lsu + D2U (dn->digits) - 1;  /* msu */
5679
  printf ("%d", *up);
5680
  for (up = up - 1; up >= dn->lsu; up--)
5681
    {
5682
      u = *up;
5683
      printf (":");
5684
      for (cut = DECDPUN - 1; cut >= 0; cut--)
5685
        {
5686
          d = u / powers[cut];
5687
          u -= d * powers[cut];
5688
          printf ("%d", d);
5689
        }                       /* cut */
5690
    }                           /* up */
5691
  if (dn->exponent != 0)
5692
    {
5693
      char esign = '+';
5694
      if (dn->exponent < 0)
5695
        esign = '-';
5696
      printf (" E%c%d", esign, abs (dn->exponent));
5697
    }
5698
  printf (" [%d]\n", dn->digits);
5699
}
5700
#endif
5701
 
5702
#if DECTRACE || DECCHECK
5703
/* ------------------------------------------------------------------ */
5704
/* decDumpAr -- display a unit array [debug aid]                      */
5705
/*   name is a single-character tag name                              */
5706
/*   ar   is the array to display                                     */
5707
/*   len  is the length of the array in Units                         */
5708
/* ------------------------------------------------------------------ */
5709
static void
5710
decDumpAr (char name, const Unit * ar, Int len)
5711
{
5712
  Int i;
5713
#if DECDPUN==4
5714
  const char *spec = "%04d ";
5715
#else
5716
  const char *spec = "%d ";
5717
#endif
5718
  printf ("  :%c: ", name);
5719
  for (i = len - 1; i >= 0; i--)
5720
    {
5721
      if (i == len - 1)
5722
        printf ("%d ", ar[i]);
5723
      else
5724
        printf (spec, ar[i]);
5725
    }
5726
  printf ("\n");
5727
  return;
5728
}
5729
#endif
5730
 
5731
#if DECCHECK
5732
/* ------------------------------------------------------------------ */
5733
/* decCheckOperands -- check operand(s) to a routine                  */
5734
/*   res is the result structure (not checked; it will be set to      */
5735
/*          quiet NaN if error found (and it is not NULL))            */
5736
/*   lhs is the first operand (may be DECUNUSED)                      */
5737
/*   rhs is the second (may be DECUNUSED)                             */
5738
/*   set is the context (may be DECUNUSED)                            */
5739
/*   returns 0 if both operands, and the context are clean, or 1      */
5740
/*     otherwise (in which case the context will show an error,       */
5741
/*     unless NULL).  Note that res is not cleaned; caller should     */
5742
/*     handle this so res=NULL case is safe.                          */
5743
/* The caller is expected to abandon immediately if 1 is returned.    */
5744
/* ------------------------------------------------------------------ */
5745
static Flag
5746
decCheckOperands (decNumber * res, const decNumber * lhs,
5747
                  const decNumber * rhs, decContext * set)
5748
{
5749
  Flag bad = 0;
5750
  if (set == NULL)
5751
    {                           /* oops; hopeless */
5752
#if DECTRACE
5753
      printf ("Context is NULL.\n");
5754
#endif
5755
      bad = 1;
5756
      return 1;
5757
    }
5758
  else if (set != DECUNUSED
5759
           && (set->digits < 1 || set->round < 0
5760
               || set->round >= DEC_ROUND_MAX))
5761
    {
5762
      bad = 1;
5763
#if DECTRACE
5764
      printf ("Bad context [digits=%d round=%d].\n", set->digits, set->round);
5765
#endif
5766
    }
5767
  else
5768
    {
5769
      if (res == NULL)
5770
        {
5771
          bad = 1;
5772
#if DECTRACE
5773
          printf ("Bad result [is NULL].\n");
5774
#endif
5775
        }
5776
      if (!bad && lhs != DECUNUSED)
5777
        bad = (decCheckNumber (lhs, set));
5778
      if (!bad && rhs != DECUNUSED)
5779
        bad = (decCheckNumber (rhs, set));
5780
    }
5781
  if (bad)
5782
    {
5783
      if (set != DECUNUSED)
5784
        decContextSetStatus (set, DEC_Invalid_operation);
5785
      if (res != DECUNUSED && res != NULL)
5786
        {
5787
          decNumberZero (res);
5788
          res->bits = DECNAN;   /* qNaN */
5789
        }
5790
    }
5791
  return bad;
5792
}
5793
 
5794
/* ------------------------------------------------------------------ */
5795
/* decCheckNumber -- check a number                                   */
5796
/*   dn is the number to check                                        */
5797
/*   set is the context (may be DECUNUSED)                            */
5798
/*   returns 0 if the number is clean, or 1 otherwise                 */
5799
/*                                                                    */
5800
/* The number is considered valid if it could be a result from some   */
5801
/* operation in some valid context (not necessarily the current one). */
5802
/* ------------------------------------------------------------------ */
5803
Flag
5804
decCheckNumber (const decNumber * dn, decContext * set)
5805
{
5806
  const Unit *up;               /* work */
5807
  uInt maxuint;                 /* .. */
5808
  Int ae, d, digits;            /* .. */
5809
  Int emin, emax;               /* .. */
5810
 
5811
  if (dn == NULL)
5812
    {                           /* hopeless */
5813
#if DECTRACE
5814
      printf ("Reference to decNumber is NULL.\n");
5815
#endif
5816
      return 1;
5817
    }
5818
 
5819
  /* check special values */
5820
  if (dn->bits & DECSPECIAL)
5821
    {
5822
      if (dn->exponent != 0)
5823
        {
5824
#if DECTRACE
5825
          printf ("Exponent %d (not 0) for a special value.\n", dn->exponent);
5826
#endif
5827
          return 1;
5828
        }
5829
 
5830
      /* 2003.09.08: NaNs may now have coefficients, so next tests Inf only */
5831
      if (decNumberIsInfinite (dn))
5832
        {
5833
          if (dn->digits != 1)
5834
            {
5835
#if DECTRACE
5836
              printf ("Digits %d (not 1) for an infinity.\n", dn->digits);
5837
#endif
5838
              return 1;
5839
            }
5840
          if (*dn->lsu != 0)
5841
            {
5842
#if DECTRACE
5843
              printf ("LSU %d (not 0) for an infinity.\n", *dn->lsu);
5844
#endif
5845
              return 1;
5846
            }
5847
        }                       /* Inf */
5848
      /* 2002.12.26: negative NaNs can now appear through proposed IEEE */
5849
      /*             concrete formats (decimal64, etc.), though they are */
5850
      /*             never visible in strings. */
5851
      return 0;
5852
 
5853
      /* if ((dn->bits & DECINF) || (dn->bits & DECNEG)==0) return 0; */
5854
      /* #if DECTRACE */
5855
      /* printf("Negative NaN in number.\n"); */
5856
      /* #endif */
5857
      /* return 1; */
5858
    }
5859
 
5860
  /* check the coefficient */
5861
  if (dn->digits < 1 || dn->digits > DECNUMMAXP)
5862
    {
5863
#if DECTRACE
5864
      printf ("Digits %d in number.\n", dn->digits);
5865
#endif
5866
      return 1;
5867
    }
5868
 
5869
  d = dn->digits;
5870
 
5871
  for (up = dn->lsu; d > 0; up++)
5872
    {
5873
      if (d > DECDPUN)
5874
        maxuint = DECDPUNMAX;
5875
      else
5876
        {                       /* we are at the msu */
5877
          maxuint = powers[d] - 1;
5878
          if (dn->digits > 1 && *up < powers[d - 1])
5879
            {
5880
#if DECTRACE
5881
              printf ("Leading 0 in number.\n");
5882
              decNumberShow (dn);
5883
#endif
5884
              return 1;
5885
            }
5886
        }
5887
      if (*up > maxuint)
5888
        {
5889
#if DECTRACE
5890
          printf ("Bad Unit [%08x] in number at offset %d [maxuint %d].\n",
5891
                  *up, up - dn->lsu, maxuint);
5892
#endif
5893
          return 1;
5894
        }
5895
      d -= DECDPUN;
5896
    }
5897
 
5898
  /* check the exponent.  Note that input operands can have exponents */
5899
  /* which are out of the set->emin/set->emax and set->digits range */
5900
  /* (just as they can have more digits than set->digits). */
5901
  ae = dn->exponent + dn->digits - 1;   /* adjusted exponent */
5902
  emax = DECNUMMAXE;
5903
  emin = DECNUMMINE;
5904
  digits = DECNUMMAXP;
5905
  if (ae < emin - (digits - 1))
5906
    {
5907
#if DECTRACE
5908
      printf ("Adjusted exponent underflow [%d].\n", ae);
5909
      decNumberShow (dn);
5910
#endif
5911
      return 1;
5912
    }
5913
  if (ae > +emax)
5914
    {
5915
#if DECTRACE
5916
      printf ("Adjusted exponent overflow [%d].\n", ae);
5917
      decNumberShow (dn);
5918
#endif
5919
      return 1;
5920
    }
5921
 
5922
  return 0;                      /* it's OK */
5923
}
5924
#endif
5925
 
5926
#if DECALLOC
5927
#undef malloc
5928
#undef free
5929
/* ------------------------------------------------------------------ */
5930
/* decMalloc -- accountable allocation routine                        */
5931
/*   n is the number of bytes to allocate                             */
5932
/*                                                                    */
5933
/* Semantics is the same as the stdlib malloc routine, but bytes      */
5934
/* allocated are accounted for globally, and corruption fences are    */
5935
/* added before and after the 'actual' storage.                       */
5936
/* ------------------------------------------------------------------ */
5937
/* This routine allocates storage with an extra twelve bytes; 8 are   */
5938
/* at the start and hold:                                             */
5939
/*   0-3 the original length requested                                */
5940
/*   4-7 buffer corruption detection fence (DECFENCE, x4)             */
5941
/* The 4 bytes at the end also hold a corruption fence (DECFENCE, x4) */
5942
/* ------------------------------------------------------------------ */
5943
static void *
5944
decMalloc (uInt n)
5945
{
5946
  uInt size = n + 12;           /* true size */
5947
  void *alloc;                  /* -> allocated storage */
5948
  uInt *j;                      /* work */
5949
  uByte *b, *b0;                /* .. */
5950
 
5951
  alloc = malloc (size);        /* -> allocated storage */
5952
  if (alloc == NULL)
5953
    return NULL;                /* out of strorage */
5954
  b0 = (uByte *) alloc;         /* as bytes */
5955
  decAllocBytes += n;           /* account for storage */
5956
  j = (uInt *) alloc;           /* -> first four bytes */
5957
  *j = n;                       /* save n */
5958
  /* printf("++ alloc(%d)\n", n); */
5959
  for (b = b0 + 4; b < b0 + 8; b++)
5960
    *b = DECFENCE;
5961
  for (b = b0 + n + 8; b < b0 + n + 12; b++)
5962
    *b = DECFENCE;
5963
  return b0 + 8;                /* -> play area */
5964
}
5965
 
5966
/* ------------------------------------------------------------------ */
5967
/* decFree -- accountable free routine                                */
5968
/*   alloc is the storage to free                                     */
5969
/*                                                                    */
5970
/* Semantics is the same as the stdlib malloc routine, except that    */
5971
/* the global storage accounting is updated and the fences are        */
5972
/* checked to ensure that no routine has written 'out of bounds'.     */
5973
/* ------------------------------------------------------------------ */
5974
/* This routine first checks that the fences have not been corrupted. */
5975
/* It then frees the storage using the 'truw' storage address (that   */
5976
/* is, offset by 8).                                                  */
5977
/* ------------------------------------------------------------------ */
5978
static void
5979
decFree (void *alloc)
5980
{
5981
  uInt *j, n;                   /* pointer, original length */
5982
  uByte *b, *b0;                /* work */
5983
 
5984
  if (alloc == NULL)
5985
    return;                     /* allowed; it's a nop */
5986
  b0 = (uByte *) alloc;         /* as bytes */
5987
  b0 -= 8;                      /* -> true start of storage */
5988
  j = (uInt *) b0;              /* -> first four bytes */
5989
  n = *j;                       /* lift */
5990
  for (b = b0 + 4; b < b0 + 8; b++)
5991
    if (*b != DECFENCE)
5992
      printf ("=== Corrupt byte [%02x] at offset %d from %d ===\n", *b,
5993
              b - b0 - 8, (Int) b0);
5994
  for (b = b0 + n + 8; b < b0 + n + 12; b++)
5995
    if (*b != DECFENCE)
5996
      printf ("=== Corrupt byte [%02x] at offset +%d from %d, n=%d ===\n", *b,
5997
              b - b0 - 8, (Int) b0, n);
5998
  free (b0);                    /* drop the storage */
5999
  decAllocBytes -= n;           /* account for storage */
6000
}
6001
#endif

powered by: WebSVN 2.1.0

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