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 38

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

powered by: WebSVN 2.1.0

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