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