qxl: check release info object
[qemu.git] / libdecnumber / decNumber.c
1 /* Decimal number arithmetic module for the decNumber C Library.
2 Copyright (C) 2005, 2007 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 /* Decimal Number arithmetic module */
33 /* ------------------------------------------------------------------ */
34 /* This module comprises the routines for General Decimal Arithmetic */
35 /* as defined in the specification which may be found on the */
36 /* http://www2.hursley.ibm.com/decimal web pages. It implements both */
37 /* the full ('extended') arithmetic and the simpler ('subset') */
38 /* arithmetic. */
39 /* */
40 /* Usage notes: */
41 /* */
42 /* 1. This code is ANSI C89 except: */
43 /* */
44 /* If DECDPUN>4 or DECUSE64=1, the C99 64-bit int64_t and */
45 /* uint64_t types may be used. To avoid these, set DECUSE64=0 */
46 /* and DECDPUN<=4 (see documentation). */
47 /* */
48 /* 2. The decNumber format which this library uses is optimized for */
49 /* efficient processing of relatively short numbers; in particular */
50 /* it allows the use of fixed sized structures and minimizes copy */
51 /* and move operations. It does, however, support arbitrary */
52 /* precision (up to 999,999,999 digits) and arbitrary exponent */
53 /* range (Emax in the range 0 through 999,999,999 and Emin in the */
54 /* range -999,999,999 through 0). Mathematical functions (for */
55 /* example decNumberExp) as identified below are restricted more */
56 /* tightly: digits, emax, and -emin in the context must be <= */
57 /* DEC_MAX_MATH (999999), and their operand(s) must be within */
58 /* these bounds. */
59 /* */
60 /* 3. Logical functions are further restricted; their operands must */
61 /* be finite, positive, have an exponent of zero, and all digits */
62 /* must be either 0 or 1. The result will only contain digits */
63 /* which are 0 or 1 (and will have exponent=0 and a sign of 0). */
64 /* */
65 /* 4. Operands to operator functions are never modified unless they */
66 /* are also specified to be the result number (which is always */
67 /* permitted). Other than that case, operands must not overlap. */
68 /* */
69 /* 5. Error handling: the type of the error is ORed into the status */
70 /* flags in the current context (decContext structure). The */
71 /* SIGFPE signal is then raised if the corresponding trap-enabler */
72 /* flag in the decContext is set (is 1). */
73 /* */
74 /* It is the responsibility of the caller to clear the status */
75 /* flags as required. */
76 /* */
77 /* The result of any routine which returns a number will always */
78 /* be a valid number (which may be a special value, such as an */
79 /* Infinity or NaN). */
80 /* */
81 /* 6. The decNumber format is not an exchangeable concrete */
82 /* representation as it comprises fields which may be machine- */
83 /* dependent (packed or unpacked, or special length, for example). */
84 /* Canonical conversions to and from strings are provided; other */
85 /* conversions are available in separate modules. */
86 /* */
87 /* 7. Normally, input operands are assumed to be valid. Set DECCHECK */
88 /* to 1 for extended operand checking (including NULL operands). */
89 /* Results are undefined if a badly-formed structure (or a NULL */
90 /* pointer to a structure) is provided, though with DECCHECK */
91 /* enabled the operator routines are protected against exceptions. */
92 /* (Except if the result pointer is NULL, which is unrecoverable.) */
93 /* */
94 /* However, the routines will never cause exceptions if they are */
95 /* given well-formed operands, even if the value of the operands */
96 /* is inappropriate for the operation and DECCHECK is not set. */
97 /* (Except for SIGFPE, as and where documented.) */
98 /* */
99 /* 8. Subset arithmetic is available only if DECSUBSET is set to 1. */
100 /* ------------------------------------------------------------------ */
101 /* Implementation notes for maintenance of this module: */
102 /* */
103 /* 1. Storage leak protection: Routines which use malloc are not */
104 /* permitted to use return for fastpath or error exits (i.e., */
105 /* they follow strict structured programming conventions). */
106 /* Instead they have a do{}while(0); construct surrounding the */
107 /* code which is protected -- break may be used to exit this. */
108 /* Other routines can safely use the return statement inline. */
109 /* */
110 /* Storage leak accounting can be enabled using DECALLOC. */
111 /* */
112 /* 2. All loops use the for(;;) construct. Any do construct does */
113 /* not loop; it is for allocation protection as just described. */
114 /* */
115 /* 3. Setting status in the context must always be the very last */
116 /* action in a routine, as non-0 status may raise a trap and hence */
117 /* the call to set status may not return (if the handler uses long */
118 /* jump). Therefore all cleanup must be done first. In general, */
119 /* to achieve this status is accumulated and is only applied just */
120 /* before return by calling decContextSetStatus (via decStatus). */
121 /* */
122 /* Routines which allocate storage cannot, in general, use the */
123 /* 'top level' routines which could cause a non-returning */
124 /* transfer of control. The decXxxxOp routines are safe (do not */
125 /* call decStatus even if traps are set in the context) and should */
126 /* be used instead (they are also a little faster). */
127 /* */
128 /* 4. Exponent checking is minimized by allowing the exponent to */
129 /* grow outside its limits during calculations, provided that */
130 /* the decFinalize function is called later. Multiplication and */
131 /* division, and intermediate calculations in exponentiation, */
132 /* require more careful checks because of the risk of 31-bit */
133 /* overflow (the most negative valid exponent is -1999999997, for */
134 /* a 999999999-digit number with adjusted exponent of -999999999). */
135 /* */
136 /* 5. Rounding is deferred until finalization of results, with any */
137 /* 'off to the right' data being represented as a single digit */
138 /* residue (in the range -1 through 9). This avoids any double- */
139 /* rounding when more than one shortening takes place (for */
140 /* example, when a result is subnormal). */
141 /* */
142 /* 6. The digits count is allowed to rise to a multiple of DECDPUN */
143 /* during many operations, so whole Units are handled and exact */
144 /* accounting of digits is not needed. The correct digits value */
145 /* is found by decGetDigits, which accounts for leading zeros. */
146 /* This must be called before any rounding if the number of digits */
147 /* is not known exactly. */
148 /* */
149 /* 7. The multiply-by-reciprocal 'trick' is used for partitioning */
150 /* numbers up to four digits, using appropriate constants. This */
151 /* is not useful for longer numbers because overflow of 32 bits */
152 /* would lead to 4 multiplies, which is almost as expensive as */
153 /* a divide (unless a floating-point or 64-bit multiply is */
154 /* assumed to be available). */
155 /* */
156 /* 8. Unusual abbreviations that may be used in the commentary: */
157 /* lhs -- left hand side (operand, of an operation) */
158 /* lsd -- least significant digit (of coefficient) */
159 /* lsu -- least significant Unit (of coefficient) */
160 /* msd -- most significant digit (of coefficient) */
161 /* msi -- most significant item (in an array) */
162 /* msu -- most significant Unit (of coefficient) */
163 /* rhs -- right hand side (operand, of an operation) */
164 /* +ve -- positive */
165 /* -ve -- negative */
166 /* ** -- raise to the power */
167 /* ------------------------------------------------------------------ */
168
169 #include "qemu/osdep.h"
170 #include "libdecnumber/dconfig.h"
171 #include "libdecnumber/decNumber.h"
172 #include "libdecnumber/decNumberLocal.h"
173
174 /* Constants */
175 /* Public lookup table used by the D2U macro */
176 const uByte d2utable[DECMAXD2U+1]=D2UTABLE;
177
178 #define DECVERB 1 /* set to 1 for verbose DECCHECK */
179 #define powers DECPOWERS /* old internal name */
180
181 /* Local constants */
182 #define DIVIDE 0x80 /* Divide operators */
183 #define REMAINDER 0x40 /* .. */
184 #define DIVIDEINT 0x20 /* .. */
185 #define REMNEAR 0x10 /* .. */
186 #define COMPARE 0x01 /* Compare operators */
187 #define COMPMAX 0x02 /* .. */
188 #define COMPMIN 0x03 /* .. */
189 #define COMPTOTAL 0x04 /* .. */
190 #define COMPNAN 0x05 /* .. [NaN processing] */
191 #define COMPSIG 0x06 /* .. [signaling COMPARE] */
192 #define COMPMAXMAG 0x07 /* .. */
193 #define COMPMINMAG 0x08 /* .. */
194
195 #define DEC_sNaN 0x40000000 /* local status: sNaN signal */
196 #define BADINT (Int)0x80000000 /* most-negative Int; error indicator */
197 /* Next two indicate an integer >= 10**6, and its parity (bottom bit) */
198 #define BIGEVEN (Int)0x80000002
199 #define BIGODD (Int)0x80000003
200
201 static Unit uarrone[1]={1}; /* Unit array of 1, used for incrementing */
202
203 /* Granularity-dependent code */
204 #if DECDPUN<=4
205 #define eInt Int /* extended integer */
206 #define ueInt uInt /* unsigned extended integer */
207 /* Constant multipliers for divide-by-power-of five using reciprocal */
208 /* multiply, after removing powers of 2 by shifting, and final shift */
209 /* of 17 [we only need up to **4] */
210 static const uInt multies[]={131073, 26215, 5243, 1049, 210};
211 /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
212 #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
213 #else
214 /* For DECDPUN>4 non-ANSI-89 64-bit types are needed. */
215 #if !DECUSE64
216 #error decNumber.c: DECUSE64 must be 1 when DECDPUN>4
217 #endif
218 #define eInt Long /* extended integer */
219 #define ueInt uLong /* unsigned extended integer */
220 #endif
221
222 /* Local routines */
223 static decNumber * decAddOp(decNumber *, const decNumber *, const decNumber *,
224 decContext *, uByte, uInt *);
225 static Flag decBiStr(const char *, const char *, const char *);
226 static uInt decCheckMath(const decNumber *, decContext *, uInt *);
227 static void decApplyRound(decNumber *, decContext *, Int, uInt *);
228 static Int decCompare(const decNumber *lhs, const decNumber *rhs, Flag);
229 static decNumber * decCompareOp(decNumber *, const decNumber *,
230 const decNumber *, decContext *,
231 Flag, uInt *);
232 static void decCopyFit(decNumber *, const decNumber *, decContext *,
233 Int *, uInt *);
234 static decNumber * decDecap(decNumber *, Int);
235 static decNumber * decDivideOp(decNumber *, const decNumber *,
236 const decNumber *, decContext *, Flag, uInt *);
237 static decNumber * decExpOp(decNumber *, const decNumber *,
238 decContext *, uInt *);
239 static void decFinalize(decNumber *, decContext *, Int *, uInt *);
240 static Int decGetDigits(Unit *, Int);
241 static Int decGetInt(const decNumber *);
242 static decNumber * decLnOp(decNumber *, const decNumber *,
243 decContext *, uInt *);
244 static decNumber * decMultiplyOp(decNumber *, const decNumber *,
245 const decNumber *, decContext *,
246 uInt *);
247 static decNumber * decNaNs(decNumber *, const decNumber *,
248 const decNumber *, decContext *, uInt *);
249 static decNumber * decQuantizeOp(decNumber *, const decNumber *,
250 const decNumber *, decContext *, Flag,
251 uInt *);
252 static void decReverse(Unit *, Unit *);
253 static void decSetCoeff(decNumber *, decContext *, const Unit *,
254 Int, Int *, uInt *);
255 static void decSetMaxValue(decNumber *, decContext *);
256 static void decSetOverflow(decNumber *, decContext *, uInt *);
257 static void decSetSubnormal(decNumber *, decContext *, Int *, uInt *);
258 static Int decShiftToLeast(Unit *, Int, Int);
259 static Int decShiftToMost(Unit *, Int, Int);
260 static void decStatus(decNumber *, uInt, decContext *);
261 static void decToString(const decNumber *, char[], Flag);
262 static decNumber * decTrim(decNumber *, decContext *, Flag, Int *);
263 static Int decUnitAddSub(const Unit *, Int, const Unit *, Int, Int,
264 Unit *, Int);
265 static Int decUnitCompare(const Unit *, Int, const Unit *, Int, Int);
266
267 #if !DECSUBSET
268 /* decFinish == decFinalize when no subset arithmetic needed */
269 #define decFinish(a,b,c,d) decFinalize(a,b,c,d)
270 #else
271 static void decFinish(decNumber *, decContext *, Int *, uInt *);
272 static decNumber * decRoundOperand(const decNumber *, decContext *, uInt *);
273 #endif
274
275 /* Local macros */
276 /* masked special-values bits */
277 #define SPECIALARG (rhs->bits & DECSPECIAL)
278 #define SPECIALARGS ((lhs->bits | rhs->bits) & DECSPECIAL)
279
280 /* Diagnostic macros, etc. */
281 #if DECALLOC
282 /* Handle malloc/free accounting. If enabled, our accountable routines */
283 /* are used; otherwise the code just goes straight to the system malloc */
284 /* and free routines. */
285 #define malloc(a) decMalloc(a)
286 #define free(a) decFree(a)
287 #define DECFENCE 0x5a /* corruption detector */
288 /* 'Our' malloc and free: */
289 static void *decMalloc(size_t);
290 static void decFree(void *);
291 uInt decAllocBytes=0; /* count of bytes allocated */
292 /* Note that DECALLOC code only checks for storage buffer overflow. */
293 /* To check for memory leaks, the decAllocBytes variable must be */
294 /* checked to be 0 at appropriate times (e.g., after the test */
295 /* harness completes a set of tests). This checking may be unreliable */
296 /* if the testing is done in a multi-thread environment. */
297 #endif
298
299 #if DECCHECK
300 /* Optional checking routines. Enabling these means that decNumber */
301 /* and decContext operands to operator routines are checked for */
302 /* correctness. This roughly doubles the execution time of the */
303 /* fastest routines (and adds 600+ bytes), so should not normally be */
304 /* used in 'production'. */
305 /* decCheckInexact is used to check that inexact results have a full */
306 /* complement of digits (where appropriate -- this is not the case */
307 /* for Quantize, for example) */
308 #define DECUNRESU ((decNumber *)(void *)0xffffffff)
309 #define DECUNUSED ((const decNumber *)(void *)0xffffffff)
310 #define DECUNCONT ((decContext *)(void *)(0xffffffff))
311 static Flag decCheckOperands(decNumber *, const decNumber *,
312 const decNumber *, decContext *);
313 static Flag decCheckNumber(const decNumber *);
314 static void decCheckInexact(const decNumber *, decContext *);
315 #endif
316
317 #if DECTRACE || DECCHECK
318 /* Optional trace/debugging routines (may or may not be used) */
319 void decNumberShow(const decNumber *); /* displays the components of a number */
320 static void decDumpAr(char, const Unit *, Int);
321 #endif
322
323 /* ================================================================== */
324 /* Conversions */
325 /* ================================================================== */
326
327 /* ------------------------------------------------------------------ */
328 /* from-int32 -- conversion from Int or uInt */
329 /* */
330 /* dn is the decNumber to receive the integer */
331 /* in or uin is the integer to be converted */
332 /* returns dn */
333 /* */
334 /* No error is possible. */
335 /* ------------------------------------------------------------------ */
336 decNumber * decNumberFromInt32(decNumber *dn, Int in) {
337 uInt unsig;
338 if (in>=0) unsig=in;
339 else { /* negative (possibly BADINT) */
340 if (in==BADINT) unsig=(uInt)1073741824*2; /* special case */
341 else unsig=-in; /* invert */
342 }
343 /* in is now positive */
344 decNumberFromUInt32(dn, unsig);
345 if (in<0) dn->bits=DECNEG; /* sign needed */
346 return dn;
347 } /* decNumberFromInt32 */
348
349 decNumber * decNumberFromUInt32(decNumber *dn, uInt uin) {
350 Unit *up; /* work pointer */
351 decNumberZero(dn); /* clean */
352 if (uin==0) return dn; /* [or decGetDigits bad call] */
353 for (up=dn->lsu; uin>0; up++) {
354 *up=(Unit)(uin%(DECDPUNMAX+1));
355 uin=uin/(DECDPUNMAX+1);
356 }
357 dn->digits=decGetDigits(dn->lsu, up-dn->lsu);
358 return dn;
359 } /* decNumberFromUInt32 */
360
361 /* ------------------------------------------------------------------ */
362 /* to-int32 -- conversion to Int or uInt */
363 /* */
364 /* dn is the decNumber to convert */
365 /* set is the context for reporting errors */
366 /* returns the converted decNumber, or 0 if Invalid is set */
367 /* */
368 /* Invalid is set if the decNumber does not have exponent==0 or if */
369 /* it is a NaN, Infinite, or out-of-range. */
370 /* ------------------------------------------------------------------ */
371 Int decNumberToInt32(const decNumber *dn, decContext *set) {
372 #if DECCHECK
373 if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
374 #endif
375
376 /* special or too many digits, or bad exponent */
377 if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0) ; /* bad */
378 else { /* is a finite integer with 10 or fewer digits */
379 Int d; /* work */
380 const Unit *up; /* .. */
381 uInt hi=0, lo; /* .. */
382 up=dn->lsu; /* -> lsu */
383 lo=*up; /* get 1 to 9 digits */
384 #if DECDPUN>1 /* split to higher */
385 hi=lo/10;
386 lo=lo%10;
387 #endif
388 up++;
389 /* collect remaining Units, if any, into hi */
390 for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1];
391 /* now low has the lsd, hi the remainder */
392 if (hi>214748364 || (hi==214748364 && lo>7)) { /* out of range? */
393 /* most-negative is a reprieve */
394 if (dn->bits&DECNEG && hi==214748364 && lo==8) return 0x80000000;
395 /* bad -- drop through */
396 }
397 else { /* in-range always */
398 Int i=X10(hi)+lo;
399 if (dn->bits&DECNEG) return -i;
400 return i;
401 }
402 } /* integer */
403 decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */
404 return 0;
405 } /* decNumberToInt32 */
406
407 uInt decNumberToUInt32(const decNumber *dn, decContext *set) {
408 #if DECCHECK
409 if (decCheckOperands(DECUNRESU, DECUNUSED, dn, set)) return 0;
410 #endif
411 /* special or too many digits, or bad exponent, or negative (<0) */
412 if (dn->bits&DECSPECIAL || dn->digits>10 || dn->exponent!=0
413 || (dn->bits&DECNEG && !ISZERO(dn))); /* bad */
414 else { /* is a finite integer with 10 or fewer digits */
415 Int d; /* work */
416 const Unit *up; /* .. */
417 uInt hi=0, lo; /* .. */
418 up=dn->lsu; /* -> lsu */
419 lo=*up; /* get 1 to 9 digits */
420 #if DECDPUN>1 /* split to higher */
421 hi=lo/10;
422 lo=lo%10;
423 #endif
424 up++;
425 /* collect remaining Units, if any, into hi */
426 for (d=DECDPUN; d<dn->digits; up++, d+=DECDPUN) hi+=*up*powers[d-1];
427
428 /* now low has the lsd, hi the remainder */
429 if (hi>429496729 || (hi==429496729 && lo>5)) ; /* no reprieve possible */
430 else return X10(hi)+lo;
431 } /* integer */
432 decContextSetStatus(set, DEC_Invalid_operation); /* [may not return] */
433 return 0;
434 } /* decNumberToUInt32 */
435
436 decNumber *decNumberFromInt64(decNumber *dn, int64_t in)
437 {
438 uint64_t unsig = in;
439 if (in < 0) {
440 unsig = -unsig;
441 }
442
443 decNumberFromUInt64(dn, unsig);
444 if (in < 0) {
445 dn->bits = DECNEG; /* sign needed */
446 }
447 return dn;
448 } /* decNumberFromInt64 */
449
450 decNumber *decNumberFromUInt64(decNumber *dn, uint64_t uin)
451 {
452 Unit *up; /* work pointer */
453 decNumberZero(dn); /* clean */
454 if (uin == 0) {
455 return dn; /* [or decGetDigits bad call] */
456 }
457 for (up = dn->lsu; uin > 0; up++) {
458 *up = (Unit)(uin % (DECDPUNMAX + 1));
459 uin = uin / (DECDPUNMAX + 1);
460 }
461 dn->digits = decGetDigits(dn->lsu, up-dn->lsu);
462 return dn;
463 } /* decNumberFromUInt64 */
464
465 /* ------------------------------------------------------------------ */
466 /* to-int64 -- conversion to int64 */
467 /* */
468 /* dn is the decNumber to convert. dn is assumed to have been */
469 /* rounded to a floating point integer value. */
470 /* set is the context for reporting errors */
471 /* returns the converted decNumber, or 0 if Invalid is set */
472 /* */
473 /* Invalid is set if the decNumber is a NaN, Infinite or is out of */
474 /* range for a signed 64 bit integer. */
475 /* ------------------------------------------------------------------ */
476
477 int64_t decNumberIntegralToInt64(const decNumber *dn, decContext *set)
478 {
479 if (decNumberIsSpecial(dn) || (dn->exponent < 0) ||
480 (dn->digits + dn->exponent > 19)) {
481 goto Invalid;
482 } else {
483 int64_t d; /* work */
484 const Unit *up; /* .. */
485 uint64_t hi = 0;
486 up = dn->lsu; /* -> lsu */
487
488 for (d = 1; d <= dn->digits; up++, d += DECDPUN) {
489 uint64_t prev = hi;
490 hi += *up * powers[d-1];
491 if ((hi < prev) || (hi > INT64_MAX)) {
492 goto Invalid;
493 }
494 }
495
496 uint64_t prev = hi;
497 hi *= (uint64_t)powers[dn->exponent];
498 if ((hi < prev) || (hi > INT64_MAX)) {
499 goto Invalid;
500 }
501 return (decNumberIsNegative(dn)) ? -((int64_t)hi) : (int64_t)hi;
502 }
503
504 Invalid:
505 decContextSetStatus(set, DEC_Invalid_operation);
506 return 0;
507 } /* decNumberIntegralToInt64 */
508
509
510 /* ------------------------------------------------------------------ */
511 /* to-scientific-string -- conversion to numeric string */
512 /* to-engineering-string -- conversion to numeric string */
513 /* */
514 /* decNumberToString(dn, string); */
515 /* decNumberToEngString(dn, string); */
516 /* */
517 /* dn is the decNumber to convert */
518 /* string is the string where the result will be laid out */
519 /* */
520 /* string must be at least dn->digits+14 characters long */
521 /* */
522 /* No error is possible, and no status can be set. */
523 /* ------------------------------------------------------------------ */
524 char * decNumberToString(const decNumber *dn, char *string){
525 decToString(dn, string, 0);
526 return string;
527 } /* DecNumberToString */
528
529 char * decNumberToEngString(const decNumber *dn, char *string){
530 decToString(dn, string, 1);
531 return string;
532 } /* DecNumberToEngString */
533
534 /* ------------------------------------------------------------------ */
535 /* to-number -- conversion from numeric string */
536 /* */
537 /* decNumberFromString -- convert string to decNumber */
538 /* dn -- the number structure to fill */
539 /* chars[] -- the string to convert ('\0' terminated) */
540 /* set -- the context used for processing any error, */
541 /* determining the maximum precision available */
542 /* (set.digits), determining the maximum and minimum */
543 /* exponent (set.emax and set.emin), determining if */
544 /* extended values are allowed, and checking the */
545 /* rounding mode if overflow occurs or rounding is */
546 /* needed. */
547 /* */
548 /* The length of the coefficient and the size of the exponent are */
549 /* checked by this routine, so the correct error (Underflow or */
550 /* Overflow) can be reported or rounding applied, as necessary. */
551 /* */
552 /* If bad syntax is detected, the result will be a quiet NaN. */
553 /* ------------------------------------------------------------------ */
554 decNumber * decNumberFromString(decNumber *dn, const char chars[],
555 decContext *set) {
556 Int exponent=0; /* working exponent [assume 0] */
557 uByte bits=0; /* working flags [assume +ve] */
558 Unit *res; /* where result will be built */
559 Unit resbuff[SD2U(DECBUFFER+9)];/* local buffer in case need temporary */
560 /* [+9 allows for ln() constants] */
561 Unit *allocres=NULL; /* -> allocated result, iff allocated */
562 Int d=0; /* count of digits found in decimal part */
563 const char *dotchar=NULL; /* where dot was found */
564 const char *cfirst=chars; /* -> first character of decimal part */
565 const char *last=NULL; /* -> last digit of decimal part */
566 const char *c; /* work */
567 Unit *up; /* .. */
568 #if DECDPUN>1
569 Int cut, out; /* .. */
570 #endif
571 Int residue; /* rounding residue */
572 uInt status=0; /* error code */
573
574 #if DECCHECK
575 if (decCheckOperands(DECUNRESU, DECUNUSED, DECUNUSED, set))
576 return decNumberZero(dn);
577 #endif
578
579 do { /* status & malloc protection */
580 for (c=chars;; c++) { /* -> input character */
581 if (*c>='0' && *c<='9') { /* test for Arabic digit */
582 last=c;
583 d++; /* count of real digits */
584 continue; /* still in decimal part */
585 }
586 if (*c=='.' && dotchar==NULL) { /* first '.' */
587 dotchar=c; /* record offset into decimal part */
588 if (c==cfirst) cfirst++; /* first digit must follow */
589 continue;}
590 if (c==chars) { /* first in string... */
591 if (*c=='-') { /* valid - sign */
592 cfirst++;
593 bits=DECNEG;
594 continue;}
595 if (*c=='+') { /* valid + sign */
596 cfirst++;
597 continue;}
598 }
599 /* *c is not a digit, or a valid +, -, or '.' */
600 break;
601 } /* c */
602
603 if (last==NULL) { /* no digits yet */
604 status=DEC_Conversion_syntax;/* assume the worst */
605 if (*c=='\0') break; /* and no more to come... */
606 #if DECSUBSET
607 /* if subset then infinities and NaNs are not allowed */
608 if (!set->extended) break; /* hopeless */
609 #endif
610 /* Infinities and NaNs are possible, here */
611 if (dotchar!=NULL) break; /* .. unless had a dot */
612 decNumberZero(dn); /* be optimistic */
613 if (decBiStr(c, "infinity", "INFINITY")
614 || decBiStr(c, "inf", "INF")) {
615 dn->bits=bits | DECINF;
616 status=0; /* is OK */
617 break; /* all done */
618 }
619 /* a NaN expected */
620 /* 2003.09.10 NaNs are now permitted to have a sign */
621 dn->bits=bits | DECNAN; /* assume simple NaN */
622 if (*c=='s' || *c=='S') { /* looks like an sNaN */
623 c++;
624 dn->bits=bits | DECSNAN;
625 }
626 if (*c!='n' && *c!='N') break; /* check caseless "NaN" */
627 c++;
628 if (*c!='a' && *c!='A') break; /* .. */
629 c++;
630 if (*c!='n' && *c!='N') break; /* .. */
631 c++;
632 /* now either nothing, or nnnn payload, expected */
633 /* -> start of integer and skip leading 0s [including plain 0] */
634 for (cfirst=c; *cfirst=='0';) cfirst++;
635 if (*cfirst=='\0') { /* "NaN" or "sNaN", maybe with all 0s */
636 status=0; /* it's good */
637 break; /* .. */
638 }
639 /* something other than 0s; setup last and d as usual [no dots] */
640 for (c=cfirst;; c++, d++) {
641 if (*c<'0' || *c>'9') break; /* test for Arabic digit */
642 last=c;
643 }
644 if (*c!='\0') break; /* not all digits */
645 if (d>set->digits-1) {
646 /* [NB: payload in a decNumber can be full length unless */
647 /* clamped, in which case can only be digits-1] */
648 if (set->clamp) break;
649 if (d>set->digits) break;
650 } /* too many digits? */
651 /* good; drop through to convert the integer to coefficient */
652 status=0; /* syntax is OK */
653 bits=dn->bits; /* for copy-back */
654 } /* last==NULL */
655
656 else if (*c!='\0') { /* more to process... */
657 /* had some digits; exponent is only valid sequence now */
658 Flag nege; /* 1=negative exponent */
659 const char *firstexp; /* -> first significant exponent digit */
660 status=DEC_Conversion_syntax;/* assume the worst */
661 if (*c!='e' && *c!='E') break;
662 /* Found 'e' or 'E' -- now process explicit exponent */
663 /* 1998.07.11: sign no longer required */
664 nege=0;
665 c++; /* to (possible) sign */
666 if (*c=='-') {nege=1; c++;}
667 else if (*c=='+') c++;
668 if (*c=='\0') break;
669
670 for (; *c=='0' && *(c+1)!='\0';) c++; /* strip insignificant zeros */
671 firstexp=c; /* save exponent digit place */
672 for (; ;c++) {
673 if (*c<'0' || *c>'9') break; /* not a digit */
674 exponent=X10(exponent)+(Int)*c-(Int)'0';
675 } /* c */
676 /* if not now on a '\0', *c must not be a digit */
677 if (*c!='\0') break;
678
679 /* (this next test must be after the syntax checks) */
680 /* if it was too long the exponent may have wrapped, so check */
681 /* carefully and set it to a certain overflow if wrap possible */
682 if (c>=firstexp+9+1) {
683 if (c>firstexp+9+1 || *firstexp>'1') exponent=DECNUMMAXE*2;
684 /* [up to 1999999999 is OK, for example 1E-1000000998] */
685 }
686 if (nege) exponent=-exponent; /* was negative */
687 status=0; /* is OK */
688 } /* stuff after digits */
689
690 /* Here when whole string has been inspected; syntax is good */
691 /* cfirst->first digit (never dot), last->last digit (ditto) */
692
693 /* strip leading zeros/dot [leave final 0 if all 0's] */
694 if (*cfirst=='0') { /* [cfirst has stepped over .] */
695 for (c=cfirst; c<last; c++, cfirst++) {
696 if (*c=='.') continue; /* ignore dots */
697 if (*c!='0') break; /* non-zero found */
698 d--; /* 0 stripped */
699 } /* c */
700 #if DECSUBSET
701 /* make a rapid exit for easy zeros if !extended */
702 if (*cfirst=='0' && !set->extended) {
703 decNumberZero(dn); /* clean result */
704 break; /* [could be return] */
705 }
706 #endif
707 } /* at least one leading 0 */
708
709 /* Handle decimal point... */
710 if (dotchar!=NULL && dotchar<last) /* non-trailing '.' found? */
711 exponent-=(last-dotchar); /* adjust exponent */
712 /* [we can now ignore the .] */
713
714 /* OK, the digits string is good. Assemble in the decNumber, or in */
715 /* a temporary units array if rounding is needed */
716 if (d<=set->digits) res=dn->lsu; /* fits into supplied decNumber */
717 else { /* rounding needed */
718 Int needbytes=D2U(d)*sizeof(Unit);/* bytes needed */
719 res=resbuff; /* assume use local buffer */
720 if (needbytes>(Int)sizeof(resbuff)) { /* too big for local */
721 allocres=(Unit *)malloc(needbytes);
722 if (allocres==NULL) {status|=DEC_Insufficient_storage; break;}
723 res=allocres;
724 }
725 }
726 /* res now -> number lsu, buffer, or allocated storage for Unit array */
727
728 /* Place the coefficient into the selected Unit array */
729 /* [this is often 70% of the cost of this function when DECDPUN>1] */
730 #if DECDPUN>1
731 out=0; /* accumulator */
732 up=res+D2U(d)-1; /* -> msu */
733 cut=d-(up-res)*DECDPUN; /* digits in top unit */
734 for (c=cfirst;; c++) { /* along the digits */
735 if (*c=='.') continue; /* ignore '.' [don't decrement cut] */
736 out=X10(out)+(Int)*c-(Int)'0';
737 if (c==last) break; /* done [never get to trailing '.'] */
738 cut--;
739 if (cut>0) continue; /* more for this unit */
740 *up=(Unit)out; /* write unit */
741 up--; /* prepare for unit below.. */
742 cut=DECDPUN; /* .. */
743 out=0; /* .. */
744 } /* c */
745 *up=(Unit)out; /* write lsu */
746
747 #else
748 /* DECDPUN==1 */
749 up=res; /* -> lsu */
750 for (c=last; c>=cfirst; c--) { /* over each character, from least */
751 if (*c=='.') continue; /* ignore . [don't step up] */
752 *up=(Unit)((Int)*c-(Int)'0');
753 up++;
754 } /* c */
755 #endif
756
757 dn->bits=bits;
758 dn->exponent=exponent;
759 dn->digits=d;
760
761 /* if not in number (too long) shorten into the number */
762 if (d>set->digits) {
763 residue=0;
764 decSetCoeff(dn, set, res, d, &residue, &status);
765 /* always check for overflow or subnormal and round as needed */
766 decFinalize(dn, set, &residue, &status);
767 }
768 else { /* no rounding, but may still have overflow or subnormal */
769 /* [these tests are just for performance; finalize repeats them] */
770 if ((dn->exponent-1<set->emin-dn->digits)
771 || (dn->exponent-1>set->emax-set->digits)) {
772 residue=0;
773 decFinalize(dn, set, &residue, &status);
774 }
775 }
776 /* decNumberShow(dn); */
777 } while(0); /* [for break] */
778
779 if (allocres!=NULL) free(allocres); /* drop any storage used */
780 if (status!=0) decStatus(dn, status, set);
781 return dn;
782 } /* decNumberFromString */
783
784 /* ================================================================== */
785 /* Operators */
786 /* ================================================================== */
787
788 /* ------------------------------------------------------------------ */
789 /* decNumberAbs -- absolute value operator */
790 /* */
791 /* This computes C = abs(A) */
792 /* */
793 /* res is C, the result. C may be A */
794 /* rhs is A */
795 /* set is the context */
796 /* */
797 /* See also decNumberCopyAbs for a quiet bitwise version of this. */
798 /* C must have space for set->digits digits. */
799 /* ------------------------------------------------------------------ */
800 /* This has the same effect as decNumberPlus unless A is negative, */
801 /* in which case it has the same effect as decNumberMinus. */
802 /* ------------------------------------------------------------------ */
803 decNumber * decNumberAbs(decNumber *res, const decNumber *rhs,
804 decContext *set) {
805 decNumber dzero; /* for 0 */
806 uInt status=0; /* accumulator */
807
808 #if DECCHECK
809 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
810 #endif
811
812 decNumberZero(&dzero); /* set 0 */
813 dzero.exponent=rhs->exponent; /* [no coefficient expansion] */
814 decAddOp(res, &dzero, rhs, set, (uByte)(rhs->bits & DECNEG), &status);
815 if (status!=0) decStatus(res, status, set);
816 #if DECCHECK
817 decCheckInexact(res, set);
818 #endif
819 return res;
820 } /* decNumberAbs */
821
822 /* ------------------------------------------------------------------ */
823 /* decNumberAdd -- add two Numbers */
824 /* */
825 /* This computes C = A + B */
826 /* */
827 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
828 /* lhs is A */
829 /* rhs is B */
830 /* set is the context */
831 /* */
832 /* C must have space for set->digits digits. */
833 /* ------------------------------------------------------------------ */
834 /* This just calls the routine shared with Subtract */
835 decNumber * decNumberAdd(decNumber *res, const decNumber *lhs,
836 const decNumber *rhs, decContext *set) {
837 uInt status=0; /* accumulator */
838 decAddOp(res, lhs, rhs, set, 0, &status);
839 if (status!=0) decStatus(res, status, set);
840 #if DECCHECK
841 decCheckInexact(res, set);
842 #endif
843 return res;
844 } /* decNumberAdd */
845
846 /* ------------------------------------------------------------------ */
847 /* decNumberAnd -- AND two Numbers, digitwise */
848 /* */
849 /* This computes C = A & B */
850 /* */
851 /* res is C, the result. C may be A and/or B (e.g., X=X&X) */
852 /* lhs is A */
853 /* rhs is B */
854 /* set is the context (used for result length and error report) */
855 /* */
856 /* C must have space for set->digits digits. */
857 /* */
858 /* Logical function restrictions apply (see above); a NaN is */
859 /* returned with Invalid_operation if a restriction is violated. */
860 /* ------------------------------------------------------------------ */
861 decNumber * decNumberAnd(decNumber *res, const decNumber *lhs,
862 const decNumber *rhs, decContext *set) {
863 const Unit *ua, *ub; /* -> operands */
864 const Unit *msua, *msub; /* -> operand msus */
865 Unit *uc, *msuc; /* -> result and its msu */
866 Int msudigs; /* digits in res msu */
867 #if DECCHECK
868 if (decCheckOperands(res, lhs, rhs, set)) return res;
869 #endif
870
871 if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs)
872 || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
873 decStatus(res, DEC_Invalid_operation, set);
874 return res;
875 }
876
877 /* operands are valid */
878 ua=lhs->lsu; /* bottom-up */
879 ub=rhs->lsu; /* .. */
880 uc=res->lsu; /* .. */
881 msua=ua+D2U(lhs->digits)-1; /* -> msu of lhs */
882 msub=ub+D2U(rhs->digits)-1; /* -> msu of rhs */
883 msuc=uc+D2U(set->digits)-1; /* -> msu of result */
884 msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */
885 for (; uc<=msuc; ua++, ub++, uc++) { /* Unit loop */
886 Unit a, b; /* extract units */
887 if (ua>msua) a=0;
888 else a=*ua;
889 if (ub>msub) b=0;
890 else b=*ub;
891 *uc=0; /* can now write back */
892 if (a|b) { /* maybe 1 bits to examine */
893 Int i, j;
894 *uc=0; /* can now write back */
895 /* This loop could be unrolled and/or use BIN2BCD tables */
896 for (i=0; i<DECDPUN; i++) {
897 if (a&b&1) *uc=*uc+(Unit)powers[i]; /* effect AND */
898 j=a%10;
899 a=a/10;
900 j|=b%10;
901 b=b/10;
902 if (j>1) {
903 decStatus(res, DEC_Invalid_operation, set);
904 return res;
905 }
906 if (uc==msuc && i==msudigs-1) break; /* just did final digit */
907 } /* each digit */
908 } /* both OK */
909 } /* each unit */
910 /* [here uc-1 is the msu of the result] */
911 res->digits=decGetDigits(res->lsu, uc-res->lsu);
912 res->exponent=0; /* integer */
913 res->bits=0; /* sign=0 */
914 return res; /* [no status to set] */
915 } /* decNumberAnd */
916
917 /* ------------------------------------------------------------------ */
918 /* decNumberCompare -- compare two Numbers */
919 /* */
920 /* This computes C = A ? B */
921 /* */
922 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
923 /* lhs is A */
924 /* rhs is B */
925 /* set is the context */
926 /* */
927 /* C must have space for one digit (or NaN). */
928 /* ------------------------------------------------------------------ */
929 decNumber * decNumberCompare(decNumber *res, const decNumber *lhs,
930 const decNumber *rhs, decContext *set) {
931 uInt status=0; /* accumulator */
932 decCompareOp(res, lhs, rhs, set, COMPARE, &status);
933 if (status!=0) decStatus(res, status, set);
934 return res;
935 } /* decNumberCompare */
936
937 /* ------------------------------------------------------------------ */
938 /* decNumberCompareSignal -- compare, signalling on all NaNs */
939 /* */
940 /* This computes C = A ? B */
941 /* */
942 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
943 /* lhs is A */
944 /* rhs is B */
945 /* set is the context */
946 /* */
947 /* C must have space for one digit (or NaN). */
948 /* ------------------------------------------------------------------ */
949 decNumber * decNumberCompareSignal(decNumber *res, const decNumber *lhs,
950 const decNumber *rhs, decContext *set) {
951 uInt status=0; /* accumulator */
952 decCompareOp(res, lhs, rhs, set, COMPSIG, &status);
953 if (status!=0) decStatus(res, status, set);
954 return res;
955 } /* decNumberCompareSignal */
956
957 /* ------------------------------------------------------------------ */
958 /* decNumberCompareTotal -- compare two Numbers, using total ordering */
959 /* */
960 /* This computes C = A ? B, under total ordering */
961 /* */
962 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
963 /* lhs is A */
964 /* rhs is B */
965 /* set is the context */
966 /* */
967 /* C must have space for one digit; the result will always be one of */
968 /* -1, 0, or 1. */
969 /* ------------------------------------------------------------------ */
970 decNumber * decNumberCompareTotal(decNumber *res, const decNumber *lhs,
971 const decNumber *rhs, decContext *set) {
972 uInt status=0; /* accumulator */
973 decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status);
974 if (status!=0) decStatus(res, status, set);
975 return res;
976 } /* decNumberCompareTotal */
977
978 /* ------------------------------------------------------------------ */
979 /* decNumberCompareTotalMag -- compare, total ordering of magnitudes */
980 /* */
981 /* This computes C = |A| ? |B|, under total ordering */
982 /* */
983 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
984 /* lhs is A */
985 /* rhs is B */
986 /* set is the context */
987 /* */
988 /* C must have space for one digit; the result will always be one of */
989 /* -1, 0, or 1. */
990 /* ------------------------------------------------------------------ */
991 decNumber * decNumberCompareTotalMag(decNumber *res, const decNumber *lhs,
992 const decNumber *rhs, decContext *set) {
993 uInt status=0; /* accumulator */
994 uInt needbytes; /* for space calculations */
995 decNumber bufa[D2N(DECBUFFER+1)];/* +1 in case DECBUFFER=0 */
996 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */
997 decNumber bufb[D2N(DECBUFFER+1)];
998 decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */
999 decNumber *a, *b; /* temporary pointers */
1000
1001 #if DECCHECK
1002 if (decCheckOperands(res, lhs, rhs, set)) return res;
1003 #endif
1004
1005 do { /* protect allocated storage */
1006 /* if either is negative, take a copy and absolute */
1007 if (decNumberIsNegative(lhs)) { /* lhs<0 */
1008 a=bufa;
1009 needbytes=sizeof(decNumber)+(D2U(lhs->digits)-1)*sizeof(Unit);
1010 if (needbytes>sizeof(bufa)) { /* need malloc space */
1011 allocbufa=(decNumber *)malloc(needbytes);
1012 if (allocbufa==NULL) { /* hopeless -- abandon */
1013 status|=DEC_Insufficient_storage;
1014 break;}
1015 a=allocbufa; /* use the allocated space */
1016 }
1017 decNumberCopy(a, lhs); /* copy content */
1018 a->bits&=~DECNEG; /* .. and clear the sign */
1019 lhs=a; /* use copy from here on */
1020 }
1021 if (decNumberIsNegative(rhs)) { /* rhs<0 */
1022 b=bufb;
1023 needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
1024 if (needbytes>sizeof(bufb)) { /* need malloc space */
1025 allocbufb=(decNumber *)malloc(needbytes);
1026 if (allocbufb==NULL) { /* hopeless -- abandon */
1027 status|=DEC_Insufficient_storage;
1028 break;}
1029 b=allocbufb; /* use the allocated space */
1030 }
1031 decNumberCopy(b, rhs); /* copy content */
1032 b->bits&=~DECNEG; /* .. and clear the sign */
1033 rhs=b; /* use copy from here on */
1034 }
1035 decCompareOp(res, lhs, rhs, set, COMPTOTAL, &status);
1036 } while(0); /* end protected */
1037
1038 if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
1039 if (allocbufb!=NULL) free(allocbufb); /* .. */
1040 if (status!=0) decStatus(res, status, set);
1041 return res;
1042 } /* decNumberCompareTotalMag */
1043
1044 /* ------------------------------------------------------------------ */
1045 /* decNumberDivide -- divide one number by another */
1046 /* */
1047 /* This computes C = A / B */
1048 /* */
1049 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
1050 /* lhs is A */
1051 /* rhs is B */
1052 /* set is the context */
1053 /* */
1054 /* C must have space for set->digits digits. */
1055 /* ------------------------------------------------------------------ */
1056 decNumber * decNumberDivide(decNumber *res, const decNumber *lhs,
1057 const decNumber *rhs, decContext *set) {
1058 uInt status=0; /* accumulator */
1059 decDivideOp(res, lhs, rhs, set, DIVIDE, &status);
1060 if (status!=0) decStatus(res, status, set);
1061 #if DECCHECK
1062 decCheckInexact(res, set);
1063 #endif
1064 return res;
1065 } /* decNumberDivide */
1066
1067 /* ------------------------------------------------------------------ */
1068 /* decNumberDivideInteger -- divide and return integer quotient */
1069 /* */
1070 /* This computes C = A # B, where # is the integer divide operator */
1071 /* */
1072 /* res is C, the result. C may be A and/or B (e.g., X=X#X) */
1073 /* lhs is A */
1074 /* rhs is B */
1075 /* set is the context */
1076 /* */
1077 /* C must have space for set->digits digits. */
1078 /* ------------------------------------------------------------------ */
1079 decNumber * decNumberDivideInteger(decNumber *res, const decNumber *lhs,
1080 const decNumber *rhs, decContext *set) {
1081 uInt status=0; /* accumulator */
1082 decDivideOp(res, lhs, rhs, set, DIVIDEINT, &status);
1083 if (status!=0) decStatus(res, status, set);
1084 return res;
1085 } /* decNumberDivideInteger */
1086
1087 /* ------------------------------------------------------------------ */
1088 /* decNumberExp -- exponentiation */
1089 /* */
1090 /* This computes C = exp(A) */
1091 /* */
1092 /* res is C, the result. C may be A */
1093 /* rhs is A */
1094 /* set is the context; note that rounding mode has no effect */
1095 /* */
1096 /* C must have space for set->digits digits. */
1097 /* */
1098 /* Mathematical function restrictions apply (see above); a NaN is */
1099 /* returned with Invalid_operation if a restriction is violated. */
1100 /* */
1101 /* Finite results will always be full precision and Inexact, except */
1102 /* when A is a zero or -Infinity (giving 1 or 0 respectively). */
1103 /* */
1104 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
1105 /* almost always be correctly rounded, but may be up to 1 ulp in */
1106 /* error in rare cases. */
1107 /* ------------------------------------------------------------------ */
1108 /* This is a wrapper for decExpOp which can handle the slightly wider */
1109 /* (double) range needed by Ln (which has to be able to calculate */
1110 /* exp(-a) where a can be the tiniest number (Ntiny). */
1111 /* ------------------------------------------------------------------ */
1112 decNumber * decNumberExp(decNumber *res, const decNumber *rhs,
1113 decContext *set) {
1114 uInt status=0; /* accumulator */
1115 #if DECSUBSET
1116 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */
1117 #endif
1118
1119 #if DECCHECK
1120 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1121 #endif
1122
1123 /* Check restrictions; these restrictions ensure that if h=8 (see */
1124 /* decExpOp) then the result will either overflow or underflow to 0. */
1125 /* Other math functions restrict the input range, too, for inverses. */
1126 /* If not violated then carry out the operation. */
1127 if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */
1128 #if DECSUBSET
1129 if (!set->extended) {
1130 /* reduce operand and set lostDigits status, as needed */
1131 if (rhs->digits>set->digits) {
1132 allocrhs=decRoundOperand(rhs, set, &status);
1133 if (allocrhs==NULL) break;
1134 rhs=allocrhs;
1135 }
1136 }
1137 #endif
1138 decExpOp(res, rhs, set, &status);
1139 } while(0); /* end protected */
1140
1141 #if DECSUBSET
1142 if (allocrhs !=NULL) free(allocrhs); /* drop any storage used */
1143 #endif
1144 /* apply significant status */
1145 if (status!=0) decStatus(res, status, set);
1146 #if DECCHECK
1147 decCheckInexact(res, set);
1148 #endif
1149 return res;
1150 } /* decNumberExp */
1151
1152 /* ------------------------------------------------------------------ */
1153 /* decNumberFMA -- fused multiply add */
1154 /* */
1155 /* This computes D = (A * B) + C with only one rounding */
1156 /* */
1157 /* res is D, the result. D may be A or B or C (e.g., X=FMA(X,X,X)) */
1158 /* lhs is A */
1159 /* rhs is B */
1160 /* fhs is C [far hand side] */
1161 /* set is the context */
1162 /* */
1163 /* Mathematical function restrictions apply (see above); a NaN is */
1164 /* returned with Invalid_operation if a restriction is violated. */
1165 /* */
1166 /* C must have space for set->digits digits. */
1167 /* ------------------------------------------------------------------ */
1168 decNumber * decNumberFMA(decNumber *res, const decNumber *lhs,
1169 const decNumber *rhs, const decNumber *fhs,
1170 decContext *set) {
1171 uInt status=0; /* accumulator */
1172 decContext dcmul; /* context for the multiplication */
1173 uInt needbytes; /* for space calculations */
1174 decNumber bufa[D2N(DECBUFFER*2+1)];
1175 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */
1176 decNumber *acc; /* accumulator pointer */
1177 decNumber dzero; /* work */
1178
1179 #if DECCHECK
1180 if (decCheckOperands(res, lhs, rhs, set)) return res;
1181 if (decCheckOperands(res, fhs, DECUNUSED, set)) return res;
1182 #endif
1183
1184 do { /* protect allocated storage */
1185 #if DECSUBSET
1186 if (!set->extended) { /* [undefined if subset] */
1187 status|=DEC_Invalid_operation;
1188 break;}
1189 #endif
1190 /* Check math restrictions [these ensure no overflow or underflow] */
1191 if ((!decNumberIsSpecial(lhs) && decCheckMath(lhs, set, &status))
1192 || (!decNumberIsSpecial(rhs) && decCheckMath(rhs, set, &status))
1193 || (!decNumberIsSpecial(fhs) && decCheckMath(fhs, set, &status))) break;
1194 /* set up context for multiply */
1195 dcmul=*set;
1196 dcmul.digits=lhs->digits+rhs->digits; /* just enough */
1197 /* [The above may be an over-estimate for subset arithmetic, but that's OK] */
1198 dcmul.emax=DEC_MAX_EMAX; /* effectively unbounded .. */
1199 dcmul.emin=DEC_MIN_EMIN; /* [thanks to Math restrictions] */
1200 /* set up decNumber space to receive the result of the multiply */
1201 acc=bufa; /* may fit */
1202 needbytes=sizeof(decNumber)+(D2U(dcmul.digits)-1)*sizeof(Unit);
1203 if (needbytes>sizeof(bufa)) { /* need malloc space */
1204 allocbufa=(decNumber *)malloc(needbytes);
1205 if (allocbufa==NULL) { /* hopeless -- abandon */
1206 status|=DEC_Insufficient_storage;
1207 break;}
1208 acc=allocbufa; /* use the allocated space */
1209 }
1210 /* multiply with extended range and necessary precision */
1211 /*printf("emin=%ld\n", dcmul.emin); */
1212 decMultiplyOp(acc, lhs, rhs, &dcmul, &status);
1213 /* Only Invalid operation (from sNaN or Inf * 0) is possible in */
1214 /* status; if either is seen than ignore fhs (in case it is */
1215 /* another sNaN) and set acc to NaN unless we had an sNaN */
1216 /* [decMultiplyOp leaves that to caller] */
1217 /* Note sNaN has to go through addOp to shorten payload if */
1218 /* necessary */
1219 if ((status&DEC_Invalid_operation)!=0) {
1220 if (!(status&DEC_sNaN)) { /* but be true invalid */
1221 decNumberZero(res); /* acc not yet set */
1222 res->bits=DECNAN;
1223 break;
1224 }
1225 decNumberZero(&dzero); /* make 0 (any non-NaN would do) */
1226 fhs=&dzero; /* use that */
1227 }
1228 #if DECCHECK
1229 else { /* multiply was OK */
1230 if (status!=0) printf("Status=%08lx after FMA multiply\n", status);
1231 }
1232 #endif
1233 /* add the third operand and result -> res, and all is done */
1234 decAddOp(res, acc, fhs, set, 0, &status);
1235 } while(0); /* end protected */
1236
1237 if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
1238 if (status!=0) decStatus(res, status, set);
1239 #if DECCHECK
1240 decCheckInexact(res, set);
1241 #endif
1242 return res;
1243 } /* decNumberFMA */
1244
1245 /* ------------------------------------------------------------------ */
1246 /* decNumberInvert -- invert a Number, digitwise */
1247 /* */
1248 /* This computes C = ~A */
1249 /* */
1250 /* res is C, the result. C may be A (e.g., X=~X) */
1251 /* rhs is A */
1252 /* set is the context (used for result length and error report) */
1253 /* */
1254 /* C must have space for set->digits digits. */
1255 /* */
1256 /* Logical function restrictions apply (see above); a NaN is */
1257 /* returned with Invalid_operation if a restriction is violated. */
1258 /* ------------------------------------------------------------------ */
1259 decNumber * decNumberInvert(decNumber *res, const decNumber *rhs,
1260 decContext *set) {
1261 const Unit *ua, *msua; /* -> operand and its msu */
1262 Unit *uc, *msuc; /* -> result and its msu */
1263 Int msudigs; /* digits in res msu */
1264 #if DECCHECK
1265 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1266 #endif
1267
1268 if (rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
1269 decStatus(res, DEC_Invalid_operation, set);
1270 return res;
1271 }
1272 /* operand is valid */
1273 ua=rhs->lsu; /* bottom-up */
1274 uc=res->lsu; /* .. */
1275 msua=ua+D2U(rhs->digits)-1; /* -> msu of rhs */
1276 msuc=uc+D2U(set->digits)-1; /* -> msu of result */
1277 msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */
1278 for (; uc<=msuc; ua++, uc++) { /* Unit loop */
1279 Unit a; /* extract unit */
1280 Int i, j; /* work */
1281 if (ua>msua) a=0;
1282 else a=*ua;
1283 *uc=0; /* can now write back */
1284 /* always need to examine all bits in rhs */
1285 /* This loop could be unrolled and/or use BIN2BCD tables */
1286 for (i=0; i<DECDPUN; i++) {
1287 if ((~a)&1) *uc=*uc+(Unit)powers[i]; /* effect INVERT */
1288 j=a%10;
1289 a=a/10;
1290 if (j>1) {
1291 decStatus(res, DEC_Invalid_operation, set);
1292 return res;
1293 }
1294 if (uc==msuc && i==msudigs-1) break; /* just did final digit */
1295 } /* each digit */
1296 } /* each unit */
1297 /* [here uc-1 is the msu of the result] */
1298 res->digits=decGetDigits(res->lsu, uc-res->lsu);
1299 res->exponent=0; /* integer */
1300 res->bits=0; /* sign=0 */
1301 return res; /* [no status to set] */
1302 } /* decNumberInvert */
1303
1304 /* ------------------------------------------------------------------ */
1305 /* decNumberLn -- natural logarithm */
1306 /* */
1307 /* This computes C = ln(A) */
1308 /* */
1309 /* res is C, the result. C may be A */
1310 /* rhs is A */
1311 /* set is the context; note that rounding mode has no effect */
1312 /* */
1313 /* C must have space for set->digits digits. */
1314 /* */
1315 /* Notable cases: */
1316 /* A<0 -> Invalid */
1317 /* A=0 -> -Infinity (Exact) */
1318 /* A=+Infinity -> +Infinity (Exact) */
1319 /* A=1 exactly -> 0 (Exact) */
1320 /* */
1321 /* Mathematical function restrictions apply (see above); a NaN is */
1322 /* returned with Invalid_operation if a restriction is violated. */
1323 /* */
1324 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
1325 /* almost always be correctly rounded, but may be up to 1 ulp in */
1326 /* error in rare cases. */
1327 /* ------------------------------------------------------------------ */
1328 /* This is a wrapper for decLnOp which can handle the slightly wider */
1329 /* (+11) range needed by Ln, Log10, etc. (which may have to be able */
1330 /* to calculate at p+e+2). */
1331 /* ------------------------------------------------------------------ */
1332 decNumber * decNumberLn(decNumber *res, const decNumber *rhs,
1333 decContext *set) {
1334 uInt status=0; /* accumulator */
1335 #if DECSUBSET
1336 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */
1337 #endif
1338
1339 #if DECCHECK
1340 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1341 #endif
1342
1343 /* Check restrictions; this is a math function; if not violated */
1344 /* then carry out the operation. */
1345 if (!decCheckMath(rhs, set, &status)) do { /* protect allocation */
1346 #if DECSUBSET
1347 if (!set->extended) {
1348 /* reduce operand and set lostDigits status, as needed */
1349 if (rhs->digits>set->digits) {
1350 allocrhs=decRoundOperand(rhs, set, &status);
1351 if (allocrhs==NULL) break;
1352 rhs=allocrhs;
1353 }
1354 /* special check in subset for rhs=0 */
1355 if (ISZERO(rhs)) { /* +/- zeros -> error */
1356 status|=DEC_Invalid_operation;
1357 break;}
1358 } /* extended=0 */
1359 #endif
1360 decLnOp(res, rhs, set, &status);
1361 } while(0); /* end protected */
1362
1363 #if DECSUBSET
1364 if (allocrhs !=NULL) free(allocrhs); /* drop any storage used */
1365 #endif
1366 /* apply significant status */
1367 if (status!=0) decStatus(res, status, set);
1368 #if DECCHECK
1369 decCheckInexact(res, set);
1370 #endif
1371 return res;
1372 } /* decNumberLn */
1373
1374 /* ------------------------------------------------------------------ */
1375 /* decNumberLogB - get adjusted exponent, by 754r rules */
1376 /* */
1377 /* This computes C = adjustedexponent(A) */
1378 /* */
1379 /* res is C, the result. C may be A */
1380 /* rhs is A */
1381 /* set is the context, used only for digits and status */
1382 /* */
1383 /* C must have space for 10 digits (A might have 10**9 digits and */
1384 /* an exponent of +999999999, or one digit and an exponent of */
1385 /* -1999999999). */
1386 /* */
1387 /* This returns the adjusted exponent of A after (in theory) padding */
1388 /* with zeros on the right to set->digits digits while keeping the */
1389 /* same value. The exponent is not limited by emin/emax. */
1390 /* */
1391 /* Notable cases: */
1392 /* A<0 -> Use |A| */
1393 /* A=0 -> -Infinity (Division by zero) */
1394 /* A=Infinite -> +Infinity (Exact) */
1395 /* A=1 exactly -> 0 (Exact) */
1396 /* NaNs are propagated as usual */
1397 /* ------------------------------------------------------------------ */
1398 decNumber * decNumberLogB(decNumber *res, const decNumber *rhs,
1399 decContext *set) {
1400 uInt status=0; /* accumulator */
1401
1402 #if DECCHECK
1403 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1404 #endif
1405
1406 /* NaNs as usual; Infinities return +Infinity; 0->oops */
1407 if (decNumberIsNaN(rhs)) decNaNs(res, rhs, NULL, set, &status);
1408 else if (decNumberIsInfinite(rhs)) decNumberCopyAbs(res, rhs);
1409 else if (decNumberIsZero(rhs)) {
1410 decNumberZero(res); /* prepare for Infinity */
1411 res->bits=DECNEG|DECINF; /* -Infinity */
1412 status|=DEC_Division_by_zero; /* as per 754r */
1413 }
1414 else { /* finite non-zero */
1415 Int ae=rhs->exponent+rhs->digits-1; /* adjusted exponent */
1416 decNumberFromInt32(res, ae); /* lay it out */
1417 }
1418
1419 if (status!=0) decStatus(res, status, set);
1420 return res;
1421 } /* decNumberLogB */
1422
1423 /* ------------------------------------------------------------------ */
1424 /* decNumberLog10 -- logarithm in base 10 */
1425 /* */
1426 /* This computes C = log10(A) */
1427 /* */
1428 /* res is C, the result. C may be A */
1429 /* rhs is A */
1430 /* set is the context; note that rounding mode has no effect */
1431 /* */
1432 /* C must have space for set->digits digits. */
1433 /* */
1434 /* Notable cases: */
1435 /* A<0 -> Invalid */
1436 /* A=0 -> -Infinity (Exact) */
1437 /* A=+Infinity -> +Infinity (Exact) */
1438 /* A=10**n (if n is an integer) -> n (Exact) */
1439 /* */
1440 /* Mathematical function restrictions apply (see above); a NaN is */
1441 /* returned with Invalid_operation if a restriction is violated. */
1442 /* */
1443 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
1444 /* almost always be correctly rounded, but may be up to 1 ulp in */
1445 /* error in rare cases. */
1446 /* ------------------------------------------------------------------ */
1447 /* This calculates ln(A)/ln(10) using appropriate precision. For */
1448 /* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the */
1449 /* requested digits and t is the number of digits in the exponent */
1450 /* (maximum 6). For ln(10) it is p + 3; this is often handled by the */
1451 /* fastpath in decLnOp. The final division is done to the requested */
1452 /* precision. */
1453 /* ------------------------------------------------------------------ */
1454 decNumber * decNumberLog10(decNumber *res, const decNumber *rhs,
1455 decContext *set) {
1456 uInt status=0, ignore=0; /* status accumulators */
1457 uInt needbytes; /* for space calculations */
1458 Int p; /* working precision */
1459 Int t; /* digits in exponent of A */
1460
1461 /* buffers for a and b working decimals */
1462 /* (adjustment calculator, same size) */
1463 decNumber bufa[D2N(DECBUFFER+2)];
1464 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */
1465 decNumber *a=bufa; /* temporary a */
1466 decNumber bufb[D2N(DECBUFFER+2)];
1467 decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */
1468 decNumber *b=bufb; /* temporary b */
1469 decNumber bufw[D2N(10)]; /* working 2-10 digit number */
1470 decNumber *w=bufw; /* .. */
1471 #if DECSUBSET
1472 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */
1473 #endif
1474
1475 decContext aset; /* working context */
1476
1477 #if DECCHECK
1478 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1479 #endif
1480
1481 /* Check restrictions; this is a math function; if not violated */
1482 /* then carry out the operation. */
1483 if (!decCheckMath(rhs, set, &status)) do { /* protect malloc */
1484 #if DECSUBSET
1485 if (!set->extended) {
1486 /* reduce operand and set lostDigits status, as needed */
1487 if (rhs->digits>set->digits) {
1488 allocrhs=decRoundOperand(rhs, set, &status);
1489 if (allocrhs==NULL) break;
1490 rhs=allocrhs;
1491 }
1492 /* special check in subset for rhs=0 */
1493 if (ISZERO(rhs)) { /* +/- zeros -> error */
1494 status|=DEC_Invalid_operation;
1495 break;}
1496 } /* extended=0 */
1497 #endif
1498
1499 decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */
1500
1501 /* handle exact powers of 10; only check if +ve finite */
1502 if (!(rhs->bits&(DECNEG|DECSPECIAL)) && !ISZERO(rhs)) {
1503 Int residue=0; /* (no residue) */
1504 uInt copystat=0; /* clean status */
1505
1506 /* round to a single digit... */
1507 aset.digits=1;
1508 decCopyFit(w, rhs, &aset, &residue, &copystat); /* copy & shorten */
1509 /* if exact and the digit is 1, rhs is a power of 10 */
1510 if (!(copystat&DEC_Inexact) && w->lsu[0]==1) {
1511 /* the exponent, conveniently, is the power of 10; making */
1512 /* this the result needs a little care as it might not fit, */
1513 /* so first convert it into the working number, and then move */
1514 /* to res */
1515 decNumberFromInt32(w, w->exponent);
1516 residue=0;
1517 decCopyFit(res, w, set, &residue, &status); /* copy & round */
1518 decFinish(res, set, &residue, &status); /* cleanup/set flags */
1519 break;
1520 } /* not a power of 10 */
1521 } /* not a candidate for exact */
1522
1523 /* simplify the information-content calculation to use 'total */
1524 /* number of digits in a, including exponent' as compared to the */
1525 /* requested digits, as increasing this will only rarely cost an */
1526 /* iteration in ln(a) anyway */
1527 t=6; /* it can never be >6 */
1528
1529 /* allocate space when needed... */
1530 p=(rhs->digits+t>set->digits?rhs->digits+t:set->digits)+3;
1531 needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
1532 if (needbytes>sizeof(bufa)) { /* need malloc space */
1533 allocbufa=(decNumber *)malloc(needbytes);
1534 if (allocbufa==NULL) { /* hopeless -- abandon */
1535 status|=DEC_Insufficient_storage;
1536 break;}
1537 a=allocbufa; /* use the allocated space */
1538 }
1539 aset.digits=p; /* as calculated */
1540 aset.emax=DEC_MAX_MATH; /* usual bounds */
1541 aset.emin=-DEC_MAX_MATH; /* .. */
1542 aset.clamp=0; /* and no concrete format */
1543 decLnOp(a, rhs, &aset, &status); /* a=ln(rhs) */
1544
1545 /* skip the division if the result so far is infinite, NaN, or */
1546 /* zero, or there was an error; note NaN from sNaN needs copy */
1547 if (status&DEC_NaNs && !(status&DEC_sNaN)) break;
1548 if (a->bits&DECSPECIAL || ISZERO(a)) {
1549 decNumberCopy(res, a); /* [will fit] */
1550 break;}
1551
1552 /* for ln(10) an extra 3 digits of precision are needed */
1553 p=set->digits+3;
1554 needbytes=sizeof(decNumber)+(D2U(p)-1)*sizeof(Unit);
1555 if (needbytes>sizeof(bufb)) { /* need malloc space */
1556 allocbufb=(decNumber *)malloc(needbytes);
1557 if (allocbufb==NULL) { /* hopeless -- abandon */
1558 status|=DEC_Insufficient_storage;
1559 break;}
1560 b=allocbufb; /* use the allocated space */
1561 }
1562 decNumberZero(w); /* set up 10... */
1563 #if DECDPUN==1
1564 w->lsu[1]=1; w->lsu[0]=0; /* .. */
1565 #else
1566 w->lsu[0]=10; /* .. */
1567 #endif
1568 w->digits=2; /* .. */
1569
1570 aset.digits=p;
1571 decLnOp(b, w, &aset, &ignore); /* b=ln(10) */
1572
1573 aset.digits=set->digits; /* for final divide */
1574 decDivideOp(res, a, b, &aset, DIVIDE, &status); /* into result */
1575 } while(0); /* [for break] */
1576
1577 if (allocbufa!=NULL) free(allocbufa); /* drop any storage used */
1578 if (allocbufb!=NULL) free(allocbufb); /* .. */
1579 #if DECSUBSET
1580 if (allocrhs !=NULL) free(allocrhs); /* .. */
1581 #endif
1582 /* apply significant status */
1583 if (status!=0) decStatus(res, status, set);
1584 #if DECCHECK
1585 decCheckInexact(res, set);
1586 #endif
1587 return res;
1588 } /* decNumberLog10 */
1589
1590 /* ------------------------------------------------------------------ */
1591 /* decNumberMax -- compare two Numbers and return the maximum */
1592 /* */
1593 /* This computes C = A ? B, returning the maximum by 754R rules */
1594 /* */
1595 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1596 /* lhs is A */
1597 /* rhs is B */
1598 /* set is the context */
1599 /* */
1600 /* C must have space for set->digits digits. */
1601 /* ------------------------------------------------------------------ */
1602 decNumber * decNumberMax(decNumber *res, const decNumber *lhs,
1603 const decNumber *rhs, decContext *set) {
1604 uInt status=0; /* accumulator */
1605 decCompareOp(res, lhs, rhs, set, COMPMAX, &status);
1606 if (status!=0) decStatus(res, status, set);
1607 #if DECCHECK
1608 decCheckInexact(res, set);
1609 #endif
1610 return res;
1611 } /* decNumberMax */
1612
1613 /* ------------------------------------------------------------------ */
1614 /* decNumberMaxMag -- compare and return the maximum by magnitude */
1615 /* */
1616 /* This computes C = A ? B, returning the maximum by 754R rules */
1617 /* */
1618 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1619 /* lhs is A */
1620 /* rhs is B */
1621 /* set is the context */
1622 /* */
1623 /* C must have space for set->digits digits. */
1624 /* ------------------------------------------------------------------ */
1625 decNumber * decNumberMaxMag(decNumber *res, const decNumber *lhs,
1626 const decNumber *rhs, decContext *set) {
1627 uInt status=0; /* accumulator */
1628 decCompareOp(res, lhs, rhs, set, COMPMAXMAG, &status);
1629 if (status!=0) decStatus(res, status, set);
1630 #if DECCHECK
1631 decCheckInexact(res, set);
1632 #endif
1633 return res;
1634 } /* decNumberMaxMag */
1635
1636 /* ------------------------------------------------------------------ */
1637 /* decNumberMin -- compare two Numbers and return the minimum */
1638 /* */
1639 /* This computes C = A ? B, returning the minimum by 754R rules */
1640 /* */
1641 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1642 /* lhs is A */
1643 /* rhs is B */
1644 /* set is the context */
1645 /* */
1646 /* C must have space for set->digits digits. */
1647 /* ------------------------------------------------------------------ */
1648 decNumber * decNumberMin(decNumber *res, const decNumber *lhs,
1649 const decNumber *rhs, decContext *set) {
1650 uInt status=0; /* accumulator */
1651 decCompareOp(res, lhs, rhs, set, COMPMIN, &status);
1652 if (status!=0) decStatus(res, status, set);
1653 #if DECCHECK
1654 decCheckInexact(res, set);
1655 #endif
1656 return res;
1657 } /* decNumberMin */
1658
1659 /* ------------------------------------------------------------------ */
1660 /* decNumberMinMag -- compare and return the minimum by magnitude */
1661 /* */
1662 /* This computes C = A ? B, returning the minimum by 754R rules */
1663 /* */
1664 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1665 /* lhs is A */
1666 /* rhs is B */
1667 /* set is the context */
1668 /* */
1669 /* C must have space for set->digits digits. */
1670 /* ------------------------------------------------------------------ */
1671 decNumber * decNumberMinMag(decNumber *res, const decNumber *lhs,
1672 const decNumber *rhs, decContext *set) {
1673 uInt status=0; /* accumulator */
1674 decCompareOp(res, lhs, rhs, set, COMPMINMAG, &status);
1675 if (status!=0) decStatus(res, status, set);
1676 #if DECCHECK
1677 decCheckInexact(res, set);
1678 #endif
1679 return res;
1680 } /* decNumberMinMag */
1681
1682 /* ------------------------------------------------------------------ */
1683 /* decNumberMinus -- prefix minus operator */
1684 /* */
1685 /* This computes C = 0 - A */
1686 /* */
1687 /* res is C, the result. C may be A */
1688 /* rhs is A */
1689 /* set is the context */
1690 /* */
1691 /* See also decNumberCopyNegate for a quiet bitwise version of this. */
1692 /* C must have space for set->digits digits. */
1693 /* ------------------------------------------------------------------ */
1694 /* Simply use AddOp for the subtract, which will do the necessary. */
1695 /* ------------------------------------------------------------------ */
1696 decNumber * decNumberMinus(decNumber *res, const decNumber *rhs,
1697 decContext *set) {
1698 decNumber dzero;
1699 uInt status=0; /* accumulator */
1700
1701 #if DECCHECK
1702 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1703 #endif
1704
1705 decNumberZero(&dzero); /* make 0 */
1706 dzero.exponent=rhs->exponent; /* [no coefficient expansion] */
1707 decAddOp(res, &dzero, rhs, set, DECNEG, &status);
1708 if (status!=0) decStatus(res, status, set);
1709 #if DECCHECK
1710 decCheckInexact(res, set);
1711 #endif
1712 return res;
1713 } /* decNumberMinus */
1714
1715 /* ------------------------------------------------------------------ */
1716 /* decNumberNextMinus -- next towards -Infinity */
1717 /* */
1718 /* This computes C = A - infinitesimal, rounded towards -Infinity */
1719 /* */
1720 /* res is C, the result. C may be A */
1721 /* rhs is A */
1722 /* set is the context */
1723 /* */
1724 /* This is a generalization of 754r NextDown. */
1725 /* ------------------------------------------------------------------ */
1726 decNumber * decNumberNextMinus(decNumber *res, const decNumber *rhs,
1727 decContext *set) {
1728 decNumber dtiny; /* constant */
1729 decContext workset=*set; /* work */
1730 uInt status=0; /* accumulator */
1731 #if DECCHECK
1732 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1733 #endif
1734
1735 /* +Infinity is the special case */
1736 if ((rhs->bits&(DECINF|DECNEG))==DECINF) {
1737 decSetMaxValue(res, set); /* is +ve */
1738 /* there is no status to set */
1739 return res;
1740 }
1741 decNumberZero(&dtiny); /* start with 0 */
1742 dtiny.lsu[0]=1; /* make number that is .. */
1743 dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */
1744 workset.round=DEC_ROUND_FLOOR;
1745 decAddOp(res, rhs, &dtiny, &workset, DECNEG, &status);
1746 status&=DEC_Invalid_operation|DEC_sNaN; /* only sNaN Invalid please */
1747 if (status!=0) decStatus(res, status, set);
1748 return res;
1749 } /* decNumberNextMinus */
1750
1751 /* ------------------------------------------------------------------ */
1752 /* decNumberNextPlus -- next towards +Infinity */
1753 /* */
1754 /* This computes C = A + infinitesimal, rounded towards +Infinity */
1755 /* */
1756 /* res is C, the result. C may be A */
1757 /* rhs is A */
1758 /* set is the context */
1759 /* */
1760 /* This is a generalization of 754r NextUp. */
1761 /* ------------------------------------------------------------------ */
1762 decNumber * decNumberNextPlus(decNumber *res, const decNumber *rhs,
1763 decContext *set) {
1764 decNumber dtiny; /* constant */
1765 decContext workset=*set; /* work */
1766 uInt status=0; /* accumulator */
1767 #if DECCHECK
1768 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1769 #endif
1770
1771 /* -Infinity is the special case */
1772 if ((rhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) {
1773 decSetMaxValue(res, set);
1774 res->bits=DECNEG; /* negative */
1775 /* there is no status to set */
1776 return res;
1777 }
1778 decNumberZero(&dtiny); /* start with 0 */
1779 dtiny.lsu[0]=1; /* make number that is .. */
1780 dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */
1781 workset.round=DEC_ROUND_CEILING;
1782 decAddOp(res, rhs, &dtiny, &workset, 0, &status);
1783 status&=DEC_Invalid_operation|DEC_sNaN; /* only sNaN Invalid please */
1784 if (status!=0) decStatus(res, status, set);
1785 return res;
1786 } /* decNumberNextPlus */
1787
1788 /* ------------------------------------------------------------------ */
1789 /* decNumberNextToward -- next towards rhs */
1790 /* */
1791 /* This computes C = A +/- infinitesimal, rounded towards */
1792 /* +/-Infinity in the direction of B, as per 754r nextafter rules */
1793 /* */
1794 /* res is C, the result. C may be A or B. */
1795 /* lhs is A */
1796 /* rhs is B */
1797 /* set is the context */
1798 /* */
1799 /* This is a generalization of 754r NextAfter. */
1800 /* ------------------------------------------------------------------ */
1801 decNumber * decNumberNextToward(decNumber *res, const decNumber *lhs,
1802 const decNumber *rhs, decContext *set) {
1803 decNumber dtiny; /* constant */
1804 decContext workset=*set; /* work */
1805 Int result; /* .. */
1806 uInt status=0; /* accumulator */
1807 #if DECCHECK
1808 if (decCheckOperands(res, lhs, rhs, set)) return res;
1809 #endif
1810
1811 if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) {
1812 decNaNs(res, lhs, rhs, set, &status);
1813 }
1814 else { /* Is numeric, so no chance of sNaN Invalid, etc. */
1815 result=decCompare(lhs, rhs, 0); /* sign matters */
1816 if (result==BADINT) status|=DEC_Insufficient_storage; /* rare */
1817 else { /* valid compare */
1818 if (result==0) decNumberCopySign(res, lhs, rhs); /* easy */
1819 else { /* differ: need NextPlus or NextMinus */
1820 uByte sub; /* add or subtract */
1821 if (result<0) { /* lhs<rhs, do nextplus */
1822 /* -Infinity is the special case */
1823 if ((lhs->bits&(DECINF|DECNEG))==(DECINF|DECNEG)) {
1824 decSetMaxValue(res, set);
1825 res->bits=DECNEG; /* negative */
1826 return res; /* there is no status to set */
1827 }
1828 workset.round=DEC_ROUND_CEILING;
1829 sub=0; /* add, please */
1830 } /* plus */
1831 else { /* lhs>rhs, do nextminus */
1832 /* +Infinity is the special case */
1833 if ((lhs->bits&(DECINF|DECNEG))==DECINF) {
1834 decSetMaxValue(res, set);
1835 return res; /* there is no status to set */
1836 }
1837 workset.round=DEC_ROUND_FLOOR;
1838 sub=DECNEG; /* subtract, please */
1839 } /* minus */
1840 decNumberZero(&dtiny); /* start with 0 */
1841 dtiny.lsu[0]=1; /* make number that is .. */
1842 dtiny.exponent=DEC_MIN_EMIN-1; /* .. smaller than tiniest */
1843 decAddOp(res, lhs, &dtiny, &workset, sub, &status); /* + or - */
1844 /* turn off exceptions if the result is a normal number */
1845 /* (including Nmin), otherwise let all status through */
1846 if (decNumberIsNormal(res, set)) status=0;
1847 } /* unequal */
1848 } /* compare OK */
1849 } /* numeric */
1850 if (status!=0) decStatus(res, status, set);
1851 return res;
1852 } /* decNumberNextToward */
1853
1854 /* ------------------------------------------------------------------ */
1855 /* decNumberOr -- OR two Numbers, digitwise */
1856 /* */
1857 /* This computes C = A | B */
1858 /* */
1859 /* res is C, the result. C may be A and/or B (e.g., X=X|X) */
1860 /* lhs is A */
1861 /* rhs is B */
1862 /* set is the context (used for result length and error report) */
1863 /* */
1864 /* C must have space for set->digits digits. */
1865 /* */
1866 /* Logical function restrictions apply (see above); a NaN is */
1867 /* returned with Invalid_operation if a restriction is violated. */
1868 /* ------------------------------------------------------------------ */
1869 decNumber * decNumberOr(decNumber *res, const decNumber *lhs,
1870 const decNumber *rhs, decContext *set) {
1871 const Unit *ua, *ub; /* -> operands */
1872 const Unit *msua, *msub; /* -> operand msus */
1873 Unit *uc, *msuc; /* -> result and its msu */
1874 Int msudigs; /* digits in res msu */
1875 #if DECCHECK
1876 if (decCheckOperands(res, lhs, rhs, set)) return res;
1877 #endif
1878
1879 if (lhs->exponent!=0 || decNumberIsSpecial(lhs) || decNumberIsNegative(lhs)
1880 || rhs->exponent!=0 || decNumberIsSpecial(rhs) || decNumberIsNegative(rhs)) {
1881 decStatus(res, DEC_Invalid_operation, set);
1882 return res;
1883 }
1884 /* operands are valid */
1885 ua=lhs->lsu; /* bottom-up */
1886 ub=rhs->lsu; /* .. */
1887 uc=res->lsu; /* .. */
1888 msua=ua+D2U(lhs->digits)-1; /* -> msu of lhs */
1889 msub=ub+D2U(rhs->digits)-1; /* -> msu of rhs */
1890 msuc=uc+D2U(set->digits)-1; /* -> msu of result */
1891 msudigs=MSUDIGITS(set->digits); /* [faster than remainder] */
1892 for (; uc<=msuc; ua++, ub++, uc++) { /* Unit loop */
1893 Unit a, b; /* extract units */
1894 if (ua>msua) a=0;
1895 else a=*ua;
1896 if (ub>msub) b=0;
1897 else b=*ub;
1898 *uc=0; /* can now write back */
1899 if (a|b) { /* maybe 1 bits to examine */
1900 Int i, j;
1901 /* This loop could be unrolled and/or use BIN2BCD tables */
1902 for (i=0; i<DECDPUN; i++) {
1903 if ((a|b)&1) *uc=*uc+(Unit)powers[i]; /* effect OR */
1904 j=a%10;
1905 a=a/10;
1906 j|=b%10;
1907 b=b/10;
1908 if (j>1) {
1909 decStatus(res, DEC_Invalid_operation, set);
1910 return res;
1911 }
1912 if (uc==msuc && i==msudigs-1) break; /* just did final digit */
1913 } /* each digit */
1914 } /* non-zero */
1915 } /* each unit */
1916 /* [here uc-1 is the msu of the result] */
1917 res->digits=decGetDigits(res->lsu, uc-res->lsu);
1918 res->exponent=0; /* integer */
1919 res->bits=0; /* sign=0 */
1920 return res; /* [no status to set] */
1921 } /* decNumberOr */
1922
1923 /* ------------------------------------------------------------------ */
1924 /* decNumberPlus -- prefix plus operator */
1925 /* */
1926 /* This computes C = 0 + A */
1927 /* */
1928 /* res is C, the result. C may be A */
1929 /* rhs is A */
1930 /* set is the context */
1931 /* */
1932 /* See also decNumberCopy for a quiet bitwise version of this. */
1933 /* C must have space for set->digits digits. */
1934 /* ------------------------------------------------------------------ */
1935 /* This simply uses AddOp; Add will take fast path after preparing A. */
1936 /* Performance is a concern here, as this routine is often used to */
1937 /* check operands and apply rounding and overflow/underflow testing. */
1938 /* ------------------------------------------------------------------ */
1939 decNumber * decNumberPlus(decNumber *res, const decNumber *rhs,
1940 decContext *set) {
1941 decNumber dzero;
1942 uInt status=0; /* accumulator */
1943 #if DECCHECK
1944 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
1945 #endif
1946
1947 decNumberZero(&dzero); /* make 0 */
1948 dzero.exponent=rhs->exponent; /* [no coefficient expansion] */
1949 decAddOp(res, &dzero, rhs, set, 0, &status);
1950 if (status!=0) decStatus(res, status, set);
1951 #if DECCHECK
1952 decCheckInexact(res, set);
1953 #endif
1954 return res;
1955 } /* decNumberPlus */
1956
1957 /* ------------------------------------------------------------------ */
1958 /* decNumberMultiply -- multiply two Numbers */
1959 /* */
1960 /* This computes C = A x B */
1961 /* */
1962 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
1963 /* lhs is A */
1964 /* rhs is B */
1965 /* set is the context */
1966 /* */
1967 /* C must have space for set->digits digits. */
1968 /* ------------------------------------------------------------------ */
1969 decNumber * decNumberMultiply(decNumber *res, const decNumber *lhs,
1970 const decNumber *rhs, decContext *set) {
1971 uInt status=0; /* accumulator */
1972 decMultiplyOp(res, lhs, rhs, set, &status);
1973 if (status!=0) decStatus(res, status, set);
1974 #if DECCHECK
1975 decCheckInexact(res, set);
1976 #endif
1977 return res;
1978 } /* decNumberMultiply */
1979
1980 /* ------------------------------------------------------------------ */
1981 /* decNumberPower -- raise a number to a power */
1982 /* */
1983 /* This computes C = A ** B */
1984 /* */
1985 /* res is C, the result. C may be A and/or B (e.g., X=X**X) */
1986 /* lhs is A */
1987 /* rhs is B */
1988 /* set is the context */
1989 /* */
1990 /* C must have space for set->digits digits. */
1991 /* */
1992 /* Mathematical function restrictions apply (see above); a NaN is */
1993 /* returned with Invalid_operation if a restriction is violated. */
1994 /* */
1995 /* However, if 1999999997<=B<=999999999 and B is an integer then the */
1996 /* restrictions on A and the context are relaxed to the usual bounds, */
1997 /* for compatibility with the earlier (integer power only) version */
1998 /* of this function. */
1999 /* */
2000 /* When B is an integer, the result may be exact, even if rounded. */
2001 /* */
2002 /* The final result is rounded according to the context; it will */
2003 /* almost always be correctly rounded, but may be up to 1 ulp in */
2004 /* error in rare cases. */
2005 /* ------------------------------------------------------------------ */
2006 decNumber * decNumberPower(decNumber *res, const decNumber *lhs,
2007 const decNumber *rhs, decContext *set) {
2008 #if DECSUBSET
2009 decNumber *alloclhs=NULL; /* non-NULL if rounded lhs allocated */
2010 decNumber *allocrhs=NULL; /* .., rhs */
2011 #endif
2012 decNumber *allocdac=NULL; /* -> allocated acc buffer, iff used */
2013 decNumber *allocinv=NULL; /* -> allocated 1/x buffer, iff used */
2014 Int reqdigits=set->digits; /* requested DIGITS */
2015 Int n; /* rhs in binary */
2016 Flag rhsint=0; /* 1 if rhs is an integer */
2017 Flag useint=0; /* 1 if can use integer calculation */
2018 Flag isoddint=0; /* 1 if rhs is an integer and odd */
2019 Int i; /* work */
2020 #if DECSUBSET
2021 Int dropped; /* .. */
2022 #endif
2023 uInt needbytes; /* buffer size needed */
2024 Flag seenbit; /* seen a bit while powering */
2025 Int residue=0; /* rounding residue */
2026 uInt status=0; /* accumulators */
2027 uByte bits=0; /* result sign if errors */
2028 decContext aset; /* working context */
2029 decNumber dnOne; /* work value 1... */
2030 /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
2031 decNumber dacbuff[D2N(DECBUFFER+9)];
2032 decNumber *dac=dacbuff; /* -> result accumulator */
2033 /* same again for possible 1/lhs calculation */
2034 decNumber invbuff[D2N(DECBUFFER+9)];
2035
2036 #if DECCHECK
2037 if (decCheckOperands(res, lhs, rhs, set)) return res;
2038 #endif
2039
2040 do { /* protect allocated storage */
2041 #if DECSUBSET
2042 if (!set->extended) { /* reduce operands and set status, as needed */
2043 if (lhs->digits>reqdigits) {
2044 alloclhs=decRoundOperand(lhs, set, &status);
2045 if (alloclhs==NULL) break;
2046 lhs=alloclhs;
2047 }
2048 if (rhs->digits>reqdigits) {
2049 allocrhs=decRoundOperand(rhs, set, &status);
2050 if (allocrhs==NULL) break;
2051 rhs=allocrhs;
2052 }
2053 }
2054 #endif
2055 /* [following code does not require input rounding] */
2056
2057 /* handle NaNs and rhs Infinity (lhs infinity is harder) */
2058 if (SPECIALARGS) {
2059 if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs)) { /* NaNs */
2060 decNaNs(res, lhs, rhs, set, &status);
2061 break;}
2062 if (decNumberIsInfinite(rhs)) { /* rhs Infinity */
2063 Flag rhsneg=rhs->bits&DECNEG; /* save rhs sign */
2064 if (decNumberIsNegative(lhs) /* lhs<0 */
2065 && !decNumberIsZero(lhs)) /* .. */
2066 status|=DEC_Invalid_operation;
2067 else { /* lhs >=0 */
2068 decNumberZero(&dnOne); /* set up 1 */
2069 dnOne.lsu[0]=1;
2070 decNumberCompare(dac, lhs, &dnOne, set); /* lhs ? 1 */
2071 decNumberZero(res); /* prepare for 0/1/Infinity */
2072 if (decNumberIsNegative(dac)) { /* lhs<1 */
2073 if (rhsneg) res->bits|=DECINF; /* +Infinity [else is +0] */
2074 }
2075 else if (dac->lsu[0]==0) { /* lhs=1 */
2076 /* 1**Infinity is inexact, so return fully-padded 1.0000 */
2077 Int shift=set->digits-1;
2078 *res->lsu=1; /* was 0, make int 1 */
2079 res->digits=decShiftToMost(res->lsu, 1, shift);
2080 res->exponent=-shift; /* make 1.0000... */
2081 status|=DEC_Inexact|DEC_Rounded; /* deemed inexact */
2082 }
2083 else { /* lhs>1 */
2084 if (!rhsneg) res->bits|=DECINF; /* +Infinity [else is +0] */
2085 }
2086 } /* lhs>=0 */
2087 break;}
2088 /* [lhs infinity drops through] */
2089 } /* specials */
2090
2091 /* Original rhs may be an integer that fits and is in range */
2092 n=decGetInt(rhs);
2093 if (n!=BADINT) { /* it is an integer */
2094 rhsint=1; /* record the fact for 1**n */
2095 isoddint=(Flag)n&1; /* [works even if big] */
2096 if (n!=BIGEVEN && n!=BIGODD) /* can use integer path? */
2097 useint=1; /* looks good */
2098 }
2099
2100 if (decNumberIsNegative(lhs) /* -x .. */
2101 && isoddint) bits=DECNEG; /* .. to an odd power */
2102
2103 /* handle LHS infinity */
2104 if (decNumberIsInfinite(lhs)) { /* [NaNs already handled] */
2105 uByte rbits=rhs->bits; /* save */
2106 decNumberZero(res); /* prepare */
2107 if (n==0) *res->lsu=1; /* [-]Inf**0 => 1 */
2108 else {
2109 /* -Inf**nonint -> error */
2110 if (!rhsint && decNumberIsNegative(lhs)) {
2111 status|=DEC_Invalid_operation; /* -Inf**nonint is error */
2112 break;}
2113 if (!(rbits & DECNEG)) bits|=DECINF; /* was not a **-n */
2114 /* [otherwise will be 0 or -0] */
2115 res->bits=bits;
2116 }
2117 break;}
2118
2119 /* similarly handle LHS zero */
2120 if (decNumberIsZero(lhs)) {
2121 if (n==0) { /* 0**0 => Error */
2122 #if DECSUBSET
2123 if (!set->extended) { /* [unless subset] */
2124 decNumberZero(res);
2125 *res->lsu=1; /* return 1 */
2126 break;}
2127 #endif
2128 status|=DEC_Invalid_operation;
2129 }
2130 else { /* 0**x */
2131 uByte rbits=rhs->bits; /* save */
2132 if (rbits & DECNEG) { /* was a 0**(-n) */
2133 #if DECSUBSET
2134 if (!set->extended) { /* [bad if subset] */
2135 status|=DEC_Invalid_operation;
2136 break;}
2137 #endif
2138 bits|=DECINF;
2139 }
2140 decNumberZero(res); /* prepare */
2141 /* [otherwise will be 0 or -0] */
2142 res->bits=bits;
2143 }
2144 break;}
2145
2146 /* here both lhs and rhs are finite; rhs==0 is handled in the */
2147 /* integer path. Next handle the non-integer cases */
2148 if (!useint) { /* non-integral rhs */
2149 /* any -ve lhs is bad, as is either operand or context out of */
2150 /* bounds */
2151 if (decNumberIsNegative(lhs)) {
2152 status|=DEC_Invalid_operation;
2153 break;}
2154 if (decCheckMath(lhs, set, &status)
2155 || decCheckMath(rhs, set, &status)) break; /* variable status */
2156
2157 decContextDefault(&aset, DEC_INIT_DECIMAL64); /* clean context */
2158 aset.emax=DEC_MAX_MATH; /* usual bounds */
2159 aset.emin=-DEC_MAX_MATH; /* .. */
2160 aset.clamp=0; /* and no concrete format */
2161
2162 /* calculate the result using exp(ln(lhs)*rhs), which can */
2163 /* all be done into the accumulator, dac. The precision needed */
2164 /* is enough to contain the full information in the lhs (which */
2165 /* is the total digits, including exponent), or the requested */
2166 /* precision, if larger, + 4; 6 is used for the exponent */
2167 /* maximum length, and this is also used when it is shorter */
2168 /* than the requested digits as it greatly reduces the >0.5 ulp */
2169 /* cases at little cost (because Ln doubles digits each */
2170 /* iteration so a few extra digits rarely causes an extra */
2171 /* iteration) */
2172 aset.digits=MAXI(lhs->digits, set->digits)+6+4;
2173 } /* non-integer rhs */
2174
2175 else { /* rhs is in-range integer */
2176 if (n==0) { /* x**0 = 1 */
2177 /* (0**0 was handled above) */
2178 decNumberZero(res); /* result=1 */
2179 *res->lsu=1; /* .. */
2180 break;}
2181 /* rhs is a non-zero integer */
2182 if (n<0) n=-n; /* use abs(n) */
2183
2184 aset=*set; /* clone the context */
2185 aset.round=DEC_ROUND_HALF_EVEN; /* internally use balanced */
2186 /* calculate the working DIGITS */
2187 aset.digits=reqdigits+(rhs->digits+rhs->exponent)+2;
2188 #if DECSUBSET
2189 if (!set->extended) aset.digits--; /* use classic precision */
2190 #endif
2191 /* it's an error if this is more than can be handled */
2192 if (aset.digits>DECNUMMAXP) {status|=DEC_Invalid_operation; break;}
2193 } /* integer path */
2194
2195 /* aset.digits is the count of digits for the accumulator needed */
2196 /* if accumulator is too long for local storage, then allocate */
2197 needbytes=sizeof(decNumber)+(D2U(aset.digits)-1)*sizeof(Unit);
2198 /* [needbytes also used below if 1/lhs needed] */
2199 if (needbytes>sizeof(dacbuff)) {
2200 allocdac=(decNumber *)malloc(needbytes);
2201 if (allocdac==NULL) { /* hopeless -- abandon */
2202 status|=DEC_Insufficient_storage;
2203 break;}
2204 dac=allocdac; /* use the allocated space */
2205 }
2206 /* here, aset is set up and accumulator is ready for use */
2207
2208 if (!useint) { /* non-integral rhs */
2209 /* x ** y; special-case x=1 here as it will otherwise always */
2210 /* reduce to integer 1; decLnOp has a fastpath which detects */
2211 /* the case of x=1 */
2212 decLnOp(dac, lhs, &aset, &status); /* dac=ln(lhs) */
2213 /* [no error possible, as lhs 0 already handled] */
2214 if (ISZERO(dac)) { /* x==1, 1.0, etc. */
2215 /* need to return fully-padded 1.0000 etc., but rhsint->1 */
2216 *dac->lsu=1; /* was 0, make int 1 */
2217 if (!rhsint) { /* add padding */
2218 Int shift=set->digits-1;
2219 dac->digits=decShiftToMost(dac->lsu, 1, shift);
2220 dac->exponent=-shift; /* make 1.0000... */
2221 status|=DEC_Inexact|DEC_Rounded; /* deemed inexact */
2222 }
2223 }
2224 else {
2225 decMultiplyOp(dac, dac, rhs, &aset, &status); /* dac=dac*rhs */
2226 decExpOp(dac, dac, &aset, &status); /* dac=exp(dac) */
2227 }
2228 /* and drop through for final rounding */
2229 } /* non-integer rhs */
2230
2231 else { /* carry on with integer */
2232 decNumberZero(dac); /* acc=1 */
2233 *dac->lsu=1; /* .. */
2234
2235 /* if a negative power the constant 1 is needed, and if not subset */
2236 /* invert the lhs now rather than inverting the result later */
2237 if (decNumberIsNegative(rhs)) { /* was a **-n [hence digits>0] */
2238 decNumber *inv=invbuff; /* assume use fixed buffer */
2239 decNumberCopy(&dnOne, dac); /* dnOne=1; [needed now or later] */
2240 #if DECSUBSET
2241 if (set->extended) { /* need to calculate 1/lhs */
2242 #endif
2243 /* divide lhs into 1, putting result in dac [dac=1/dac] */
2244 decDivideOp(dac, &dnOne, lhs, &aset, DIVIDE, &status);
2245 /* now locate or allocate space for the inverted lhs */
2246 if (needbytes>sizeof(invbuff)) {
2247 allocinv=(decNumber *)malloc(needbytes);
2248 if (allocinv==NULL) { /* hopeless -- abandon */
2249 status|=DEC_Insufficient_storage;
2250 break;}
2251 inv=allocinv; /* use the allocated space */
2252 }
2253 /* [inv now points to big-enough buffer or allocated storage] */
2254 decNumberCopy(inv, dac); /* copy the 1/lhs */
2255 decNumberCopy(dac, &dnOne); /* restore acc=1 */
2256 lhs=inv; /* .. and go forward with new lhs */
2257 #if DECSUBSET
2258 }
2259 #endif
2260 }
2261
2262 /* Raise-to-the-power loop... */
2263 seenbit=0; /* set once a 1-bit is encountered */
2264 for (i=1;;i++){ /* for each bit [top bit ignored] */
2265 /* abandon if had overflow or terminal underflow */
2266 if (status & (DEC_Overflow|DEC_Underflow)) { /* interesting? */
2267 if (status&DEC_Overflow || ISZERO(dac)) break;
2268 }
2269 /* [the following two lines revealed an optimizer bug in a C++ */
2270 /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
2271 n=n<<1; /* move next bit to testable position */
2272 if (n<0) { /* top bit is set */
2273 seenbit=1; /* OK, significant bit seen */
2274 decMultiplyOp(dac, dac, lhs, &aset, &status); /* dac=dac*x */
2275 }
2276 if (i==31) break; /* that was the last bit */
2277 if (!seenbit) continue; /* no need to square 1 */
2278 decMultiplyOp(dac, dac, dac, &aset, &status); /* dac=dac*dac [square] */
2279 } /*i*/ /* 32 bits */
2280
2281 /* complete internal overflow or underflow processing */
2282 if (status & (DEC_Overflow|DEC_Underflow)) {
2283 #if DECSUBSET
2284 /* If subset, and power was negative, reverse the kind of -erflow */
2285 /* [1/x not yet done] */
2286 if (!set->extended && decNumberIsNegative(rhs)) {
2287 if (status & DEC_Overflow)
2288 status^=DEC_Overflow | DEC_Underflow | DEC_Subnormal;
2289 else { /* trickier -- Underflow may or may not be set */
2290 status&=~(DEC_Underflow | DEC_Subnormal); /* [one or both] */
2291 status|=DEC_Overflow;
2292 }
2293 }
2294 #endif
2295 dac->bits=(dac->bits & ~DECNEG) | bits; /* force correct sign */
2296 /* round subnormals [to set.digits rather than aset.digits] */
2297 /* or set overflow result similarly as required */
2298 decFinalize(dac, set, &residue, &status);
2299 decNumberCopy(res, dac); /* copy to result (is now OK length) */
2300 break;
2301 }
2302
2303 #if DECSUBSET
2304 if (!set->extended && /* subset math */
2305 decNumberIsNegative(rhs)) { /* was a **-n [hence digits>0] */
2306 /* so divide result into 1 [dac=1/dac] */
2307 decDivideOp(dac, &dnOne, dac, &aset, DIVIDE, &status);
2308 }
2309 #endif
2310 } /* rhs integer path */
2311
2312 /* reduce result to the requested length and copy to result */
2313 decCopyFit(res, dac, set, &residue, &status);
2314 decFinish(res, set, &residue, &status); /* final cleanup */
2315 #if DECSUBSET
2316 if (!set->extended) decTrim(res, set, 0, &dropped); /* trailing zeros */
2317 #endif
2318 } while(0); /* end protected */
2319
2320 if (allocdac!=NULL) free(allocdac); /* drop any storage used */
2321 if (allocinv!=NULL) free(allocinv); /* .. */
2322 #if DECSUBSET
2323 if (alloclhs!=NULL) free(alloclhs); /* .. */
2324 if (allocrhs!=NULL) free(allocrhs); /* .. */
2325 #endif
2326 if (status!=0) decStatus(res, status, set);
2327 #if DECCHECK
2328 decCheckInexact(res, set);
2329 #endif
2330 return res;
2331 } /* decNumberPower */
2332
2333 /* ------------------------------------------------------------------ */
2334 /* decNumberQuantize -- force exponent to requested value */
2335 /* */
2336 /* This computes C = op(A, B), where op adjusts the coefficient */
2337 /* of C (by rounding or shifting) such that the exponent (-scale) */
2338 /* of C has exponent of B. The numerical value of C will equal A, */
2339 /* except for the effects of any rounding that occurred. */
2340 /* */
2341 /* res is C, the result. C may be A or B */
2342 /* lhs is A, the number to adjust */
2343 /* rhs is B, the number with exponent to match */
2344 /* set is the context */
2345 /* */
2346 /* C must have space for set->digits digits. */
2347 /* */
2348 /* Unless there is an error or the result is infinite, the exponent */
2349 /* after the operation is guaranteed to be equal to that of B. */
2350 /* ------------------------------------------------------------------ */
2351 decNumber * decNumberQuantize(decNumber *res, const decNumber *lhs,
2352 const decNumber *rhs, decContext *set) {
2353 uInt status=0; /* accumulator */
2354 decQuantizeOp(res, lhs, rhs, set, 1, &status);
2355 if (status!=0) decStatus(res, status, set);
2356 return res;
2357 } /* decNumberQuantize */
2358
2359 /* ------------------------------------------------------------------ */
2360 /* decNumberReduce -- remove trailing zeros */
2361 /* */
2362 /* This computes C = 0 + A, and normalizes the result */
2363 /* */
2364 /* res is C, the result. C may be A */
2365 /* rhs is A */
2366 /* set is the context */
2367 /* */
2368 /* C must have space for set->digits digits. */
2369 /* ------------------------------------------------------------------ */
2370 /* Previously known as Normalize */
2371 decNumber * decNumberNormalize(decNumber *res, const decNumber *rhs,
2372 decContext *set) {
2373 return decNumberReduce(res, rhs, set);
2374 } /* decNumberNormalize */
2375
2376 decNumber * decNumberReduce(decNumber *res, const decNumber *rhs,
2377 decContext *set) {
2378 #if DECSUBSET
2379 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */
2380 #endif
2381 uInt status=0; /* as usual */
2382 Int residue=0; /* as usual */
2383 Int dropped; /* work */
2384
2385 #if DECCHECK
2386 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
2387 #endif
2388
2389 do { /* protect allocated storage */
2390 #if DECSUBSET
2391 if (!set->extended) {
2392 /* reduce operand and set lostDigits status, as needed */
2393 if (rhs->digits>set->digits) {
2394 allocrhs=decRoundOperand(rhs, set, &status);
2395 if (allocrhs==NULL) break;
2396 rhs=allocrhs;
2397 }
2398 }
2399 #endif
2400 /* [following code does not require input rounding] */
2401
2402 /* Infinities copy through; NaNs need usual treatment */
2403 if (decNumberIsNaN(rhs)) {
2404 decNaNs(res, rhs, NULL, set, &status);
2405 break;
2406 }
2407
2408 /* reduce result to the requested length and copy to result */
2409 decCopyFit(res, rhs, set, &residue, &status); /* copy & round */
2410 decFinish(res, set, &residue, &status); /* cleanup/set flags */
2411 decTrim(res, set, 1, &dropped); /* normalize in place */
2412 } while(0); /* end protected */
2413
2414 #if DECSUBSET
2415 if (allocrhs !=NULL) free(allocrhs); /* .. */
2416 #endif
2417 if (status!=0) decStatus(res, status, set);/* then report status */
2418 return res;
2419 } /* decNumberReduce */
2420
2421 /* ------------------------------------------------------------------ */
2422 /* decNumberRescale -- force exponent to requested value */
2423 /* */
2424 /* This computes C = op(A, B), where op adjusts the coefficient */
2425 /* of C (by rounding or shifting) such that the exponent (-scale) */
2426 /* of C has the value B. The numerical value of C will equal A, */
2427 /* except for the effects of any rounding that occurred. */
2428 /* */
2429 /* res is C, the result. C may be A or B */
2430 /* lhs is A, the number to adjust */
2431 /* rhs is B, the requested exponent */
2432 /* set is the context */
2433 /* */
2434 /* C must have space for set->digits digits. */
2435 /* */
2436 /* Unless there is an error or the result is infinite, the exponent */
2437 /* after the operation is guaranteed to be equal to B. */
2438 /* ------------------------------------------------------------------ */
2439 decNumber * decNumberRescale(decNumber *res, const decNumber *lhs,
2440 const decNumber *rhs, decContext *set) {
2441 uInt status=0; /* accumulator */
2442 decQuantizeOp(res, lhs, rhs, set, 0, &status);
2443 if (status!=0) decStatus(res, status, set);
2444 return res;
2445 } /* decNumberRescale */
2446
2447 /* ------------------------------------------------------------------ */
2448 /* decNumberRemainder -- divide and return remainder */
2449 /* */
2450 /* This computes C = A % B */
2451 /* */
2452 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
2453 /* lhs is A */
2454 /* rhs is B */
2455 /* set is the context */
2456 /* */
2457 /* C must have space for set->digits digits. */
2458 /* ------------------------------------------------------------------ */
2459 decNumber * decNumberRemainder(decNumber *res, const decNumber *lhs,
2460 const decNumber *rhs, decContext *set) {
2461 uInt status=0; /* accumulator */
2462 decDivideOp(res, lhs, rhs, set, REMAINDER, &status);
2463 if (status!=0) decStatus(res, status, set);
2464 #if DECCHECK
2465 decCheckInexact(res, set);
2466 #endif
2467 return res;
2468 } /* decNumberRemainder */
2469
2470 /* ------------------------------------------------------------------ */
2471 /* decNumberRemainderNear -- divide and return remainder from nearest */
2472 /* */
2473 /* This computes C = A % B, where % is the IEEE remainder operator */
2474 /* */
2475 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
2476 /* lhs is A */
2477 /* rhs is B */
2478 /* set is the context */
2479 /* */
2480 /* C must have space for set->digits digits. */
2481 /* ------------------------------------------------------------------ */
2482 decNumber * decNumberRemainderNear(decNumber *res, const decNumber *lhs,
2483 const decNumber *rhs, decContext *set) {
2484 uInt status=0; /* accumulator */
2485 decDivideOp(res, lhs, rhs, set, REMNEAR, &status);
2486 if (status!=0) decStatus(res, status, set);
2487 #if DECCHECK
2488 decCheckInexact(res, set);
2489 #endif
2490 return res;
2491 } /* decNumberRemainderNear */
2492
2493 /* ------------------------------------------------------------------ */
2494 /* decNumberRotate -- rotate the coefficient of a Number left/right */
2495 /* */
2496 /* This computes C = A rot B (in base ten and rotating set->digits */
2497 /* digits). */
2498 /* */
2499 /* res is C, the result. C may be A and/or B (e.g., X=XrotX) */
2500 /* lhs is A */
2501 /* rhs is B, the number of digits to rotate (-ve to right) */
2502 /* set is the context */
2503 /* */
2504 /* The digits of the coefficient of A are rotated to the left (if B */
2505 /* is positive) or to the right (if B is negative) without adjusting */
2506 /* the exponent or the sign of A. If lhs->digits is less than */
2507 /* set->digits the coefficient is padded with zeros on the left */
2508 /* before the rotate. Any leading zeros in the result are removed */
2509 /* as usual. */
2510 /* */
2511 /* B must be an integer (q=0) and in the range -set->digits through */
2512 /* +set->digits. */
2513 /* C must have space for set->digits digits. */
2514 /* NaNs are propagated as usual. Infinities are unaffected (but */
2515 /* B must be valid). No status is set unless B is invalid or an */
2516 /* operand is an sNaN. */
2517 /* ------------------------------------------------------------------ */
2518 decNumber * decNumberRotate(decNumber *res, const decNumber *lhs,
2519 const decNumber *rhs, decContext *set) {
2520 uInt status=0; /* accumulator */
2521 Int rotate; /* rhs as an Int */
2522
2523 #if DECCHECK
2524 if (decCheckOperands(res, lhs, rhs, set)) return res;
2525 #endif
2526
2527 /* NaNs propagate as normal */
2528 if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
2529 decNaNs(res, lhs, rhs, set, &status);
2530 /* rhs must be an integer */
2531 else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
2532 status=DEC_Invalid_operation;
2533 else { /* both numeric, rhs is an integer */
2534 rotate=decGetInt(rhs); /* [cannot fail] */
2535 if (rotate==BADINT /* something bad .. */
2536 || rotate==BIGODD || rotate==BIGEVEN /* .. very big .. */
2537 || abs(rotate)>set->digits) /* .. or out of range */
2538 status=DEC_Invalid_operation;
2539 else { /* rhs is OK */
2540 decNumberCopy(res, lhs);
2541 /* convert -ve rotate to equivalent positive rotation */
2542 if (rotate<0) rotate=set->digits+rotate;
2543 if (rotate!=0 && rotate!=set->digits /* zero or full rotation */
2544 && !decNumberIsInfinite(res)) { /* lhs was infinite */
2545 /* left-rotate to do; 0 < rotate < set->digits */
2546 uInt units, shift; /* work */
2547 uInt msudigits; /* digits in result msu */
2548 Unit *msu=res->lsu+D2U(res->digits)-1; /* current msu */
2549 Unit *msumax=res->lsu+D2U(set->digits)-1; /* rotation msu */
2550 for (msu++; msu<=msumax; msu++) *msu=0; /* ensure high units=0 */
2551 res->digits=set->digits; /* now full-length */
2552 msudigits=MSUDIGITS(res->digits); /* actual digits in msu */
2553
2554 /* rotation here is done in-place, in three steps */
2555 /* 1. shift all to least up to one unit to unit-align final */
2556 /* lsd [any digits shifted out are rotated to the left, */
2557 /* abutted to the original msd (which may require split)] */
2558 /* */
2559 /* [if there are no whole units left to rotate, the */
2560 /* rotation is now complete] */
2561 /* */
2562 /* 2. shift to least, from below the split point only, so that */
2563 /* the final msd is in the right place in its Unit [any */
2564 /* digits shifted out will fit exactly in the current msu, */
2565 /* left aligned, no split required] */
2566 /* */
2567 /* 3. rotate all the units by reversing left part, right */
2568 /* part, and then whole */
2569 /* */
2570 /* example: rotate right 8 digits (2 units + 2), DECDPUN=3. */
2571 /* */
2572 /* start: 00a bcd efg hij klm npq */
2573 /* */
2574 /* 1a 000 0ab cde fgh|ijk lmn [pq saved] */
2575 /* 1b 00p qab cde fgh|ijk lmn */
2576 /* */
2577 /* 2a 00p qab cde fgh|00i jkl [mn saved] */
2578 /* 2b mnp qab cde fgh|00i jkl */
2579 /* */
2580 /* 3a fgh cde qab mnp|00i jkl */
2581 /* 3b fgh cde qab mnp|jkl 00i */
2582 /* 3c 00i jkl mnp qab cde fgh */
2583
2584 /* Step 1: amount to shift is the partial right-rotate count */
2585 rotate=set->digits-rotate; /* make it right-rotate */
2586 units=rotate/DECDPUN; /* whole units to rotate */
2587 shift=rotate%DECDPUN; /* left-over digits count */
2588 if (shift>0) { /* not an exact number of units */
2589 uInt save=res->lsu[0]%powers[shift]; /* save low digit(s) */
2590 decShiftToLeast(res->lsu, D2U(res->digits), shift);
2591 if (shift>msudigits) { /* msumax-1 needs >0 digits */
2592 uInt rem=save%powers[shift-msudigits];/* split save */
2593 *msumax=(Unit)(save/powers[shift-msudigits]); /* and insert */
2594 *(msumax-1)=*(msumax-1)
2595 +(Unit)(rem*powers[DECDPUN-(shift-msudigits)]); /* .. */
2596 }
2597 else { /* all fits in msumax */
2598 *msumax=*msumax+(Unit)(save*powers[msudigits-shift]); /* [maybe *1] */
2599 }
2600 } /* digits shift needed */
2601
2602 /* If whole units to rotate... */
2603 if (units>0) { /* some to do */
2604 /* Step 2: the units to touch are the whole ones in rotate, */
2605 /* if any, and the shift is DECDPUN-msudigits (which may be */
2606 /* 0, again) */
2607 shift=DECDPUN-msudigits;
2608 if (shift>0) { /* not an exact number of units */
2609 uInt save=res->lsu[0]%powers[shift]; /* save low digit(s) */
2610 decShiftToLeast(res->lsu, units, shift);
2611 *msumax=*msumax+(Unit)(save*powers[msudigits]);
2612 } /* partial shift needed */
2613
2614 /* Step 3: rotate the units array using triple reverse */
2615 /* (reversing is easy and fast) */
2616 decReverse(res->lsu+units, msumax); /* left part */
2617 decReverse(res->lsu, res->lsu+units-1); /* right part */
2618 decReverse(res->lsu, msumax); /* whole */
2619 } /* whole units to rotate */
2620 /* the rotation may have left an undetermined number of zeros */
2621 /* on the left, so true length needs to be calculated */
2622 res->digits=decGetDigits(res->lsu, msumax-res->lsu+1);
2623 } /* rotate needed */
2624 } /* rhs OK */
2625 } /* numerics */
2626 if (status!=0) decStatus(res, status, set);
2627 return res;
2628 } /* decNumberRotate */
2629
2630 /* ------------------------------------------------------------------ */
2631 /* decNumberSameQuantum -- test for equal exponents */
2632 /* */
2633 /* res is the result number, which will contain either 0 or 1 */
2634 /* lhs is a number to test */
2635 /* rhs is the second (usually a pattern) */
2636 /* */
2637 /* No errors are possible and no context is needed. */
2638 /* ------------------------------------------------------------------ */
2639 decNumber * decNumberSameQuantum(decNumber *res, const decNumber *lhs,
2640 const decNumber *rhs) {
2641 Unit ret=0; /* return value */
2642
2643 #if DECCHECK
2644 if (decCheckOperands(res, lhs, rhs, DECUNCONT)) return res;
2645 #endif
2646
2647 if (SPECIALARGS) {
2648 if (decNumberIsNaN(lhs) && decNumberIsNaN(rhs)) ret=1;
2649 else if (decNumberIsInfinite(lhs) && decNumberIsInfinite(rhs)) ret=1;
2650 /* [anything else with a special gives 0] */
2651 }
2652 else if (lhs->exponent==rhs->exponent) ret=1;
2653
2654 decNumberZero(res); /* OK to overwrite an operand now */
2655 *res->lsu=ret;
2656 return res;
2657 } /* decNumberSameQuantum */
2658
2659 /* ------------------------------------------------------------------ */
2660 /* decNumberScaleB -- multiply by a power of 10 */
2661 /* */
2662 /* This computes C = A x 10**B where B is an integer (q=0) with */
2663 /* maximum magnitude 2*(emax+digits) */
2664 /* */
2665 /* res is C, the result. C may be A or B */
2666 /* lhs is A, the number to adjust */
2667 /* rhs is B, the requested power of ten to use */
2668 /* set is the context */
2669 /* */
2670 /* C must have space for set->digits digits. */
2671 /* */
2672 /* The result may underflow or overflow. */
2673 /* ------------------------------------------------------------------ */
2674 decNumber * decNumberScaleB(decNumber *res, const decNumber *lhs,
2675 const decNumber *rhs, decContext *set) {
2676 Int reqexp; /* requested exponent change [B] */
2677 uInt status=0; /* accumulator */
2678 Int residue; /* work */
2679
2680 #if DECCHECK
2681 if (decCheckOperands(res, lhs, rhs, set)) return res;
2682 #endif
2683
2684 /* Handle special values except lhs infinite */
2685 if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
2686 decNaNs(res, lhs, rhs, set, &status);
2687 /* rhs must be an integer */
2688 else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
2689 status=DEC_Invalid_operation;
2690 else {
2691 /* lhs is a number; rhs is a finite with q==0 */
2692 reqexp=decGetInt(rhs); /* [cannot fail] */
2693 if (reqexp==BADINT /* something bad .. */
2694 || reqexp==BIGODD || reqexp==BIGEVEN /* .. very big .. */
2695 || abs(reqexp)>(2*(set->digits+set->emax))) /* .. or out of range */
2696 status=DEC_Invalid_operation;
2697 else { /* rhs is OK */
2698 decNumberCopy(res, lhs); /* all done if infinite lhs */
2699 if (!decNumberIsInfinite(res)) { /* prepare to scale */
2700 res->exponent+=reqexp; /* adjust the exponent */
2701 residue=0;
2702 decFinalize(res, set, &residue, &status); /* .. and check */
2703 } /* finite LHS */
2704 } /* rhs OK */
2705 } /* rhs finite */
2706 if (status!=0) decStatus(res, status, set);
2707 return res;
2708 } /* decNumberScaleB */
2709
2710 /* ------------------------------------------------------------------ */
2711 /* decNumberShift -- shift the coefficient of a Number left or right */
2712 /* */
2713 /* This computes C = A << B or C = A >> -B (in base ten). */
2714 /* */
2715 /* res is C, the result. C may be A and/or B (e.g., X=X<<X) */
2716 /* lhs is A */
2717 /* rhs is B, the number of digits to shift (-ve to right) */
2718 /* set is the context */
2719 /* */
2720 /* The digits of the coefficient of A are shifted to the left (if B */
2721 /* is positive) or to the right (if B is negative) without adjusting */
2722 /* the exponent or the sign of A. */
2723 /* */
2724 /* B must be an integer (q=0) and in the range -set->digits through */
2725 /* +set->digits. */
2726 /* C must have space for set->digits digits. */
2727 /* NaNs are propagated as usual. Infinities are unaffected (but */
2728 /* B must be valid). No status is set unless B is invalid or an */
2729 /* operand is an sNaN. */
2730 /* ------------------------------------------------------------------ */
2731 decNumber * decNumberShift(decNumber *res, const decNumber *lhs,
2732 const decNumber *rhs, decContext *set) {
2733 uInt status=0; /* accumulator */
2734 Int shift; /* rhs as an Int */
2735
2736 #if DECCHECK
2737 if (decCheckOperands(res, lhs, rhs, set)) return res;
2738 #endif
2739
2740 /* NaNs propagate as normal */
2741 if (decNumberIsNaN(lhs) || decNumberIsNaN(rhs))
2742 decNaNs(res, lhs, rhs, set, &status);
2743 /* rhs must be an integer */
2744 else if (decNumberIsInfinite(rhs) || rhs->exponent!=0)
2745 status=DEC_Invalid_operation;
2746 else { /* both numeric, rhs is an integer */
2747 shift=decGetInt(rhs); /* [cannot fail] */
2748 if (shift==BADINT /* something bad .. */
2749 || shift==BIGODD || shift==BIGEVEN /* .. very big .. */
2750 || abs(shift)>set->digits) /* .. or out of range */
2751 status=DEC_Invalid_operation;
2752 else { /* rhs is OK */
2753 decNumberCopy(res, lhs);
2754 if (shift!=0 && !decNumberIsInfinite(res)) { /* something to do */
2755 if (shift>0) { /* to left */
2756 if (shift==set->digits) { /* removing all */
2757 *res->lsu=0; /* so place 0 */
2758 res->digits=1; /* .. */
2759 }
2760 else { /* */
2761 /* first remove leading digits if necessary */
2762 if (res->digits+shift>set->digits) {
2763 decDecap(res, res->digits+shift-set->digits);
2764 /* that updated res->digits; may have gone to 1 (for a */
2765 /* single digit or for zero */
2766 }
2767 if (res->digits>1 || *res->lsu) /* if non-zero.. */
2768 res->digits=decShiftToMost(res->lsu, res->digits, shift);
2769 } /* partial left */
2770 } /* left */
2771 else { /* to right */
2772 if (-shift>=res->digits) { /* discarding all */
2773 *res->lsu=0; /* so place 0 */
2774 res->digits=1; /* .. */
2775 }
2776 else {
2777 decShiftToLeast(res->lsu, D2U(res->digits), -shift);
2778 res->digits-=(-shift);
2779 }
2780 } /* to right */
2781 } /* non-0 non-Inf shift */
2782 } /* rhs OK */
2783 } /* numerics */
2784 if (status!=0) decStatus(res, status, set);
2785 return res;
2786 } /* decNumberShift */
2787
2788 /* ------------------------------------------------------------------ */
2789 /* decNumberSquareRoot -- square root operator */
2790 /* */
2791 /* This computes C = squareroot(A) */
2792 /* */
2793 /* res is C, the result. C may be A */
2794 /* rhs is A */
2795 /* set is the context; note that rounding mode has no effect */
2796 /* */
2797 /* C must have space for set->digits digits. */
2798 /* ------------------------------------------------------------------ */
2799 /* This uses the following varying-precision algorithm in: */
2800 /* */
2801 /* Properly Rounded Variable Precision Square Root, T. E. Hull and */
2802 /* A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
2803 /* pp229-237, ACM, September 1985. */
2804 /* */
2805 /* The square-root is calculated using Newton's method, after which */
2806 /* a check is made to ensure the result is correctly rounded. */
2807 /* */
2808 /* % [Reformatted original Numerical Turing source code follows.] */
2809 /* function sqrt(x : real) : real */
2810 /* % sqrt(x) returns the properly rounded approximation to the square */
2811 /* % root of x, in the precision of the calling environment, or it */
2812 /* % fails if x < 0. */
2813 /* % t e hull and a abrham, august, 1984 */
2814 /* if x <= 0 then */
2815 /* if x < 0 then */
2816 /* assert false */
2817 /* else */
2818 /* result 0 */
2819 /* end if */
2820 /* end if */
2821 /* var f := setexp(x, 0) % fraction part of x [0.1 <= x < 1] */
2822 /* var e := getexp(x) % exponent part of x */
2823 /* var approx : real */
2824 /* if e mod 2 = 0 then */
2825 /* approx := .259 + .819 * f % approx to root of f */
2826 /* else */
2827 /* f := f/l0 % adjustments */
2828 /* e := e + 1 % for odd */
2829 /* approx := .0819 + 2.59 * f % exponent */
2830 /* end if */
2831 /* */
2832 /* var p:= 3 */
2833 /* const maxp := currentprecision + 2 */
2834 /* loop */
2835 /* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */
2836 /* precision p */
2837 /* approx := .5 * (approx + f/approx) */
2838 /* exit when p = maxp */
2839 /* end loop */
2840 /* */
2841 /* % approx is now within 1 ulp of the properly rounded square root */
2842 /* % of f; to ensure proper rounding, compare squares of (approx - */
2843 /* % l/2 ulp) and (approx + l/2 ulp) with f. */
2844 /* p := currentprecision */
2845 /* begin */
2846 /* precision p + 2 */
2847 /* const approxsubhalf := approx - setexp(.5, -p) */
2848 /* if mulru(approxsubhalf, approxsubhalf) > f then */
2849 /* approx := approx - setexp(.l, -p + 1) */
2850 /* else */
2851 /* const approxaddhalf := approx + setexp(.5, -p) */
2852 /* if mulrd(approxaddhalf, approxaddhalf) < f then */
2853 /* approx := approx + setexp(.l, -p + 1) */
2854 /* end if */
2855 /* end if */
2856 /* end */
2857 /* result setexp(approx, e div 2) % fix exponent */
2858 /* end sqrt */
2859 /* ------------------------------------------------------------------ */
2860 decNumber * decNumberSquareRoot(decNumber *res, const decNumber *rhs,
2861 decContext *set) {
2862 decContext workset, approxset; /* work contexts */
2863 decNumber dzero; /* used for constant zero */
2864 Int maxp; /* largest working precision */
2865 Int workp; /* working precision */
2866 Int residue=0; /* rounding residue */
2867 uInt status=0, ignore=0; /* status accumulators */
2868 uInt rstatus; /* .. */
2869 Int exp; /* working exponent */
2870 Int ideal; /* ideal (preferred) exponent */
2871 Int needbytes; /* work */
2872 Int dropped; /* .. */
2873
2874 #if DECSUBSET
2875 decNumber *allocrhs=NULL; /* non-NULL if rounded rhs allocated */
2876 #endif
2877 /* buffer for f [needs +1 in case DECBUFFER 0] */
2878 decNumber buff[D2N(DECBUFFER+1)];
2879 /* buffer for a [needs +2 to match likely maxp] */
2880 decNumber bufa[D2N(DECBUFFER+2)];
2881 /* buffer for temporary, b [must be same size as a] */
2882 decNumber bufb[D2N(DECBUFFER+2)];
2883 decNumber *allocbuff=NULL; /* -> allocated buff, iff allocated */
2884 decNumber *allocbufa=NULL; /* -> allocated bufa, iff allocated */
2885 decNumber *allocbufb=NULL; /* -> allocated bufb, iff allocated */
2886 decNumber *f=buff; /* reduced fraction */
2887 decNumber *a=bufa; /* approximation to result */
2888 decNumber *b=bufb; /* intermediate result */
2889 /* buffer for temporary variable, up to 3 digits */
2890 decNumber buft[D2N(3)];
2891 decNumber *t=buft; /* up-to-3-digit constant or work */
2892
2893 #if DECCHECK
2894 if (decCheckOperands(res, DECUNUSED, rhs, set)) return res;
2895 #endif
2896
2897 do { /* protect allocated storage */
2898 #if DECSUBSET
2899 if (!set->extended) {
2900 /* reduce operand and set lostDigits status, as needed */
2901 if (rhs->digits>set->digits) {
2902 allocrhs=decRoundOperand(rhs, set, &status);
2903 if (allocrhs==NULL) break;
2904 /* [Note: 'f' allocation below could reuse this buffer if */
2905 /* used, but as this is rare they are kept separate for clarity.] */
2906 rhs=allocrhs;
2907 }
2908 }
2909 #endif
2910 /* [following code does not require input rounding] */
2911
2912 /* handle infinities and NaNs */
2913 if (SPECIALARG) {
2914 if (decNumberIsInfinite(rhs)) { /* an infinity */
2915 if (decNumberIsNegative(rhs)) status|=DEC_Invalid_operation;
2916 else decNumberCopy(res, rhs); /* +Infinity */
2917 }
2918 else decNaNs(res, rhs, NULL, set, &status); /* a NaN */
2919 break;
2920 }
2921
2922 /* calculate the ideal (preferred) exponent [floor(exp/2)] */
2923 /* [We would like to write: ideal=rhs->exponent>>1, but this */
2924 /* generates a compiler warning. Generated code is the same.] */
2925 ideal=(rhs->exponent&~1)/2; /* target */
2926
2927 /* handle zeros */
2928 if (ISZERO(rhs)) {
2929 decNumberCopy(res, rhs); /* could be 0 or -0 */
2930 res->exponent=ideal; /* use the ideal [safe] */
2931 /* use decFinish to clamp any out-of-range exponent, etc. */
2932 decFinish(res, set, &residue, &status);
2933 break;
2934 }
2935
2936 /* any other -x is an oops */
2937 if (decNumberIsNegative(rhs)) {
2938 status|=DEC_Invalid_operation;
2939 break;
2940 }
2941
2942 /* space is needed for three working variables */
2943 /* f -- the same precision as the RHS, reduced to 0.01->0.99... */
2944 /* a -- Hull's approximation -- precision, when assigned, is */
2945 /* currentprecision+1 or the input argument precision, */
2946 /* whichever is larger (+2 for use as temporary) */
2947 /* b -- intermediate temporary result (same size as a) */
2948 /* if any is too long for local storage, then allocate */
2949 workp=MAXI(set->digits+1, rhs->digits); /* actual rounding precision */
2950 maxp=workp+2; /* largest working precision */
2951
2952 needbytes=sizeof(decNumber)+(D2U(rhs->digits)-1)*sizeof(Unit);
2953 if (needbytes>(Int)sizeof(buff)) {
2954 allocbuff=(decNumber *)malloc(needbytes);
2955 if (allocbuff==NULL) { /* hopeless -- abandon */
2956 status|=DEC_Insufficient_storage;
2957 break;}
2958 f=allocbuff; /* use the allocated space */
2959 }
2960 /* a and b both need to be able to hold a maxp-length number */
2961 needbytes=sizeof(decNumber)+(D2U(maxp)-1)*sizeof(Unit);
2962 if (needbytes>(Int)sizeof(bufa)) { /* [same applies to b] */
2963 allocbufa=(decNumber *)malloc(needbytes);
2964 allocbufb=(decNumber *)malloc(needbytes);
2965 if (allocbufa==NULL || allocbufb==NULL) { /* hopeless */
2966 status|=DEC_Insufficient_storage;
2967 break;}
2968 a=allocbufa; /* use the allocated spaces */
2969 b=allocbufb; /* .. */
2970 }
2971
2972 /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
2973 decNumberCopy(f, rhs);
2974 exp=f->exponent+f->digits; /* adjusted to Hull rules */
2975 f->exponent=-(f->digits); /* to range */
2976
2977 /* set up working context */
2978 decContextDefault(&workset, DEC_INIT_DECIMAL64);
2979
2980 /* [Until further notice, no error is possible and status bits */
2981 /* (Rounded, etc.) should be ignored, not accumulated.] */
2982
2983 /* Calculate initial approximation, and allow for odd exponent */
2984 workset.digits=workp; /* p for initial calculation */
2985 t->bits=0; t->digits=3;
2986 a->bits=0; a->digits=3;
2987 if ((exp & 1)==0) { /* even exponent */
2988 /* Set t=0.259, a=0.819 */
2989 t->exponent=-3;
2990 a->exponent=-3;
2991 #if DECDPUN>=3
2992 t->lsu[0]=259;
2993 a->lsu[0]=819;
2994 #elif DECDPUN==2
2995 t->lsu[0]=59; t->lsu[1]=2;
2996 a->lsu[0]=19; a->lsu[1]=8;
2997 #else
2998 t->lsu[0]=9; t->lsu[1]=5; t->lsu[2]=2;
2999 a->lsu[0]=9; a->lsu[1]=1; a->lsu[2]=8;
3000 #endif
3001 }
3002 else { /* odd exponent */
3003 /* Set t=0.0819, a=2.59 */
3004 f->exponent--; /* f=f/10 */
3005 exp++; /* e=e+1 */
3006 t->exponent=-4;
3007 a->exponent=-2;
3008 #if DECDPUN>=3
3009 t->lsu[0]=819;
3010 a->lsu[0]=259;
3011 #elif DECDPUN==2
3012 t->lsu[0]=19; t->lsu[1]=8;
3013 a->lsu[0]=59; a->lsu[1]=2;
3014 #else
3015 t->lsu[0]=9; t->lsu[1]=1; t->lsu[2]=8;
3016 a->lsu[0]=9; a->lsu[1]=5; a->lsu[2]=2;
3017 #endif
3018 }
3019 decMultiplyOp(a, a, f, &workset, &ignore); /* a=a*f */
3020 decAddOp(a, a, t, &workset, 0, &ignore); /* ..+t */
3021 /* [a is now the initial approximation for sqrt(f), calculated with */
3022 /* currentprecision, which is also a's precision.] */
3023
3024 /* the main calculation loop */
3025 decNumberZero(&dzero); /* make 0 */
3026 decNumberZero(t); /* set t = 0.5 */
3027 t->lsu[0]=5; /* .. */
3028 t->exponent=-1; /* .. */
3029 workset.digits=3; /* initial p */
3030 for (;;) {
3031 /* set p to min(2*p - 2, maxp) [hence 3; or: 4, 6, 10, ... , maxp] */
3032 workset.digits=workset.digits*2-2;
3033 if (workset.digits>maxp) workset.digits=maxp;
3034 /* a = 0.5 * (a + f/a) */
3035 /* [calculated at p then rounded to currentprecision] */
3036 decDivideOp(b, f, a, &workset, DIVIDE, &ignore); /* b=f/a */
3037 decAddOp(b, b, a, &workset, 0, &ignore); /* b=b+a */
3038 decMultiplyOp(a, b, t, &workset, &ignore); /* a=b*0.5 */
3039 if (a->digits==maxp) break; /* have required digits */
3040 } /* loop */
3041
3042 /* Here, 0.1 <= a < 1 [Hull], and a has maxp digits */
3043 /* now reduce to length, etc.; this needs to be done with a */
3044 /* having the correct exponent so as to handle subnormals */
3045 /* correctly */
3046 approxset=*set; /* get emin, emax, etc. */
3047 approxset.round=DEC_ROUND_HALF_EVEN;
3048 a->exponent+=exp/2; /* set correct exponent */
3049
3050 rstatus=0; /* clear status */
3051 residue=0; /* .. and accumulator */
3052 decCopyFit(a, a, &approxset, &residue, &rstatus); /* reduce (if needed) */
3053 decFinish(a, &approxset, &residue, &rstatus); /* clean and finalize */
3054
3055 /* Overflow was possible if the input exponent was out-of-range, */
3056 /* in which case quit */
3057 if (rstatus&DEC_Overflow) {
3058 status=rstatus; /* use the status as-is */
3059 decNumberCopy(res, a); /* copy to result */
3060 break;
3061 }
3062
3063 /* Preserve status except Inexact/Rounded */
3064 status|=(rstatus & ~(DEC_Rounded|DEC_Inexact));
3065
3066 /* Carry out the Hull correction */
3067 a->exponent-=exp/2; /* back to 0.1->1 */
3068
3069 /* a is now at final precision and within 1 ulp of the properly */
3070 /* rounded square root of f; to ensure proper rounding, compare */
3071 /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
3072 /* Here workset.digits=maxp and t=0.5, and a->digits determines */
3073 /* the ulp */
3074 workset.digits--; /* maxp-1 is OK now */
3075 t->exponent=-a->digits-1; /* make 0.5 ulp */
3076 decAddOp(b, a, t, &workset, DECNEG, &ignore); /* b = a - 0.5 ulp */
3077 workset.round=DEC_ROUND_UP;
3078 decMultiplyOp(b, b, b, &workset, &ignore