| blob | 1ffe458ad8377e38773d0b79f7d289f84c654fbe |
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.
5 This file is part of GCC.
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.
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.)
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.
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. */
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 /* ------------------------------------------------------------------ */
174 /* Constants */
175 /* Public lookup table used by the D2U macro */
181 /* Local constants */
197 /* Next two indicate an integer >= 10**6, and its parity (bottom bit) */
198 #define BIGEVEN (Int)0x80000002
199 #define BIGODD (Int)0x80000003
203 /* Granularity-dependent code */
204 #if DECDPUN<=4
207 /* Constant multipliers for divide-by-power-of five using reciprocal */
208 /* multiply, after removing powers of 2 by shifting, and final shift */
209 /* of 17 [we only need up to **4] */
211 /* QUOT10 -- macro to return the quotient of unit u divided by 10**n */
212 #define QUOT10(u, n) ((((uInt)(u)>>(n))*multies[n])>>17)
213 #else
214 /* For DECDPUN>4 non-ANSI-89 64-bit types are needed. */
215 #if !DECUSE64
216 #error decNumber.c: DECUSE64 must be 1 when DECDPUN>4
217 #endif
220 #endif
222 /* Local routines */
246 uInt *);
251 uInt *);
267 #if !DECSUBSET
268 /* decFinish == decFinalize when no subset arithmetic needed */
269 #define decFinish(a,b,c,d) decFinalize(a,b,c,d)
270 #else
273 #endif
275 /* Local macros */
276 /* masked special-values bits */
277 #define SPECIALARG (rhs->bits & DECSPECIAL)
278 #define SPECIALARGS ((lhs->bits | rhs->bits) & DECSPECIAL)
280 /* Diagnostic macros, etc. */
281 #if DECALLOC
282 /* Handle malloc/free accounting. If enabled, our accountable routines */
283 /* are used; otherwise the code just goes straight to the system malloc */
284 /* and free routines. */
285 #define malloc(a) decMalloc(a)
286 #define free(a) decFree(a)
288 /* 'Our' malloc and free: */
292 /* Note that DECALLOC code only checks for storage buffer overflow. */
293 /* To check for memory leaks, the decAllocBytes variable must be */
294 /* checked to be 0 at appropriate times (e.g., after the test */
295 /* harness completes a set of tests). This checking may be unreliable */
296 /* if the testing is done in a multi-thread environment. */
297 #endif
299 #if DECCHECK
300 /* Optional checking routines. Enabling these means that decNumber */
301 /* and decContext operands to operator routines are checked for */
302 /* correctness. This roughly doubles the execution time of the */
303 /* fastest routines (and adds 600+ bytes), so should not normally be */
304 /* used in 'production'. */
305 /* decCheckInexact is used to check that inexact results have a full */
306 /* complement of digits (where appropriate -- this is not the case */
307 /* for Quantize, for example) */
308 #define DECUNRESU ((decNumber *)(void *)0xffffffff)
309 #define DECUNUSED ((const decNumber *)(void *)0xffffffff)
310 #define DECUNCONT ((decContext *)(void *)(0xffffffff))
315 #endif
317 #if DECTRACE || DECCHECK
318 /* Optional trace/debugging routines (may or may not be used) */
321 #endif
323 /* ================================================================== */
324 /* Conversions */
325 /* ================================================================== */
327 /* ------------------------------------------------------------------ */
328 /* from-int32 -- conversion from Int or uInt */
329 /* */
330 /* dn is the decNumber to receive the integer */
331 /* in or uin is the integer to be converted */
332 /* returns dn */
333 /* */
334 /* No error is possible. */
335 /* ------------------------------------------------------------------ */
337 uInt unsig;
342 }
343 /* in is now positive */
356 }
361 /* ------------------------------------------------------------------ */
362 /* to-int32 -- conversion to Int or uInt */
363 /* */
364 /* dn is the decNumber to convert */
365 /* set is the context for reporting errors */
366 /* returns the converted decNumber, or 0 if Invalid is set */
367 /* */
368 /* Invalid is set if the decNumber does not have exponent==0 or if */
369 /* it is a NaN, Infinite, or out-of-range. */
370 /* ------------------------------------------------------------------ */
372 #if DECCHECK
374 #endif
376 /* special or too many digits, or bad exponent */
387 #endif
388 up++;
389 /* collect remaining Units, if any, into hi */
391 /* now low has the lsd, hi the remainder */
393 /* most-negative is a reprieve */
395 /* bad -- drop through */
396 }
401 }
408 #if DECCHECK
410 #endif
411 /* special or too many digits, or bad exponent, or negative (<0) */
423 #endif
424 up++;
425 /* collect remaining Units, if any, into hi */
428 /* now low has the lsd, hi the remainder */
437 {
441 }
446 }
451 {
456 }
460 }
465 /* ------------------------------------------------------------------ */
466 /* to-int64 -- conversion to int64 */
467 /* */
468 /* dn is the decNumber to convert. dn is assumed to have been */
469 /* rounded to a floating point integer value. */
470 /* set is the context for reporting errors */
471 /* returns the converted decNumber, or 0 if Invalid is set */
472 /* */
473 /* Invalid is set if the decNumber is a NaN, Infinite or is out of */
474 /* range for a signed 64 bit integer. */
475 /* ------------------------------------------------------------------ */
478 {
493 }
494 }
500 }
502 }
504 Invalid:
510 /* ------------------------------------------------------------------ */
511 /* to-scientific-string -- conversion to numeric string */
512 /* to-engineering-string -- conversion to numeric string */
513 /* */
514 /* decNumberToString(dn, string); */
515 /* decNumberToEngString(dn, string); */
516 /* */
517 /* dn is the decNumber to convert */
518 /* string is the string where the result will be laid out */
519 /* */
520 /* string must be at least dn->digits+14 characters long */
521 /* */
522 /* No error is possible, and no status can be set. */
523 /* ------------------------------------------------------------------ */
534 /* ------------------------------------------------------------------ */
535 /* to-number -- conversion from numeric string */
536 /* */
537 /* decNumberFromString -- convert string to decNumber */
538 /* dn -- the number structure to fill */
539 /* chars[] -- the string to convert ('\0' terminated) */
540 /* set -- the context used for processing any error, */
541 /* determining the maximum precision available */
542 /* (set.digits), determining the maximum and minimum */
543 /* exponent (set.emax and set.emin), determining if */
544 /* extended values are allowed, and checking the */
545 /* rounding mode if overflow occurs or rounding is */
546 /* needed. */
547 /* */
548 /* The length of the coefficient and the size of the exponent are */
549 /* checked by this routine, so the correct error (Underflow or */
550 /* Overflow) can be reported or rounding applied, as necessary. */
551 /* */
552 /* If bad syntax is detected, the result will be a quiet NaN. */
553 /* ------------------------------------------------------------------ */
560 /* [+9 allows for ln() constants] */
568 #if DECDPUN>1
570 #endif
574 #if DECCHECK
577 #endif
585 }
592 cfirst++;
596 cfirst++;
598 }
599 /* *c is not a digit, or a valid +, -, or '.' */
606 #if DECSUBSET
607 /* if subset then infinities and NaNs are not allowed */
609 #endif
610 /* Infinities and NaNs are possible, here */
618 }
619 /* a NaN expected */
620 /* 2003.09.10 NaNs are now permitted to have a sign */
623 c++;
625 }
627 c++;
629 c++;
631 c++;
632 /* now either nothing, or nnnn payload, expected */
633 /* -> start of integer and skip leading 0s [including plain 0] */
638 }
639 /* something other than 0s; setup last and d as usual [no dots] */
643 }
646 /* [NB: payload in a decNumber can be full length unless */
647 /* clamped, in which case can only be digits-1] */
651 /* good; drop through to convert the integer to coefficient */
657 /* had some digits; exponent is only valid sequence now */
662 /* Found 'e' or 'E' -- now process explicit exponent */
663 /* 1998.07.11: sign no longer required */
676 /* if not now on a '\0', *c must not be a digit */
679 /* (this next test must be after the syntax checks) */
680 /* if it was too long the exponent may have wrapped, so check */
681 /* carefully and set it to a certain overflow if wrap possible */
684 /* [up to 1999999999 is OK, for example 1E-1000000998] */
685 }
690 /* Here when whole string has been inspected; syntax is good */
691 /* cfirst->first digit (never dot), last->last digit (ditto) */
693 /* strip leading zeros/dot [leave final 0 if all 0's] */
700 #if DECSUBSET
701 /* make a rapid exit for easy zeros if !extended */
705 }
706 #endif
709 /* Handle decimal point... */
712 /* [we can now ignore the .] */
714 /* OK, the digits string is good. Assemble in the decNumber, or in */
715 /* a temporary units array if rounding is needed */
724 }
725 }
726 /* res now -> number lsu, buffer, or allocated storage for Unit array */
728 /* Place the coefficient into the selected Unit array */
729 /* [this is often 70% of the cost of this function when DECDPUN>1] */
730 #if DECDPUN>1
738 cut--;
747 #else
748 /* DECDPUN==1 */
753 up++;
755 #endif
761 /* if not in number (too long) shorten into the number */
765 /* always check for overflow or subnormal and round as needed */
767 }
769 /* [these tests are just for performance; finalize repeats them] */
774 }
775 }
776 /* decNumberShow(dn); */
784 /* ================================================================== */
785 /* Operators */
786 /* ================================================================== */
788 /* ------------------------------------------------------------------ */
789 /* decNumberAbs -- absolute value operator */
790 /* */
791 /* This computes C = abs(A) */
792 /* */
793 /* res is C, the result. C may be A */
794 /* rhs is A */
795 /* set is the context */
796 /* */
797 /* See also decNumberCopyAbs for a quiet bitwise version of this. */
798 /* C must have space for set->digits digits. */
799 /* ------------------------------------------------------------------ */
800 /* This has the same effect as decNumberPlus unless A is negative, */
801 /* in which case it has the same effect as decNumberMinus. */
802 /* ------------------------------------------------------------------ */
808 #if DECCHECK
810 #endif
816 #if DECCHECK
818 #endif
822 /* ------------------------------------------------------------------ */
823 /* decNumberAdd -- add two Numbers */
824 /* */
825 /* This computes C = A + B */
826 /* */
827 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
828 /* lhs is A */
829 /* rhs is B */
830 /* set is the context */
831 /* */
832 /* C must have space for set->digits digits. */
833 /* ------------------------------------------------------------------ */
834 /* This just calls the routine shared with Subtract */
840 #if DECCHECK
842 #endif
846 /* ------------------------------------------------------------------ */
847 /* decNumberAnd -- AND two Numbers, digitwise */
848 /* */
849 /* This computes C = A & B */
850 /* */
851 /* res is C, the result. C may be A and/or B (e.g., X=X&X) */
852 /* lhs is A */
853 /* rhs is B */
854 /* set is the context (used for result length and error report) */
855 /* */
856 /* C must have space for set->digits digits. */
857 /* */
858 /* Logical function restrictions apply (see above); a NaN is */
859 /* returned with Invalid_operation if a restriction is violated. */
860 /* ------------------------------------------------------------------ */
867 #if DECCHECK
869 #endif
875 }
877 /* operands are valid */
895 /* This loop could be unrolled and/or use BIN2BCD tables */
905 }
910 /* [here uc-1 is the msu of the result] */
917 /* ------------------------------------------------------------------ */
918 /* decNumberCompare -- compare two Numbers */
919 /* */
920 /* This computes C = A ? B */
921 /* */
922 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
923 /* lhs is A */
924 /* rhs is B */
925 /* set is the context */
926 /* */
927 /* C must have space for one digit (or NaN). */
928 /* ------------------------------------------------------------------ */
937 /* ------------------------------------------------------------------ */
938 /* decNumberCompareSignal -- compare, signalling on all NaNs */
939 /* */
940 /* This computes C = A ? B */
941 /* */
942 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
943 /* lhs is A */
944 /* rhs is B */
945 /* set is the context */
946 /* */
947 /* C must have space for one digit (or NaN). */
948 /* ------------------------------------------------------------------ */
957 /* ------------------------------------------------------------------ */
958 /* decNumberCompareTotal -- compare two Numbers, using total ordering */
959 /* */
960 /* This computes C = A ? B, under total ordering */
961 /* */
962 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
963 /* lhs is A */
964 /* rhs is B */
965 /* set is the context */
966 /* */
967 /* C must have space for one digit; the result will always be one of */
968 /* -1, 0, or 1. */
969 /* ------------------------------------------------------------------ */
978 /* ------------------------------------------------------------------ */
979 /* decNumberCompareTotalMag -- compare, total ordering of magnitudes */
980 /* */
981 /* This computes C = |A| ? |B|, under total ordering */
982 /* */
983 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
984 /* lhs is A */
985 /* rhs is B */
986 /* set is the context */
987 /* */
988 /* C must have space for one digit; the result will always be one of */
989 /* -1, 0, or 1. */
990 /* ------------------------------------------------------------------ */
1001 #if DECCHECK
1003 #endif
1006 /* if either is negative, take a copy and absolute */
1016 }
1020 }
1030 }
1034 }
1044 /* ------------------------------------------------------------------ */
1045 /* decNumberDivide -- divide one number by another */
1046 /* */
1047 /* This computes C = A / B */
1048 /* */
1049 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
1050 /* lhs is A */
1051 /* rhs is B */
1052 /* set is the context */
1053 /* */
1054 /* C must have space for set->digits digits. */
1055 /* ------------------------------------------------------------------ */
1061 #if DECCHECK
1063 #endif
1067 /* ------------------------------------------------------------------ */
1068 /* decNumberDivideInteger -- divide and return integer quotient */
1069 /* */
1070 /* This computes C = A # B, where # is the integer divide operator */
1071 /* */
1072 /* res is C, the result. C may be A and/or B (e.g., X=X#X) */
1073 /* lhs is A */
1074 /* rhs is B */
1075 /* set is the context */
1076 /* */
1077 /* C must have space for set->digits digits. */
1078 /* ------------------------------------------------------------------ */
1087 /* ------------------------------------------------------------------ */
1088 /* decNumberExp -- exponentiation */
1089 /* */
1090 /* This computes C = exp(A) */
1091 /* */
1092 /* res is C, the result. C may be A */
1093 /* rhs is A */
1094 /* set is the context; note that rounding mode has no effect */
1095 /* */
1096 /* C must have space for set->digits digits. */
1097 /* */
1098 /* Mathematical function restrictions apply (see above); a NaN is */
1099 /* returned with Invalid_operation if a restriction is violated. */
1100 /* */
1101 /* Finite results will always be full precision and Inexact, except */
1102 /* when A is a zero or -Infinity (giving 1 or 0 respectively). */
1103 /* */
1104 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
1105 /* almost always be correctly rounded, but may be up to 1 ulp in */
1106 /* error in rare cases. */
1107 /* ------------------------------------------------------------------ */
1108 /* This is a wrapper for decExpOp which can handle the slightly wider */
1109 /* (double) range needed by Ln (which has to be able to calculate */
1110 /* exp(-a) where a can be the tiniest number (Ntiny). */
1111 /* ------------------------------------------------------------------ */
1115 #if DECSUBSET
1117 #endif
1119 #if DECCHECK
1121 #endif
1123 /* Check restrictions; these restrictions ensure that if h=8 (see */
1124 /* decExpOp) then the result will either overflow or underflow to 0. */
1125 /* Other math functions restrict the input range, too, for inverses. */
1126 /* If not violated then carry out the operation. */
1128 #if DECSUBSET
1130 /* reduce operand and set lostDigits status, as needed */
1135 }
1136 }
1137 #endif
1141 #if DECSUBSET
1143 #endif
1144 /* apply significant status */
1146 #if DECCHECK
1148 #endif
1152 /* ------------------------------------------------------------------ */
1153 /* decNumberFMA -- fused multiply add */
1154 /* */
1155 /* This computes D = (A * B) + C with only one rounding */
1156 /* */
1157 /* res is D, the result. D may be A or B or C (e.g., X=FMA(X,X,X)) */
1158 /* lhs is A */
1159 /* rhs is B */
1160 /* fhs is C [far hand side] */
1161 /* set is the context */
1162 /* */
1163 /* Mathematical function restrictions apply (see above); a NaN is */
1164 /* returned with Invalid_operation if a restriction is violated. */
1165 /* */
1166 /* C must have space for set->digits digits. */
1167 /* ------------------------------------------------------------------ */
1179 #if DECCHECK
1182 #endif
1185 #if DECSUBSET
1189 #endif
1190 /* Check math restrictions [these ensure no overflow or underflow] */
1194 /* set up context for multiply */
1197 /* [The above may be an over-estimate for subset arithmetic, but that's OK] */
1200 /* set up decNumber space to receive the result of the multiply */
1209 }
1210 /* multiply with extended range and necessary precision */
1211 /*printf("emin=%ld\n", dcmul.emin); */
1213 /* Only Invalid operation (from sNaN or Inf * 0) is possible in */
1214 /* status; if either is seen than ignore fhs (in case it is */
1215 /* another sNaN) and set acc to NaN unless we had an sNaN */
1216 /* [decMultiplyOp leaves that to caller] */
1217 /* Note sNaN has to go through addOp to shorten payload if */
1218 /* necessary */
1224 }
1227 }
1228 #if DECCHECK
1231 }
1232 #endif
1233 /* add the third operand and result -> res, and all is done */
1239 #if DECCHECK
1241 #endif
1245 /* ------------------------------------------------------------------ */
1246 /* decNumberInvert -- invert a Number, digitwise */
1247 /* */
1248 /* This computes C = ~A */
1249 /* */
1250 /* res is C, the result. C may be A (e.g., X=~X) */
1251 /* rhs is A */
1252 /* set is the context (used for result length and error report) */
1253 /* */
1254 /* C must have space for set->digits digits. */
1255 /* */
1256 /* Logical function restrictions apply (see above); a NaN is */
1257 /* returned with Invalid_operation if a restriction is violated. */
1258 /* ------------------------------------------------------------------ */
1264 #if DECCHECK
1266 #endif
1271 }
1272 /* operand is valid */
1284 /* always need to examine all bits in rhs */
1285 /* This loop could be unrolled and/or use BIN2BCD tables */
1293 }
1297 /* [here uc-1 is the msu of the result] */
1304 /* ------------------------------------------------------------------ */
1305 /* decNumberLn -- natural logarithm */
1306 /* */
1307 /* This computes C = ln(A) */
1308 /* */
1309 /* res is C, the result. C may be A */
1310 /* rhs is A */
1311 /* set is the context; note that rounding mode has no effect */
1312 /* */
1313 /* C must have space for set->digits digits. */
1314 /* */
1315 /* Notable cases: */
1316 /* A<0 -> Invalid */
1317 /* A=0 -> -Infinity (Exact) */
1318 /* A=+Infinity -> +Infinity (Exact) */
1319 /* A=1 exactly -> 0 (Exact) */
1320 /* */
1321 /* Mathematical function restrictions apply (see above); a NaN is */
1322 /* returned with Invalid_operation if a restriction is violated. */
1323 /* */
1324 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
1325 /* almost always be correctly rounded, but may be up to 1 ulp in */
1326 /* error in rare cases. */
1327 /* ------------------------------------------------------------------ */
1328 /* This is a wrapper for decLnOp which can handle the slightly wider */
1329 /* (+11) range needed by Ln, Log10, etc. (which may have to be able */
1330 /* to calculate at p+e+2). */
1331 /* ------------------------------------------------------------------ */
1335 #if DECSUBSET
1337 #endif
1339 #if DECCHECK
1341 #endif
1343 /* Check restrictions; this is a math function; if not violated */
1344 /* then carry out the operation. */
1346 #if DECSUBSET
1348 /* reduce operand and set lostDigits status, as needed */
1353 }
1354 /* special check in subset for rhs=0 */
1359 #endif
1363 #if DECSUBSET
1365 #endif
1366 /* apply significant status */
1368 #if DECCHECK
1370 #endif
1374 /* ------------------------------------------------------------------ */
1375 /* decNumberLogB - get adjusted exponent, by 754r rules */
1376 /* */
1377 /* This computes C = adjustedexponent(A) */
1378 /* */
1379 /* res is C, the result. C may be A */
1380 /* rhs is A */
1381 /* set is the context, used only for digits and status */
1382 /* */
1383 /* C must have space for 10 digits (A might have 10**9 digits and */
1384 /* an exponent of +999999999, or one digit and an exponent of */
1385 /* -1999999999). */
1386 /* */
1387 /* This returns the adjusted exponent of A after (in theory) padding */
1388 /* with zeros on the right to set->digits digits while keeping the */
1389 /* same value. The exponent is not limited by emin/emax. */
1390 /* */
1391 /* Notable cases: */
1392 /* A<0 -> Use |A| */
1393 /* A=0 -> -Infinity (Division by zero) */
1394 /* A=Infinite -> +Infinity (Exact) */
1395 /* A=1 exactly -> 0 (Exact) */
1396 /* NaNs are propagated as usual */
1397 /* ------------------------------------------------------------------ */
1402 #if DECCHECK
1404 #endif
1406 /* NaNs as usual; Infinities return +Infinity; 0->oops */
1413 }
1417 }
1423 /* ------------------------------------------------------------------ */
1424 /* decNumberLog10 -- logarithm in base 10 */
1425 /* */
1426 /* This computes C = log10(A) */
1427 /* */
1428 /* res is C, the result. C may be A */
1429 /* rhs is A */
1430 /* set is the context; note that rounding mode has no effect */
1431 /* */
1432 /* C must have space for set->digits digits. */
1433 /* */
1434 /* Notable cases: */
1435 /* A<0 -> Invalid */
1436 /* A=0 -> -Infinity (Exact) */
1437 /* A=+Infinity -> +Infinity (Exact) */
1438 /* A=10**n (if n is an integer) -> n (Exact) */
1439 /* */
1440 /* Mathematical function restrictions apply (see above); a NaN is */
1441 /* returned with Invalid_operation if a restriction is violated. */
1442 /* */
1443 /* An Inexact result is rounded using DEC_ROUND_HALF_EVEN; it will */
1444 /* almost always be correctly rounded, but may be up to 1 ulp in */
1445 /* error in rare cases. */
1446 /* ------------------------------------------------------------------ */
1447 /* This calculates ln(A)/ln(10) using appropriate precision. For */
1448 /* ln(A) this is the max(p, rhs->digits + t) + 3, where p is the */
1449 /* requested digits and t is the number of digits in the exponent */
1450 /* (maximum 6). For ln(10) it is p + 3; this is often handled by the */
1451 /* fastpath in decLnOp. The final division is done to the requested */
1452 /* precision. */
1453 /* ------------------------------------------------------------------ */
1461 /* buffers for a and b working decimals */
1462 /* (adjustment calculator, same size) */
1471 #if DECSUBSET
1473 #endif
1477 #if DECCHECK
1479 #endif
1481 /* Check restrictions; this is a math function; if not violated */
1482 /* then carry out the operation. */
1484 #if DECSUBSET
1486 /* reduce operand and set lostDigits status, as needed */
1491 }
1492 /* special check in subset for rhs=0 */
1497 #endif
1501 /* handle exact powers of 10; only check if +ve finite */
1506 /* round to a single digit... */
1509 /* if exact and the digit is 1, rhs is a power of 10 */
1511 /* the exponent, conveniently, is the power of 10; making */
1512 /* this the result needs a little care as it might not fit, */
1513 /* so first convert it into the working number, and then move */
1514 /* to res */
1523 /* simplify the information-content calculation to use 'total */
1524 /* number of digits in a, including exponent' as compared to the */
1525 /* requested digits, as increasing this will only rarely cost an */
1526 /* iteration in ln(a) anyway */
1529 /* allocate space when needed... */
1538 }
1545 /* skip the division if the result so far is infinite, NaN, or */
1546 /* zero, or there was an error; note NaN from sNaN needs copy */
1552 /* for ln(10) an extra 3 digits of precision are needed */
1561 }
1563 #if DECDPUN==1
1565 #else
1567 #endif
1579 #if DECSUBSET
1581 #endif
1582 /* apply significant status */
1584 #if DECCHECK
1586 #endif
1590 /* ------------------------------------------------------------------ */
1591 /* decNumberMax -- compare two Numbers and return the maximum */
1592 /* */
1593 /* This computes C = A ? B, returning the maximum by 754R rules */
1594 /* */
1595 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1596 /* lhs is A */
1597 /* rhs is B */
1598 /* set is the context */
1599 /* */
1600 /* C must have space for set->digits digits. */
1601 /* ------------------------------------------------------------------ */
1607 #if DECCHECK
1609 #endif
1613 /* ------------------------------------------------------------------ */
1614 /* decNumberMaxMag -- compare and return the maximum by magnitude */
1615 /* */
1616 /* This computes C = A ? B, returning the maximum by 754R rules */
1617 /* */
1618 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1619 /* lhs is A */
1620 /* rhs is B */
1621 /* set is the context */
1622 /* */
1623 /* C must have space for set->digits digits. */
1624 /* ------------------------------------------------------------------ */
1630 #if DECCHECK
1632 #endif
1636 /* ------------------------------------------------------------------ */
1637 /* decNumberMin -- compare two Numbers and return the minimum */
1638 /* */
1639 /* This computes C = A ? B, returning the minimum by 754R rules */
1640 /* */
1641 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1642 /* lhs is A */
1643 /* rhs is B */
1644 /* set is the context */
1645 /* */
1646 /* C must have space for set->digits digits. */
1647 /* ------------------------------------------------------------------ */
1653 #if DECCHECK
1655 #endif
1659 /* ------------------------------------------------------------------ */
1660 /* decNumberMinMag -- compare and return the minimum by magnitude */
1661 /* */
1662 /* This computes C = A ? B, returning the minimum by 754R rules */
1663 /* */
1664 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
1665 /* lhs is A */
1666 /* rhs is B */
1667 /* set is the context */
1668 /* */
1669 /* C must have space for set->digits digits. */
1670 /* ------------------------------------------------------------------ */
1676 #if DECCHECK
1678 #endif
1682 /* ------------------------------------------------------------------ */
1683 /* decNumberMinus -- prefix minus operator */
1684 /* */
1685 /* This computes C = 0 - A */
1686 /* */
1687 /* res is C, the result. C may be A */
1688 /* rhs is A */
1689 /* set is the context */
1690 /* */
1691 /* See also decNumberCopyNegate for a quiet bitwise version of this. */
1692 /* C must have space for set->digits digits. */
1693 /* ------------------------------------------------------------------ */
1694 /* Simply use AddOp for the subtract, which will do the necessary. */
1695 /* ------------------------------------------------------------------ */
1698 decNumber dzero;
1701 #if DECCHECK
1703 #endif
1709 #if DECCHECK
1711 #endif
1715 /* ------------------------------------------------------------------ */
1716 /* decNumberNextMinus -- next towards -Infinity */
1717 /* */
1718 /* This computes C = A - infinitesimal, rounded towards -Infinity */
1719 /* */
1720 /* res is C, the result. C may be A */
1721 /* rhs is A */
1722 /* set is the context */
1723 /* */
1724 /* This is a generalization of 754r NextDown. */
1725 /* ------------------------------------------------------------------ */
1731 #if DECCHECK
1733 #endif
1735 /* +Infinity is the special case */
1738 /* there is no status to set */
1740 }
1751 /* ------------------------------------------------------------------ */
1752 /* decNumberNextPlus -- next towards +Infinity */
1753 /* */
1754 /* This computes C = A + infinitesimal, rounded towards +Infinity */
1755 /* */
1756 /* res is C, the result. C may be A */
1757 /* rhs is A */
1758 /* set is the context */
1759 /* */
1760 /* This is a generalization of 754r NextUp. */
1761 /* ------------------------------------------------------------------ */
1767 #if DECCHECK
1769 #endif
1771 /* -Infinity is the special case */
1775 /* there is no status to set */
1777 }
1788 /* ------------------------------------------------------------------ */
1789 /* decNumberNextToward -- next towards rhs */
1790 /* */
1791 /* This computes C = A +/- infinitesimal, rounded towards */
1792 /* +/-Infinity in the direction of B, as per 754r nextafter rules */
1793 /* */
1794 /* res is C, the result. C may be A or B. */
1795 /* lhs is A */
1796 /* rhs is B */
1797 /* set is the context */
1798 /* */
1799 /* This is a generalization of 754r NextAfter. */
1800 /* ------------------------------------------------------------------ */
1807 #if DECCHECK
1809 #endif
1813 }
1822 /* -Infinity is the special case */
1827 }
1832 /* +Infinity is the special case */
1836 }
1844 /* turn off exceptions if the result is a normal number */
1845 /* (including Nmin), otherwise let all status through */
1854 /* ------------------------------------------------------------------ */
1855 /* decNumberOr -- OR two Numbers, digitwise */
1856 /* */
1857 /* This computes C = A | B */
1858 /* */
1859 /* res is C, the result. C may be A and/or B (e.g., X=X|X) */
1860 /* lhs is A */
1861 /* rhs is B */
1862 /* set is the context (used for result length and error report) */
1863 /* */
1864 /* C must have space for set->digits digits. */
1865 /* */
1866 /* Logical function restrictions apply (see above); a NaN is */
1867 /* returned with Invalid_operation if a restriction is violated. */
1868 /* ------------------------------------------------------------------ */
1875 #if DECCHECK
1877 #endif
1883 }
1884 /* operands are valid */
1901 /* This loop could be unrolled and/or use BIN2BCD tables */
1911 }
1916 /* [here uc-1 is the msu of the result] */
1923 /* ------------------------------------------------------------------ */
1924 /* decNumberPlus -- prefix plus operator */
1925 /* */
1926 /* This computes C = 0 + A */
1927 /* */
1928 /* res is C, the result. C may be A */
1929 /* rhs is A */
1930 /* set is the context */
1931 /* */
1932 /* See also decNumberCopy for a quiet bitwise version of this. */
1933 /* C must have space for set->digits digits. */
1934 /* ------------------------------------------------------------------ */
1935 /* This simply uses AddOp; Add will take fast path after preparing A. */
1936 /* Performance is a concern here, as this routine is often used to */
1937 /* check operands and apply rounding and overflow/underflow testing. */
1938 /* ------------------------------------------------------------------ */
1941 decNumber dzero;
1943 #if DECCHECK
1945 #endif
1951 #if DECCHECK
1953 #endif
1957 /* ------------------------------------------------------------------ */
1958 /* decNumberMultiply -- multiply two Numbers */
1959 /* */
1960 /* This computes C = A x B */
1961 /* */
1962 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
1963 /* lhs is A */
1964 /* rhs is B */
1965 /* set is the context */
1966 /* */
1967 /* C must have space for set->digits digits. */
1968 /* ------------------------------------------------------------------ */
1974 #if DECCHECK
1976 #endif
1980 /* ------------------------------------------------------------------ */
1981 /* decNumberPower -- raise a number to a power */
1982 /* */
1983 /* This computes C = A ** B */
1984 /* */
1985 /* res is C, the result. C may be A and/or B (e.g., X=X**X) */
1986 /* lhs is A */
1987 /* rhs is B */
1988 /* set is the context */
1989 /* */
1990 /* C must have space for set->digits digits. */
1991 /* */
1992 /* Mathematical function restrictions apply (see above); a NaN is */
1993 /* returned with Invalid_operation if a restriction is violated. */
1994 /* */
1995 /* However, if 1999999997<=B<=999999999 and B is an integer then the */
1996 /* restrictions on A and the context are relaxed to the usual bounds, */
1997 /* for compatibility with the earlier (integer power only) version */
1998 /* of this function. */
1999 /* */
2000 /* When B is an integer, the result may be exact, even if rounded. */
2001 /* */
2002 /* The final result is rounded according to the context; it will */
2003 /* almost always be correctly rounded, but may be up to 1 ulp in */
2004 /* error in rare cases. */
2005 /* ------------------------------------------------------------------ */
2008 #if DECSUBSET
2011 #endif
2020 #if DECSUBSET
2022 #endif
2030 /* local accumulator buffer [a decNumber, with digits+elength+1 digits] */
2033 /* same again for possible 1/lhs calculation */
2036 #if DECCHECK
2038 #endif
2041 #if DECSUBSET
2047 }
2052 }
2053 }
2054 #endif
2055 /* [following code does not require input rounding] */
2057 /* handle NaNs and rhs Infinity (lhs infinity is harder) */
2074 }
2076 /* 1**Infinity is inexact, so return fully-padded 1.0000 */
2082 }
2085 }
2088 /* [lhs infinity drops through] */
2091 /* Original rhs may be an integer that fits and is in range */
2098 }
2103 /* handle LHS infinity */
2109 /* -Inf**nonint -> error */
2114 /* [otherwise will be 0 or -0] */
2116 }
2119 /* similarly handle LHS zero */
2122 #if DECSUBSET
2127 #endif
2129 }
2133 #if DECSUBSET
2137 #endif
2139 }
2141 /* [otherwise will be 0 or -0] */
2143 }
2146 /* here both lhs and rhs are finite; rhs==0 is handled in the */
2147 /* integer path. Next handle the non-integer cases */
2149 /* any -ve lhs is bad, as is either operand or context out of */
2150 /* bounds */
2162 /* calculate the result using exp(ln(lhs)*rhs), which can */
2163 /* all be done into the accumulator, dac. The precision needed */
2164 /* is enough to contain the full information in the lhs (which */
2165 /* is the total digits, including exponent), or the requested */
2166 /* precision, if larger, + 4; 6 is used for the exponent */
2167 /* maximum length, and this is also used when it is shorter */
2168 /* than the requested digits as it greatly reduces the >0.5 ulp */
2169 /* cases at little cost (because Ln doubles digits each */
2170 /* iteration so a few extra digits rarely causes an extra */
2171 /* iteration) */
2177 /* (0**0 was handled above) */
2181 /* rhs is a non-zero integer */
2186 /* calculate the working DIGITS */
2188 #if DECSUBSET
2190 #endif
2191 /* it's an error if this is more than can be handled */
2195 /* aset.digits is the count of digits for the accumulator needed */
2196 /* if accumulator is too long for local storage, then allocate */
2198 /* [needbytes also used below if 1/lhs needed] */
2205 }
2206 /* here, aset is set up and accumulator is ready for use */
2209 /* x ** y; special-case x=1 here as it will otherwise always */
2210 /* reduce to integer 1; decLnOp has a fastpath which detects */
2211 /* the case of x=1 */
2213 /* [no error possible, as lhs 0 already handled] */
2215 /* need to return fully-padded 1.0000 etc., but rhsint->1 */
2222 }
2223 }
2227 }
2228 /* and drop through for final rounding */
2235 /* if a negative power the constant 1 is needed, and if not subset */
2236 /* invert the lhs now rather than inverting the result later */
2240 #if DECSUBSET
2242 #endif
2243 /* divide lhs into 1, putting result in dac [dac=1/dac] */
2245 /* now locate or allocate space for the inverted lhs */
2252 }
2253 /* [inv now points to big-enough buffer or allocated storage] */
2257 #if DECSUBSET
2258 }
2259 #endif
2260 }
2262 /* Raise-to-the-power loop... */
2265 /* abandon if had overflow or terminal underflow */
2268 }
2269 /* [the following two lines revealed an optimizer bug in a C++ */
2270 /* compiler, with symptom: 5**3 -> 25, when n=n+n was used] */
2275 }
2281 /* complete internal overflow or underflow processing */
2283 #if DECSUBSET
2284 /* If subset, and power was negative, reverse the kind of -erflow */
2285 /* [1/x not yet done] */
2292 }
2293 }
2294 #endif
2296 /* round subnormals [to set.digits rather than aset.digits] */
2297 /* or set overflow result similarly as required */
2301 }
2303 #if DECSUBSET
2306 /* so divide result into 1 [dac=1/dac] */
2308 }
2309 #endif
2312 /* reduce result to the requested length and copy to result */
2315 #if DECSUBSET
2317 #endif
2322 #if DECSUBSET
2325 #endif
2327 #if DECCHECK
2329 #endif
2333 /* ------------------------------------------------------------------ */
2334 /* decNumberQuantize -- force exponent to requested value */
2335 /* */
2336 /* This computes C = op(A, B), where op adjusts the coefficient */
2337 /* of C (by rounding or shifting) such that the exponent (-scale) */
2338 /* of C has exponent of B. The numerical value of C will equal A, */
2339 /* except for the effects of any rounding that occurred. */
2340 /* */
2341 /* res is C, the result. C may be A or B */
2342 /* lhs is A, the number to adjust */
2343 /* rhs is B, the number with exponent to match */
2344 /* set is the context */
2345 /* */
2346 /* C must have space for set->digits digits. */
2347 /* */
2348 /* Unless there is an error or the result is infinite, the exponent */
2349 /* after the operation is guaranteed to be equal to that of B. */
2350 /* ------------------------------------------------------------------ */
2359 /* ------------------------------------------------------------------ */
2360 /* decNumberReduce -- remove trailing zeros */
2361 /* */
2362 /* This computes C = 0 + A, and normalizes the result */
2363 /* */
2364 /* res is C, the result. C may be A */
2365 /* rhs is A */
2366 /* set is the context */
2367 /* */
2368 /* C must have space for set->digits digits. */
2369 /* ------------------------------------------------------------------ */
2370 /* Previously known as Normalize */
2378 #if DECSUBSET
2380 #endif
2385 #if DECCHECK
2387 #endif
2390 #if DECSUBSET
2392 /* reduce operand and set lostDigits status, as needed */
2397 }
2398 }
2399 #endif
2400 /* [following code does not require input rounding] */
2402 /* Infinities copy through; NaNs need usual treatment */
2406 }
2408 /* reduce result to the requested length and copy to result */
2414 #if DECSUBSET
2416 #endif
2421 /* ------------------------------------------------------------------ */
2422 /* decNumberRescale -- force exponent to requested value */
2423 /* */
2424 /* This computes C = op(A, B), where op adjusts the coefficient */
2425 /* of C (by rounding or shifting) such that the exponent (-scale) */
2426 /* of C has the value B. The numerical value of C will equal A, */
2427 /* except for the effects of any rounding that occurred. */
2428 /* */
2429 /* res is C, the result. C may be A or B */
2430 /* lhs is A, the number to adjust */
2431 /* rhs is B, the requested exponent */
2432 /* set is the context */
2433 /* */
2434 /* C must have space for set->digits digits. */
2435 /* */
2436 /* Unless there is an error or the result is infinite, the exponent */
2437 /* after the operation is guaranteed to be equal to B. */
2438 /* ------------------------------------------------------------------ */
2447 /* ------------------------------------------------------------------ */
2448 /* decNumberRemainder -- divide and return remainder */
2449 /* */
2450 /* This computes C = A % B */
2451 /* */
2452 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
2453 /* lhs is A */
2454 /* rhs is B */
2455 /* set is the context */
2456 /* */
2457 /* C must have space for set->digits digits. */
2458 /* ------------------------------------------------------------------ */
2464 #if DECCHECK
2466 #endif
2470 /* ------------------------------------------------------------------ */
2471 /* decNumberRemainderNear -- divide and return remainder from nearest */
2472 /* */
2473 /* This computes C = A % B, where % is the IEEE remainder operator */
2474 /* */
2475 /* res is C, the result. C may be A and/or B (e.g., X=X%X) */
2476 /* lhs is A */
2477 /* rhs is B */
2478 /* set is the context */
2479 /* */
2480 /* C must have space for set->digits digits. */
2481 /* ------------------------------------------------------------------ */
2487 #if DECCHECK
2489 #endif
2493 /* ------------------------------------------------------------------ */
2494 /* decNumberRotate -- rotate the coefficient of a Number left/right */
2495 /* */
2496 /* This computes C = A rot B (in base ten and rotating set->digits */
2497 /* digits). */
2498 /* */
2499 /* res is C, the result. C may be A and/or B (e.g., X=XrotX) */
2500 /* lhs is A */
2501 /* rhs is B, the number of digits to rotate (-ve to right) */
2502 /* set is the context */
2503 /* */
2504 /* The digits of the coefficient of A are rotated to the left (if B */
2505 /* is positive) or to the right (if B is negative) without adjusting */
2506 /* the exponent or the sign of A. If lhs->digits is less than */
2507 /* set->digits the coefficient is padded with zeros on the left */
2508 /* before the rotate. Any leading zeros in the result are removed */
2509 /* as usual. */
2510 /* */
2511 /* B must be an integer (q=0) and in the range -set->digits through */
2512 /* +set->digits. */
2513 /* C must have space for set->digits digits. */
2514 /* NaNs are propagated as usual. Infinities are unaffected (but */
2515 /* B must be valid). No status is set unless B is invalid or an */
2516 /* operand is an sNaN. */
2517 /* ------------------------------------------------------------------ */
2523 #if DECCHECK
2525 #endif
2527 /* NaNs propagate as normal */
2530 /* rhs must be an integer */
2541 /* convert -ve rotate to equivalent positive rotation */
2545 /* left-rotate to do; 0 < rotate < set->digits */
2554 /* rotation here is done in-place, in three steps */
2555 /* 1. shift all to least up to one unit to unit-align final */
2556 /* lsd [any digits shifted out are rotated to the left, */
2557 /* abutted to the original msd (which may require split)] */
2558 /* */
2559 /* [if there are no whole units left to rotate, the */
2560 /* rotation is now complete] */
2561 /* */
2562 /* 2. shift to least, from below the split point only, so that */
2563 /* the final msd is in the right place in its Unit [any */
2564 /* digits shifted out will fit exactly in the current msu, */
2565 /* left aligned, no split required] */
2566 /* */
2567 /* 3. rotate all the units by reversing left part, right */
2568 /* part, and then whole */
2569 /* */
2570 /* example: rotate right 8 digits (2 units + 2), DECDPUN=3. */
2571 /* */
2572 /* start: 00a bcd efg hij klm npq */
2573 /* */
2574 /* 1a 000 0ab cde fgh|ijk lmn [pq saved] */
2575 /* 1b 00p qab cde fgh|ijk lmn */
2576 /* */
2577 /* 2a 00p qab cde fgh|00i jkl [mn saved] */
2578 /* 2b mnp qab cde fgh|00i jkl */
2579 /* */
2580 /* 3a fgh cde qab mnp|00i jkl */
2581 /* 3b fgh cde qab mnp|jkl 00i */
2582 /* 3c 00i jkl mnp qab cde fgh */
2584 /* Step 1: amount to shift is the partial right-rotate count */
2596 }
2599 }
2602 /* If whole units to rotate... */
2604 /* Step 2: the units to touch are the whole ones in rotate, */
2605 /* if any, and the shift is DECDPUN-msudigits (which may be */
2606 /* 0, again) */
2614 /* Step 3: rotate the units array using triple reverse */
2615 /* (reversing is easy and fast) */
2620 /* the rotation may have left an undetermined number of zeros */
2621 /* on the left, so true length needs to be calculated */
2630 /* ------------------------------------------------------------------ */
2631 /* decNumberSameQuantum -- test for equal exponents */
2632 /* */
2633 /* res is the result number, which will contain either 0 or 1 */
2634 /* lhs is a number to test */
2635 /* rhs is the second (usually a pattern) */
2636 /* */
2637 /* No errors are possible and no context is needed. */
2638 /* ------------------------------------------------------------------ */
2643 #if DECCHECK
2645 #endif
2650 /* [anything else with a special gives 0] */
2651 }
2659 /* ------------------------------------------------------------------ */
2660 /* decNumberScaleB -- multiply by a power of 10 */
2661 /* */
2662 /* This computes C = A x 10**B where B is an integer (q=0) with */
2663 /* maximum magnitude 2*(emax+digits) */
2664 /* */
2665 /* res is C, the result. C may be A or B */
2666 /* lhs is A, the number to adjust */
2667 /* rhs is B, the requested power of ten to use */
2668 /* set is the context */
2669 /* */
2670 /* C must have space for set->digits digits. */
2671 /* */
2672 /* The result may underflow or overflow. */
2673 /* ------------------------------------------------------------------ */
2680 #if DECCHECK
2682 #endif
2684 /* Handle special values except lhs infinite */
2687 /* rhs must be an integer */
2691 /* lhs is a number; rhs is a finite with q==0 */
2710 /* ------------------------------------------------------------------ */
2711 /* decNumberShift -- shift the coefficient of a Number left or right */
2712 /* */
2713 /* This computes C = A << B or C = A >> -B (in base ten). */
2714 /* */
2715 /* res is C, the result. C may be A and/or B (e.g., X=X<<X) */
2716 /* lhs is A */
2717 /* rhs is B, the number of digits to shift (-ve to right) */
2718 /* set is the context */
2719 /* */
2720 /* The digits of the coefficient of A are shifted to the left (if B */
2721 /* is positive) or to the right (if B is negative) without adjusting */
2722 /* the exponent or the sign of A. */
2723 /* */
2724 /* B must be an integer (q=0) and in the range -set->digits through */
2725 /* +set->digits. */
2726 /* C must have space for set->digits digits. */
2727 /* NaNs are propagated as usual. Infinities are unaffected (but */
2728 /* B must be valid). No status is set unless B is invalid or an */
2729 /* operand is an sNaN. */
2730 /* ------------------------------------------------------------------ */
2736 #if DECCHECK
2738 #endif
2740 /* NaNs propagate as normal */
2743 /* rhs must be an integer */
2759 }
2761 /* first remove leading digits if necessary */
2764 /* that updated res->digits; may have gone to 1 (for a */
2765 /* single digit or for zero */
2766 }
2775 }
2779 }
2788 /* ------------------------------------------------------------------ */
2789 /* decNumberSquareRoot -- square root operator */
2790 /* */
2791 /* This computes C = squareroot(A) */
2792 /* */
2793 /* res is C, the result. C may be A */
2794 /* rhs is A */
2795 /* set is the context; note that rounding mode has no effect */
2796 /* */
2797 /* C must have space for set->digits digits. */
2798 /* ------------------------------------------------------------------ */
2799 /* This uses the following varying-precision algorithm in: */
2800 /* */
2801 /* Properly Rounded Variable Precision Square Root, T. E. Hull and */
2802 /* A. Abrham, ACM Transactions on Mathematical Software, Vol 11 #3, */
2803 /* pp229-237, ACM, September 1985. */
2804 /* */
2805 /* The square-root is calculated using Newton's method, after which */
2806 /* a check is made to ensure the result is correctly rounded. */
2807 /* */
2808 /* % [Reformatted original Numerical Turing source code follows.] */
2809 /* function sqrt(x : real) : real */
2810 /* % sqrt(x) returns the properly rounded approximation to the square */
2811 /* % root of x, in the precision of the calling environment, or it */
2812 /* % fails if x < 0. */
2813 /* % t e hull and a abrham, august, 1984 */
2814 /* if x <= 0 then */
2815 /* if x < 0 then */
2816 /* assert false */
2817 /* else */
2818 /* result 0 */
2819 /* end if */
2820 /* end if */
2821 /* var f := setexp(x, 0) % fraction part of x [0.1 <= x < 1] */
2822 /* var e := getexp(x) % exponent part of x */
2823 /* var approx : real */
2824 /* if e mod 2 = 0 then */
2825 /* approx := .259 + .819 * f % approx to root of f */
2826 /* else */
2827 /* f := f/l0 % adjustments */
2828 /* e := e + 1 % for odd */
2829 /* approx := .0819 + 2.59 * f % exponent */
2830 /* end if */
2831 /* */
2832 /* var p:= 3 */
2833 /* const maxp := currentprecision + 2 */
2834 /* loop */
2835 /* p := min(2*p - 2, maxp) % p = 4,6,10, . . . , maxp */
2836 /* precision p */
2837 /* approx := .5 * (approx + f/approx) */
2838 /* exit when p = maxp */
2839 /* end loop */
2840 /* */
2841 /* % approx is now within 1 ulp of the properly rounded square root */
2842 /* % of f; to ensure proper rounding, compare squares of (approx - */
2843 /* % l/2 ulp) and (approx + l/2 ulp) with f. */
2844 /* p := currentprecision */
2845 /* begin */
2846 /* precision p + 2 */
2847 /* const approxsubhalf := approx - setexp(.5, -p) */
2848 /* if mulru(approxsubhalf, approxsubhalf) > f then */
2849 /* approx := approx - setexp(.l, -p + 1) */
2850 /* else */
2851 /* const approxaddhalf := approx + setexp(.5, -p) */
2852 /* if mulrd(approxaddhalf, approxaddhalf) < f then */
2853 /* approx := approx + setexp(.l, -p + 1) */
2854 /* end if */
2855 /* end if */
2856 /* end */
2857 /* result setexp(approx, e div 2) % fix exponent */
2858 /* end sqrt */
2859 /* ------------------------------------------------------------------ */
2874 #if DECSUBSET
2876 #endif
2877 /* buffer for f [needs +1 in case DECBUFFER 0] */
2879 /* buffer for a [needs +2 to match likely maxp] */
2881 /* buffer for temporary, b [must be same size as a] */
2889 /* buffer for temporary variable, up to 3 digits */
2893 #if DECCHECK
2895 #endif
2898 #if DECSUBSET
2900 /* reduce operand and set lostDigits status, as needed */
2904 /* [Note: 'f' allocation below could reuse this buffer if */
2905 /* used, but as this is rare they are kept separate for clarity.] */
2907 }
2908 }
2909 #endif
2910 /* [following code does not require input rounding] */
2912 /* handle infinities and NaNs */
2917 }
2920 }
2922 /* calculate the ideal (preferred) exponent [floor(exp/2)] */
2923 /* [We would like to write: ideal=rhs->exponent>>1, but this */
2924 /* generates a compiler warning. Generated code is the same.] */
2927 /* handle zeros */
2931 /* use decFinish to clamp any out-of-range exponent, etc. */
2934 }
2936 /* any other -x is an oops */
2940 }
2942 /* space is needed for three working variables */
2943 /* f -- the same precision as the RHS, reduced to 0.01->0.99... */
2944 /* a -- Hull's approximation -- precision, when assigned, is */
2945 /* currentprecision+1 or the input argument precision, */
2946 /* whichever is larger (+2 for use as temporary) */
2947 /* b -- intermediate temporary result (same size as a) */
2948 /* if any is too long for local storage, then allocate */
2959 }
2960 /* a and b both need to be able to hold a maxp-length number */
2970 }
2972 /* copy rhs -> f, save exponent, and reduce so 0.1 <= f < 1 */
2977 /* set up working context */
2980 /* [Until further notice, no error is possible and status bits */
2981 /* (Rounded, etc.) should be ignored, not accumulated.] */
2983 /* Calculate initial approximation, and allow for odd exponent */
2988 /* Set t=0.259, a=0.819 */
2991 #if DECDPUN>=3
2994 #elif DECDPUN==2
2997 #else
3000 #endif
3001 }
3003 /* Set t=0.0819, a=2.59 */
3008 #if DECDPUN>=3
3011 #elif DECDPUN==2
3014 #else
3017 #endif
3018 }
3021 /* [a is now the initial approximation for sqrt(f), calculated with */
3022 /* currentprecision, which is also a's precision.] */
3024 /* the main calculation loop */
3031 /* set p to min(2*p - 2, maxp) [hence 3; or: 4, 6, 10, ... , maxp] */
3034 /* a = 0.5 * (a + f/a) */
3035 /* [calculated at p then rounded to currentprecision] */
3042 /* Here, 0.1 <= a < 1 [Hull], and a has maxp digits */
3043 /* now reduce to length, etc.; this needs to be done with a */
3044 /* having the correct exponent so as to handle subnormals */
3045 /* correctly */
3055 /* Overflow was possible if the input exponent was out-of-range, */
3056 /* in which case quit */
3061 }
3063 /* Preserve status except Inexact/Rounded */
3066 /* Carry out the Hull correction */
3069 /* a is now at final precision and within 1 ulp of the properly */
3070 /* rounded square root of f; to ensure proper rounding, compare */
3071 /* squares of (a - l/2 ulp) and (a + l/2 ulp) with f. */
3072 /* Here workset.digits=maxp and t=0.5, and a->digits determines */
3073 /* the ulp */
3081 /* this is the more common adjustment, though both are rare */
3085 /* assign to approx [round to length] */
3089 }
3099 /* assign to approx [round to length] */
3103 }
3104 }
3105 /* [no errors are possible in the above, and rounding/inexact during */
3106 /* estimation are irrelevant, so status was not accumulated] */
3108 /* Here, 0.1 <= a < 1 (still), so adjust back */
3111 /* count droppable zeros [after any subnormal rounding] by */
3112 /* trimming a copy */
3116 /* Set Inexact and Rounded. The answer can only be exact if */
3117 /* it is short enough so that squaring it could fit in workp digits, */
3118 /* and it cannot have trailing zeros due to clamping, so these are */
3119 /* the only (relatively rare) conditions a careful check is needed */
3122 }
3128 }
3133 /* here, dropped is the count of trailing zeros in 'a' */
3134 /* use closest exponent to ideal... */
3141 }
3146 }
3147 }
3148 }
3149 }
3150 }
3152 /* double-check Underflow, as perhaps the result could not have */
3153 /* been subnormal (initial argument too big), or it is now Exact */
3156 /* check if truly subnormal */
3159 #else
3161 #endif
3162 /* check if truly inexact */
3164 }
3172 #if DECSUBSET
3174 #endif
3176 #if DECCHECK
3178 #endif
3182 /* ------------------------------------------------------------------ */
3183 /* decNumberSubtract -- subtract two Numbers */
3184 /* */
3185 /* This computes C = A - B */
3186 /* */
3187 /* res is C, the result. C may be A and/or B (e.g., X=X-X) */
3188 /* lhs is A */
3189 /* rhs is B */
3190 /* set is the context */
3191 /* */
3192 /* C must have space for set->digits digits. */
3193 /* ------------------------------------------------------------------ */
3200 #if DECCHECK
3202 #endif
3206 /* ------------------------------------------------------------------ */
3207 /* decNumberToIntegralExact -- round-to-integral-value with InExact */
3208 /* decNumberToIntegralValue -- round-to-integral-value */
3209 /* */
3210 /* res is the result */
3211 /* rhs is input number */
3212 /* set is the context */
3213 /* */
3214 /* res must have space for any value of rhs. */
3215 /* */
3216 /* This implements the IEEE special operators and therefore treats */
3217 /* special values as valid. For finite numbers it returns */
3218 /* rescale(rhs, 0) if rhs->exponent is <0. */
3219 /* Otherwise the result is rhs (so no error is possible, except for */
3220 /* sNaN). */
3221 /* */
3222 /* The context is used for rounding mode and status after sNaN, but */
3223 /* the digits setting is ignored. The Exact version will signal */
3224 /* Inexact if the result differs numerically from rhs; the other */
3225 /* never signals Inexact. */
3226 /* ------------------------------------------------------------------ */
3229 decNumber dn;
3233 #if DECCHECK
3235 #endif
3237 /* handle infinities and NaNs */
3241 }
3243 /* have a finite number; no error possible (res must be big enough) */
3245 /* that was easy, but if negative exponent there is work to do... */
3252 }
3262 /* this never affects set, except for sNaNs; NaN will have been set */
3263 /* or propagated already, so no need to call decStatus */
3268 /* ------------------------------------------------------------------ */
3269 /* decNumberXor -- XOR two Numbers, digitwise */
3270 /* */
3271 /* This computes C = A ^ B */
3272 /* */
3273 /* res is C, the result. C may be A and/or B (e.g., X=X^X) */
3274 /* lhs is A */
3275 /* rhs is B */
3276 /* set is the context (used for result length and error report) */
3277 /* */
3278 /* C must have space for set->digits digits. */
3279 /* */
3280 /* Logical function restrictions apply (see above); a NaN is */
3281 /* returned with Invalid_operation if a restriction is violated. */
3282 /* ------------------------------------------------------------------ */
3289 #if DECCHECK
3291 #endif
3297 }
3298 /* operands are valid */
3315 /* This loop could be unrolled and/or use BIN2BCD tables */
3325 }
3330 /* [here uc-1 is the msu of the result] */
3338 /* ================================================================== */
3339 /* Utility routines */
3340 /* ================================================================== */
3342 /* ------------------------------------------------------------------ */
3343 /* decNumberClass -- return the decClass of a decNumber */
3344 /* dn -- the decNumber to test */
3345 /* set -- the context to use for Emin */
3346 /* returns the decClass enum */
3347 /* ------------------------------------------------------------------ */
3352 /* must be an infinity */
3355 }
3356 /* is finite */
3360 }
3361 /* is subnormal or zero */
3365 }
3370 /* ------------------------------------------------------------------ */
3371 /* decNumberClassToString -- convert decClass to a string */
3372 /* */
3373 /* eclass is a valid decClass */
3374 /* returns a constant string describing the class (max 13+1 chars) */
3375 /* ------------------------------------------------------------------ */
3390 /* ------------------------------------------------------------------ */
3391 /* decNumberCopy -- copy a number */
3392 /* */
3393 /* dest is the target decNumber */
3394 /* src is the source decNumber */
3395 /* returns dest */
3396 /* */
3397 /* (dest==src is allowed and is a no-op) */
3398 /* All fields are updated as required. This is a utility operation, */
3399 /* so special values are unchanged and no error is possible. */
3400 /* ------------------------------------------------------------------ */
3403 #if DECCHECK
3405 #endif
3409 /* Use explicit assignments here as structure assignment could copy */
3410 /* more than just the lsu (for small DECDPUN). This would not affect */
3411 /* the value of the results, but could disturb test harness spill */
3412 /* checking. */
3420 /* memcpy for the remaining Units would be safe as they cannot */
3421 /* overlap. However, this explicit loop is faster in short cases. */
3425 }
3429 /* ------------------------------------------------------------------ */
3430 /* decNumberCopyAbs -- quiet absolute value operator */
3431 /* */
3432 /* This sets C = abs(A) */
3433 /* */
3434 /* res is C, the result. C may be A */
3435 /* rhs is A */
3436 /* */
3437 /* C must have space for set->digits digits. */
3438 /* No exception or error can occur; this is a quiet bitwise operation.*/
3439 /* See also decNumberAbs for a checking version of this. */
3440 /* ------------------------------------------------------------------ */
3442 #if DECCHECK
3444 #endif
3450 /* ------------------------------------------------------------------ */
3451 /* decNumberCopyNegate -- quiet negate value operator */
3452 /* */
3453 /* This sets C = negate(A) */
3454 /* */
3455 /* res is C, the result. C may be A */
3456 /* rhs is A */
3457 /* */
3458 /* C must have space for set->digits digits. */
3459 /* No exception or error can occur; this is a quiet bitwise operation.*/
3460 /* See also decNumberMinus for a checking version of this. */
3461 /* ------------------------------------------------------------------ */
3463 #if DECCHECK
3465 #endif
3471 /* ------------------------------------------------------------------ */
3472 /* decNumberCopySign -- quiet copy and set sign operator */
3473 /* */
3474 /* This sets C = A with the sign of B */
3475 /* */
3476 /* res is C, the result. C may be A */
3477 /* lhs is A */
3478 /* rhs is B */
3479 /* */
3480 /* C must have space for set->digits digits. */
3481 /* No exception or error can occur; this is a quiet bitwise operation.*/
3482 /* ------------------------------------------------------------------ */
3486 #if DECCHECK
3488 #endif
3496 /* ------------------------------------------------------------------ */
3497 /* decNumberGetBCD -- get the coefficient in BCD8 */
3498 /* dn is the source decNumber */
3499 /* bcd is the uInt array that will receive dn->digits BCD bytes, */
3500 /* most-significant at offset 0 */
3501 /* returns bcd */
3502 /* */
3503 /* bcd must have at least dn->digits bytes. No error is possible; if */
3504 /* dn is a NaN or Infinite, digits must be 1 and the coefficient 0. */
3505 /* ------------------------------------------------------------------ */
3518 cut--;
3520 up++;
3523 }
3524 #endif
3528 /* ------------------------------------------------------------------ */
3529 /* decNumberSetBCD -- set (replace) the coefficient from BCD8 */
3530 /* dn is the target decNumber */
3531 /* bcd is the uInt array that will source n BCD bytes, most- */
3532 /* significant at offset 0 */
3533 /* n is the number of digits in the source BCD array (bcd) */
3534 /* returns dn */
3535 /* */
3536 /* dn must have space for at least n digits. No error is possible; */
3537 /* if dn is a NaN, or Infinite, or is to become a zero, n must be 1 */
3538 /* and bcd[0] zero. */
3539 /* ------------------------------------------------------------------ */
3547 /* calculate how many digits in msu, and hence first cut */
3553 }
3554 #endif
3559 /* ------------------------------------------------------------------ */
3560 /* decNumberIsNormal -- test normality of a decNumber */
3561 /* dn is the decNumber to test */
3562 /* set is the context to use for Emin */
3563 /* returns 1 if |dn| is finite and >=Nmin, 0 otherwise */
3564 /* ------------------------------------------------------------------ */
3567 #if DECCHECK
3569 #endif
3579 /* ------------------------------------------------------------------ */
3580 /* decNumberIsSubnormal -- test subnormality of a decNumber */
3581 /* dn is the decNumber to test */
3582 /* set is the context to use for Emin */
3583 /* returns 1 if |dn| is finite, non-zero, and <Nmin, 0 otherwise */
3584 /* ------------------------------------------------------------------ */
3587 #if DECCHECK
3589 #endif
3599 /* ------------------------------------------------------------------ */
3600 /* decNumberTrim -- remove insignificant zeros */
3601 /* */
3602 /* dn is the number to trim */
3603 /* returns dn */
3604 /* */
3605 /* All fields are updated as required. This is a utility operation, */
3606 /* so special values are unchanged and no error is possible. */
3607 /* ------------------------------------------------------------------ */
3611 #if DECCHECK
3613 #endif
3618 /* ------------------------------------------------------------------ */
3619 /* decNumberVersion -- return the name and version of this module */
3620 /* */
3621 /* No error is possible. */
3622 /* ------------------------------------------------------------------ */
3627 /* ------------------------------------------------------------------ */
3628 /* decNumberZero -- set a number to 0 */
3629 /* */
3630 /* dn is the number to set, with space for one digit */
3631 /* returns dn */
3632 /* */
3633 /* No error is possible. */
3634 /* ------------------------------------------------------------------ */
3635 /* Memset is not used as it is much slower in some environments. */
3638 #if DECCHECK
3640 #endif
3649 /* ================================================================== */
3650 /* Local routines */
3651 /* ================================================================== */
3653 /* ------------------------------------------------------------------ */
3654 /* decToString -- lay out a number into a string */
3655 /* */
3656 /* dn is the number to lay out */
3657 /* string is where to lay out the number */
3658 /* eng is 1 if Engineering, 0 if Scientific */
3659 /* */
3660 /* string must be at least dn->digits+14 characters long */
3661 /* No error is possible. */
3662 /* */
3663 /* Note that this routine can generate a -0 or 0.000. These are */
3664 /* never generated in subset to-number or arithmetic, but can occur */
3665 /* in non-subset arithmetic (e.g., -1*0 or 1.234-1.234). */
3666 /* ------------------------------------------------------------------ */
3667 /* If DECCHECK is enabled the string "?" is returned if a number is */
3668 /* invalid. */
3678 #if DECCHECK
3682 #endif
3686 c++;
3687 }
3693 /* a NaN */
3696 c++;
3697 }
3700 /* if not a clean non-zero coefficient, that's all there is in a */
3701 /* NaN string */
3703 /* [drop through to add integer] */
3704 }
3706 /* calculate how many digits in msu, and hence first cut */
3715 }
3719 /* non-0 exponent -- assume plain form */
3727 /* The C remainder operator is undefined for negative numbers, so */
3728 /* a positive remainder calculation must be used here */
3732 }
3735 }
3737 /* if dealing with zero still produce an exponent which is a */
3738 /* multiple of three, as expected, but there will only be the */
3739 /* one zero before the E, still. Otherwise note the padding. */
3745 }
3750 /* lay out the digits of the coefficient, adding 0s and . as needed */
3757 up--;
3760 }
3762 }
3768 up--;
3771 }
3773 }
3774 }
3776 }
3784 up--;
3787 }
3789 }
3790 }
3792 /* Finally add the E-part, if needed. It will never be 0, has a
3793 base maximum and minimum of +999999999 through -999999999, but
3794 could range down to -1999999998 for anormal numbers */
3803 }
3804 /* lay out the exponent [_itoa or equivalent is not ANSI C] */
3811 }
3816 /* ------------------------------------------------------------------ */
3817 /* decAddOp -- add/subtract operation */
3818 /* */
3819 /* This computes C = A + B */
3820 /* */
3821 /* res is C, the result. C may be A and/or B (e.g., X=X+X) */
3822 /* lhs is A */
3823 /* rhs is B */
3824 /* set is the context */
3825 /* negate is DECNEG if rhs should be negated, or 0 otherwise */
3826 /* status accumulates status for the caller */
3827 /* */
3828 /* C must have space for set->digits digits. */
3829 /* Inexact in status must be 0 for correct Exact zero sign in result */
3830 /* ------------------------------------------------------------------ */
3831 /* If possible, the coefficient is calculated directly into C. */
3832 /* However, if: */
3833 /* -- a digits+1 calculation is needed because the numbers are */
3834 /* unaligned and span more than set->digits digits */
3835 /* -- a carry to digits+1 digits looks possible */
3836 /* -- C is the same as A or B, and the result would destructively */
3837 /* overlap the A or B coefficient */
3838 /* then the result must be calculated into a temporary buffer. In */
3839 /* this case a local (stack) buffer is used if possible, and only if */
3840 /* too long for that does malloc become the final resort. */
3841 /* */
3842 /* Misalignment is handled as follows: */
3843 /* Apad: (AExp>BExp) Swap operands and proceed as for BExp>AExp. */
3844 /* BPad: Apply the padding by a combination of shifting (whole */
3845 /* units) and multiplication (part units). */
3846 /* */
3847 /* Addition, especially x=x+1, is speed-critical. */
3848 /* The static buffer is larger than might be expected to allow for */
3849 /* calls from higher-level functions (notably exp). */
3850 /* ------------------------------------------------------------------ */
3854 #if DECSUBSET
3857 #endif
3866 /* allocations when called from */
3867 /* other operations, notable exp] */
3872 #if DECCHECK
3874 #endif
3877 #if DECSUBSET
3879 /* reduce operands and set lostDigits status, as needed */
3884 }
3889 }
3890 }
3891 #endif
3892 /* [following code does not require input rounding] */
3894 /* note whether signs differ [used all paths] */
3897 /* handle infinities and NaNs */
3903 /* two infinities with different signs is invalid */
3907 }
3909 }
3916 }
3918 /* Quick exit for add 0s; return the non-0, modified as need be */
3926 #if DECSUBSET
3928 #endif
3929 /* exponent will be the lower of the two */
3933 /* 0-0 gives +0 unless rounding to -infinity, and -0-0 gives -0 */
3937 }
3938 }
3944 }
3947 }
3949 #if DECSUBSET
3951 #endif
3961 #if DECSUBSET
3963 #endif
3964 /* exponent will be the lower of the two */
3965 /* [0-0 case handled above] */
3971 }
3974 }
3975 #if DECSUBSET
3977 #endif
3981 /* [NB: both fastpath and mainpath code below assume these cases */
3982 /* (notably 0-0) have already been handled] */
3984 /* calculate the padding needed to align the operands */
3987 /* Fastpath cases where the numbers are aligned and normal, the RHS */
3988 /* is all in one unit, no operand rounding is needed, and no carry, */
3989 /* lengthening, or borrow is needed */
4005 }
4006 /* else drop out for careful add */
4007 }
4013 /* this could have reduced digits [but result>0] */
4016 }
4017 /* else drop out for careful subtract */
4018 }
4019 }
4021 /* Now align (pad) the lhs or rhs so they can be added or */
4022 /* subtracted, as necessary. If one number is much larger than */
4023 /* the other (that is, if in plain form there is a least one */
4024 /* digit between the lowest digit of one and the highest of the */
4025 /* other) padding with up to DIGITS-1 trailing zeros may be */
4026 /* needed; then apply rounding (as exotic rounding modes may be */
4027 /* affected by the residue). */
4032 /* [if padding==0 the operands are aligned; no padding is needed] */
4034 /* some padding needed; always pad the RHS, as any required */
4035 /* padding can then be effected by a simple combination of */
4036 /* shifts and a multiply */
4044 }
4046 /* If, after pad, rhs would be longer than lhs by digits+1 or */
4047 /* more then lhs cannot affect the answer, except as a residue, */
4048 /* so only need to pad up to a length of DIGITS+1. */
4050 /* The RHS is sufficient */
4051 /* for residue use the relative sign indication... */
4055 /* copy, shortening if necessary */
4057 /* if it was already shorter, then need to pad with zeros */
4061 }
4062 /* flip the result sign if unswapped and rhs was negated */
4067 /* LHS digits may affect result */
4074 /* determine the longer operand */
4078 /* Decide on the result buffer to use; if possible place directly */
4079 /* into result. */
4081 /* If destructive overlap, or the number is too long, or a carry or */
4082 /* borrow to DIGITS+1 might be possible, a buffer must be used. */
4083 /* [Might be worth more sophisticated tests when maxdigits==reqdigits] */
4086 /* buffer needed, choose it; units for maxdigits digits will be */
4087 /* needed, +1 Unit for carry or borrow */
4091 /* printf("malloc add %ld %ld\n", need, sizeof(accbuff)); */
4097 }
4098 }
4103 #if DECTRACE
4107 #endif
4109 /* add [A+B*m] or subtract [A+B*(-m)] */
4117 }
4118 #if DECTRACE
4120 #endif
4122 /* If a buffer was used the result must be copied back, possibly */
4123 /* shortening. (If no buffer was used then the result must have */
4124 /* fit, so can't need rounding and residue must be 0.) */
4127 #if DECSUBSET
4129 #endif
4130 /* remove leading zeros that were added due to rounding up to */
4131 /* integral Units -- before the test for rounding. */
4135 #if DECSUBSET
4136 }
4138 /* May have an underestimate. This only occurs when both */
4139 /* numbers fit in DECDPUN digits and are padding with a */
4140 /* negative multiple (-10, -100...) and the top digit(s) become */
4141 /* 0. (This only matters when using X3.274 rules where the */
4142 /* leading zero could be included in the rounding.) */
4146 }
4148 /* remove leading zeros that added due to rounding up to */
4149 /* integral Units (but only those in excess of the original */
4150 /* maxdigits length, unless extended) before test for rounding. */
4154 }
4155 }
4157 /* Now apply rounding if needed before removing leading zeros. */
4158 /* This is safe because subnormals are not a possibility */
4162 }
4164 #endif
4167 /* strip leading zeros [these were left on in case of subset subtract] */
4170 /* apply checks and rounding */
4173 /* "When the sum of two operands with opposite signs is exactly */
4174 /* zero, the sign of that sum shall be '+' in all rounding modes */
4175 /* except round toward -Infinity, in which mode that sign shall be */
4176 /* '-'." [Subset zeros also never have '-', set by decFinish.] */
4178 #if DECSUBSET
4180 #endif
4184 }
4188 #if DECSUBSET
4191 #endif
4195 /* ------------------------------------------------------------------ */
4196 /* decDivideOp -- division operation */
4197 /* */
4198 /* This routine performs the calculations for all four division */
4199 /* operators (divide, divideInteger, remainder, remainderNear). */
4200 /* */
4201 /* C=A op B */
4202 /* */
4203 /* res is C, the result. C may be A and/or B (e.g., X=X/X) */
4204 /* lhs is A */
4205 /* rhs is B */
4206 /* set is the context */
4207 /* op is DIVIDE, DIVIDEINT, REMAINDER, or REMNEAR respectively. */
4208 /* status is the usual accumulator */
4209 /* */
4210 /* C must have space for set->digits digits. */
4211 /* */
4212 /* ------------------------------------------------------------------ */
4213 /* The underlying algorithm of this routine is the same as in the */
4214 /* 1981 S/370 implementation, that is, non-restoring long division */
4215 /* with bi-unit (rather than bi-digit) estimation for each unit */
4216 /* multiplier. In this pseudocode overview, complications for the */
4217 /* Remainder operators and division residues for exact rounding are */
4218 /* omitted for clarity. */
4219 /* */
4220 /* Prepare operands and handle special values */
4221 /* Test for x/0 and then 0/x */
4222 /* Exp =Exp1 - Exp2 */
4223 /* Exp =Exp +len(var1) -len(var2) */
4224 /* Sign=Sign1 * Sign2 */
4225 /* Pad accumulator (Var1) to double-length with 0's (pad1) */
4226 /* Pad Var2 to same length as Var1 */
4227 /* msu2pair/plus=1st 2 or 1 units of var2, +1 to allow for round */
4228 /* have=0 */
4229 /* Do until (have=digits+1 OR residue=0) */
4230 /* if exp<0 then if integer divide/residue then leave */
4231 /* this_unit=0 */
4232 /* Do forever */
4233 /* compare numbers */
4234 /* if <0 then leave inner_loop */
4235 /* if =0 then (* quick exit without subtract *) do */
4236 /* this_unit=this_unit+1; output this_unit */
4237 /* leave outer_loop; end */
4238 /* Compare lengths of numbers (mantissae): */
4239 /* If same then tops2=msu2pair -- {units 1&2 of var2} */
4240 /* else tops2=msu2plus -- {0, unit 1 of var2} */
4241 /* tops1=first_unit_of_Var1*10**DECDPUN +second_unit_of_var1 */
4242 /* mult=tops1/tops2 -- Good and safe guess at divisor */
4243 /* if mult=0 then mult=1 */
4244 /* this_unit=this_unit+mult */
4245 /* subtract */
4246 /* end inner_loop */
4247 /* if have\=0 | this_unit\=0 then do */
4248 /* output this_unit */
4249 /* have=have+1; end */
4250 /* var2=var2/10 */
4251 /* exp=exp-1 */
4252 /* end outer_loop */
4253 /* exp=exp+1 -- set the proper exponent */
4254 /* if have=0 then generate answer=0 */
4255 /* Return (Result is defined by Var1) */
4256 /* */
4257 /* ------------------------------------------------------------------ */
4258 /* Two working buffers are needed during the division; one (digits+ */
4259 /* 1) to accumulate the result, and the other (up to 2*digits+1) for */
4260 /* long subtractions. These are acc and var1 respectively. */
4261 /* var1 is a copy of the lhs coefficient, var2 is the rhs coefficient.*/
4262 /* The static buffers may be larger than might be expected to allow */
4263 /* for calls from higher-level functions (notably exp). */
4264 /* ------------------------------------------------------------------ */
4268 #if DECSUBSET
4271 #endif
4305 #if DECSUBSET
4307 #endif
4309 #if DECCHECK
4311 #endif
4314 #if DECSUBSET
4316 /* reduce operands and set lostDigits status, as needed */
4321 }
4326 }
4327 }
4328 #endif
4329 /* [following code does not require input rounding] */
4333 /* handle infinities and NaNs */
4338 }
4339 /* one or two infinities */
4345 }
4346 /* [Note that infinity/0 raises no exceptions] */
4350 }
4354 /* result is [finished clone of] lhs */
4356 }
4360 /* for DIVIDEINT the exponent is always 0. For DIVIDE, result */
4361 /* is a 0 with infinitely negative exponent, clamped to minimum */
4365 }
4366 }
4369 }
4370 }
4372 /* handle 0 rhs (x/0) */
4377 }
4384 }
4385 }
4388 /* handle 0 lhs (0/x) */
4390 #if DECSUBSET
4393 #endif
4401 }
4405 }
4410 }
4411 #if DECSUBSET
4412 }
4413 #endif
4416 /* Precalculate exponent. This starts off adjusted (and hence fits */
4417 /* in 31 bits) and becomes the usual unadjusted exponent as the */
4418 /* division proceeds. The order of evaluation is important, here, */
4419 /* to avoid wrap. */
4422 /* If the working exponent is -ve, then some quick exits are */
4423 /* possible because the quotient is known to be <1 */
4424 /* [for REMNEAR, it needs to be < -1, as -0.5 could need work] */
4428 #if DECSUBSET
4430 #endif
4433 /* fastpath remainders so long as the lhs has the smaller */
4434 /* (or equal) exponent */
4437 /* It is REMAINDER or safe REMNEAR; result is [finished */
4438 /* clone of] lhs (r = x - 0*y) */
4443 }
4444 /* [unsafe REMNEAR drops through] */
4445 }
4448 /* Long (slow) division is needed; roll up the sleeves... */
4450 /* The accumulator will hold the quotient of the division. */
4451 /* If it needs to be too long for stack storage, then allocate. */
4454 /* printf("malloc dvacc %ld units\n", acclength); */
4460 }
4462 /* var1 is the padded LHS ready for subtractions. */
4463 /* If it needs to be too long for stack storage, then allocate. */
4464 /* The maximum units needed for var1 (long subtraction) is: */
4465 /* Enough for */
4466 /* (rhs->digits+reqdigits-1) -- to allow full slide to right */
4467 /* or (lhs->digits) -- to allow for long lhs */
4468 /* whichever is larger */
4469 /* +1 -- for rounding of slide to right */
4470 /* +1 -- for leading 0s */
4471 /* +1 -- for pre-adjust if a remainder or DIVIDEINT */
4472 /* [Note: unused units do not participate in decUnitAddSub data] */
4476 /* allocate a guard unit above msu1 for REMAINDERNEAR */
4479 /* printf("malloc dvvar %ld units\n", var1units+1); */
4485 }
4487 /* Extend the lhs and rhs to full long subtraction length. The lhs */
4488 /* is truly extended into the var1 buffer, with 0 padding, so a */
4489 /* subtract in place is always possible. The rhs (var2) has */
4490 /* virtual padding (implemented by decUnitAddSub). */
4491 /* One guard unit was allocated above msu1 for rem=rem+rem in */
4492 /* REMAINDERNEAR. */
4498 /* rhs (var2) is left-aligned with var1 at the start */
4503 /* now set up the variables which will be used for estimating the */
4504 /* multiplication factor. If these variables are not exact, add */
4505 /* 1 to make sure that the multiplier is never overestimated. */
4512 }
4514 /* The calculation is working in units, which may have leading zeros, */
4515 /* but the exponent was calculated on the assumption that they are */
4516 /* both left-aligned. Adjust the exponent to compensate: add the */
4517 /* number of leading zeros in var1 msu and subtract those in var2 msu. */
4518 /* [This is actually done by counting the digits and negating, as */
4519 /* lead1=DECDPUN-digits1, and similarly for lead2.] */
4523 /* Now, if doing an integer divide or remainder, ensure that */
4524 /* the result will be Unit-aligned. To do this, shift the var1 */
4525 /* accumulator towards least if need be. (It's much easier to */
4526 /* do this now than to reassemble the residue afterwards, if */
4527 /* doing a remainder.) Also ensure the exponent is not negative. */
4530 /* save the initial 'false' padding of var1, in digits */
4532 /* Determine the shift to do. */
4538 /* clean any most-significant units which were just emptied */
4543 /* optimization: if the first iteration will just produce 0, */
4544 /* preadjust to skip it [valid for DIVIDE only] */
4548 }
4549 }
4551 /* ---- start the long-division loops ------------------------------ */
4557 /* find the next unit */
4559 /* strip leading zero units [from either pre-adjust or from */
4560 /* subtract last time around]. Leave at least one unit. */
4565 /* compare the two numbers, from msu */
4570 /* v1=*pv1 -- always OK */
4575 }
4576 /* here when all inspected or a difference seen */
4579 /* reach here if var1 and var2 are identical; subtraction */
4580 /* would increase digit by one, and the residue will be 0 so */
4581 /* the calculation is done; leave the loop with residue=0. */
4587 /* *pv1>v2. Prepare for real subtraction; the lengths are equal */
4588 /* Estimate the multiplier (there's always a msu1-1)... */
4589 /* Bring in two units of var2 to provide a good estimate. */
4593 /* The var2 msu is one unit towards the lsu of the var1 msu, */
4594 /* so only one unit for var2 can be used. */
4596 }
4598 /* subtraction needed; var1 is > var2 */
4600 /* subtract var1-var2, into var1; only the overlap needs */
4601 /* processing, as this is an in-place calculation */
4603 #if DECTRACE
4607 #endif
4611 #if DECTRACE
4613 #endif
4614 /* var1 now probably has leading zeros; these are removed at the */
4615 /* top of the inner loop. */
4618 /* The next unit has been calculated in full; unless it's a */
4619 /* leading zero, add to acc */
4622 /* account exactly for the new digits */
4626 }
4631 }
4633 /* if the residue is zero, the operation is done (unless divide */
4634 /* or divideInteger and still not enough digits yet) */
4638 /* [drop through if divideInteger] */
4639 }
4640 /* also done enough if calculating remainder or integer */
4641 /* divide and just did the last ('units') unit */
4644 /* to get here, var1 is less than var2, so divide var2 by the per- */
4645 /* Unit power of ten and go for the next digit */
4650 /* ---- division is complete --------------------------------------- */
4651 /* here: acc has at least reqdigits+1 of good results (or fewer */
4652 /* if early stop), starting at accnext+1 (its lsu) */
4653 /* var1 has any residue at the stopping point */
4654 /* accunits is the number of digits collected in acc */
4659 }
4661 /* accnext now -> lowest unit of result */
4665 /* record the presence of any residue, for rounding */
4668 /* Had an exact division; clean up spurious trailing 0s. */
4669 /* There will be at most DECDPUN-1, from the final multiply, */
4670 /* and then only if the result is non-0 (and even) and the */
4671 /* exponent is 'loose'. */
4672 #if DECDPUN>1
4675 /* count the trailing zeros */
4679 #if DECDPUN<=4
4682 #else
4684 #endif
4685 exponent++;
4686 }
4691 /* [exponent was adjusted in the loop] */
4692 }
4694 #endif
4698 /* check for coefficient overflow */
4702 }
4704 /* [Here, the exponent will be 0, because var1 was adjusted */
4705 /* appropriately.] */
4713 /* Fastpath when residue is truly 0 is worthwhile [and */
4714 /* simplifies the code below] */
4719 #if DECSUBSET
4721 #endif
4726 }
4727 /* note if the quotient was odd */
4732 /* treat the residue, in var1, as the value to return, via acc */
4733 /* calculate the unused zero digits. This is the smaller of: */
4734 /* var1 initial padding (saved above) */
4735 /* var2 residual padding, which happens to be given by: */
4737 /* [the 'exponent' term accounts for the shifts during divide] */
4740 /* shift var1 the requested amount, and adjust its digits */
4749 /* Now correct the result if doing remainderNear; if it */
4750 /* (looking just at coefficients) is > rhs/2, or == rhs/2 and */
4751 /* the integer was odd then the result should be rem-rhs. */
4755 /* calculate remainder*2 into the var1 buffer (which has */
4756 /* 'headroom' of an extra unit and hence enough space) */
4757 /* [a dedicated 'double' loop would be faster, here] */
4760 /* decDumpAr('r', accnext, tarunits); */
4762 /* Here, accnext (var1) holds tarunits Units with twice the */
4763 /* remainder's coefficient, which must now be compared to the */
4764 /* RHS. The remainder's exponent may be smaller than the RHS's. */
4771 /* now restore the remainder by dividing by two; the lsu */
4772 /* is known to be even. */
4779 }
4780 /* [accunits still describes the original remainder length] */
4784 /* This is effectively causing round-up of the quotient, */
4785 /* so if it was the rare case where it was full and all */
4786 /* nines, it would overflow and hence division-impossible */
4787 /* should be raised */
4793 }
4797 }
4805 /* rem-rhs is needed; the sign will invert. Again, var1 */
4806 /* can safely be used for the working Units array. */
4808 /* Calculate units and remainder from exponent. */
4811 /* subtract [A+B*(-m)]; the result will always be negative */
4817 /* [exponent is as for original remainder] */
4819 }
4824 /* Set exponent and bits */
4828 /* Now the coefficient. */
4833 #if DECSUBSET
4834 /* If a divide then strip trailing zeros if subset [after round] */
4836 #endif
4841 #if DECSUBSET
4844 #endif
4848 /* ------------------------------------------------------------------ */
4849 /* decMultiplyOp -- multiplication operation */
4850 /* */
4851 /* This routine performs the multiplication C=A x B. */
4852 /* */
4853 /* res is C, the result. C may be A and/or B (e.g., X=X*X) */
4854 /* lhs is A */
4855 /* rhs is B */
4856 /* set is the context */
4857 /* status is the usual accumulator */
4858 /* */
4859 /* C must have space for set->digits digits. */
4860 /* */
4861 /* ------------------------------------------------------------------ */
4862 /* 'Classic' multiplication is used rather than Karatsuba, as the */
4863 /* latter would give only a minor improvement for the short numbers */
4864 /* expected to be handled most (and uses much more memory). */
4865 /* */
4866 /* There are two major paths here: the general-purpose ('old code') */
4867 /* path which handles all DECDPUN values, and a fastpath version */
4868 /* which is used if 64-bit ints are available, DECDPUN<=4, and more */
4869 /* than two calls to decUnitAddSub would be made. */
4870 /* */
4871 /* The fastpath version lumps units together into 8-digit or 9-digit */
4872 /* chunks, and also uses a lazy carry strategy to minimise expensive */
4873 /* 64-bit divisions. The chunks are then broken apart again into */
4874 /* units for continuing processing. Despite this overhead, the */
4875 /* fastpath can speed up some 16-digit operations by 10x (and much */
4876 /* more for higher-precision calculations). */
4877 /* */
4878 /* A buffer always has to be used for the accumulator; in the */
4879 /* fastpath, buffers are also always needed for the chunked copies of */
4880 /* of the operand coefficients. */
4881 /* Static buffers are larger than needed just for multiply, to allow */
4882 /* for calls from other operations (notably exp). */
4883 /* ------------------------------------------------------------------ */
4884 #define FASTMUL (DECUSE64 && DECDPUN<5)
4896 /* *4 for calls from other operations) */
4901 #if FASTMUL
4902 /* if DECDPUN is 1 or 3 work in base 10**9, otherwise */
4903 /* (DECDPUN is 2 or 4) then work in base 10**8 */
4908 #else
4909 #define FASTBASE 100000000
4910 #define FASTDIGS 8
4912 #endif
4913 /* three buffers are used, two for chunked copies of the operands */
4914 /* (base 10**8 or base 10**9) and one base 2**64 accumulator with */
4915 /* lazy carry evaluation */
4923 /* [allocacc is shared for both paths, as only one will run] */
4925 #if DECDPUN==1
4927 #endif
4939 #endif
4941 #if DECSUBSET
4944 #endif
4946 #if DECCHECK
4948 #endif
4950 /* precalculate result sign */
4953 /* handle infinities and NaNs */
4958 /* one or two infinities; Infinity * 0 is invalid */
4967 /* For best speed, as in DMSRCN [the original Rexx numerics */
4968 /* module], use the shorter number as the multiplier (rhs) and */
4969 /* the longer as the multiplicand (lhs) to minimise the number of */
4970 /* adds (partial products) */
4975 }
4978 #if DECSUBSET
4980 /* reduce operands and set lostDigits status, as needed */
4985 }
4990 }
4991 }
4992 #endif
4993 /* [following code does not require input rounding] */
4996 /* use the fast path if there are enough digits in the shorter */
4997 /* operand to make the setup and takedown worthwhile */
5000 /* calculate the number of elements in each array */
5005 /* allocate buffers if required, as usual */
5015 /* Allocating the accumulator space needs a special case when */
5016 /* DECDPUN=1 because when converting the accumulator to Units */
5017 /* after the multiplication each 8-byte item becomes 9 1-byte */
5018 /* units. Therefore iacc extra bytes are needed at the front */
5019 /* (rounded up to a multiple of 8 bytes), and the uLong */
5020 /* accumulator starts offset the appropriate number of units */
5021 /* to the right to avoid overwrite during the unchunking. */
5023 #if DECDPUN==1
5026 #endif
5035 #if DECDPUN==1
5037 #endif
5039 /* assemble the chunked copies of the left and right sides */
5051 /* zero the accumulator */
5054 /* Start the multiplication */
5055 /* Resolving carries can dominate the cost of accumulating the */
5056 /* partial products, so this is only done when necessary. */
5057 /* Each uLong item in the accumulator can hold values up to */
5058 /* 2**64-1, and each partial product can be as large as */
5059 /* (10**FASTDIGS-1)**2. When FASTDIGS=9, this can be added to */
5060 /* itself 18.4 times in a uLong without overflowing, so during */
5061 /* the main calculation resolution is carried out every 18th */
5062 /* add -- every 162 digits. Similarly, when FASTDIGS=8, the */
5063 /* partial products can be added to themselves 1844.6 times in */
5064 /* a uLong without overflowing, so intermediate carry */
5065 /* resolution occurs only every 14752 digits. Hence for common */
5066 /* short numbers usually only the one final carry resolution */
5067 /* occurs. */
5068 /* (The count is set via FASTLAZY to simplify experiments to */
5069 /* measure the value of this approach: a 35% improvement on a */
5070 /* [34x34] multiply.) */
5077 lazy--;
5080 /* spin up the accumulator resolving overflows */
5084 /* lcarry can exceed 2**32-1, so check again; this check */
5085 /* and occasional extra divide (slow) is well worth it, as */
5086 /* it allows FASTLAZY to be increased to 18 rather than 4 */
5087 /* in the FASTDIGS=9 case */
5094 }
5100 /* The multiplication is complete; time to convert back into */
5101 /* units. This can be done in-place in the accumulator and in */
5102 /* 32-bit operations, because carries were resolved after the */
5103 /* final add. This needs N-1 divides and multiplies for */
5104 /* each item in the accumulator (which will become up to N */
5105 /* units, where 2<=N<=9). */
5116 }
5118 #endif
5120 /* if accumulator will be too long for local storage, then allocate */
5127 }
5129 /* Now the main long multiplication loop */
5130 /* Unlike the equivalent in the IBM Java implementation, there */
5131 /* is no advantage in calculating from msu to lsu. So, do it */
5132 /* by the book, as it were. */
5133 /* Each iteration calculates ACC=ACC+MULTAND*MULT */
5141 /* Here, *mer is the next Unit in the multiplier to use */
5142 /* If non-zero [optimization] add it... */
5149 accunits++;
5150 }
5151 /* multiply multiplicand by 10**DECDPUN for next Unit to left */
5154 #if FASTMUL
5156 #endif
5157 /* common end-path */
5158 #if DECTRACE
5160 #endif
5162 /* acc now contains the exact result of the multiplication, */
5163 /* possibly with a leading zero unit; build the decNumber from */
5164 /* it, noting if any residue */
5168 /* There can be a 31-bit wrap in calculating the exponent. */
5169 /* This can only happen if both input exponents are negative and */
5170 /* both their magnitudes are large. If there was a wrap, set a */
5171 /* safe very negative exponent, from which decFinalize() will */
5172 /* raise a hard underflow shortly. */
5179 /* Set the coefficient. If any rounding, residue records */
5185 #if DECSUBSET
5188 #endif
5189 #if FASTMUL
5192 #endif
5196 /* ------------------------------------------------------------------ */
5197 /* decExpOp -- effect exponentiation */
5198 /* */
5199 /* This computes C = exp(A) */
5200 /* */
5201 /* res is C, the result. C may be A */
5202 /* rhs is A */
5203 /* set is the context; note that rounding mode has no effect */
5204 /* */
5205 /* C must have space for set->digits digits. status is updated but */
5206 /* not set. */
5207 /* */
5208 /* Restrictions: */
5209 /* */
5210 /* digits, emax, and -emin in the context must be less than */
5211 /* 2*DEC_MAX_MATH (1999998), and the rhs must be within these */
5212 /* bounds or a zero. This is an internal routine, so these */
5213 /* restrictions are contractual and not enforced. */
5214 /* */
5215 /* A finite result is rounded using DEC_ROUND_HALF_EVEN; it will */
5216 /* almost always be correctly rounded, but may be up to 1 ulp in */
5217 /* error in rare cases. */
5218 /* */
5219 /* Finite results will always be full precision and Inexact, except */
5220 /* when A is a zero or -Infinity (giving 1 or 0 respectively). */
5221 /* ------------------------------------------------------------------ */
5222 /* This approach used here is similar to the algorithm described in */
5223 /* */
5224 /* Variable Precision Exponential Function, T. E. Hull and */
5225 /* A. Abrham, ACM Transactions on Mathematical Software, Vol 12 #2, */
5226 /* pp79-91, ACM, June 1986. */
5227 /* */
5228 /* with the main difference being that the iterations in the series */
5229 /* evaluation are terminated dynamically (which does not require the */
5230 /* extra variable-precision variables which are expensive in this */
5231 /* context). */
5232 /* */
5233 /* The error analysis in Hull & Abrham's paper applies except for the */
5234 /* round-off error accumulation during the series evaluation. This */
5235 /* code does not precalculate the number of iterations and so cannot */
5236 /* use Horner's scheme. Instead, the accumulation is done at double- */
5237 /* precision, which ensures that the additions of the terms are exact */
5238 /* and do not accumulate round-off (and any round-off errors in the */
5239 /* terms themselves move 'to the right' faster than they can */
5240 /* accumulate). This code also extends the calculation by allowing, */
5241 /* in the spirit of other decNumber operators, the input to be more */
5242 /* precise than the result (the precision used is based on the more */
5243 /* precise of the input or requested result). */
5244 /* */
5245 /* Implementation notes: */
5246 /* */
5247 /* 1. This is separated out as decExpOp so it can be called from */
5248 /* other Mathematical functions (notably Ln) with a wider range */
5249 /* than normal. In particular, it can handle the slightly wider */
5250 /* (double) range needed by Ln (which has to be able to calculate */
5251 /* exp(-x) where x can be the tiniest number (Ntiny). */
5252 /* */
5253 /* 2. Normalizing x to be <=0.1 (instead of <=1) reduces loop */
5254 /* iterations by approximately a third with additional (although */
5255 /* diminishing) returns as the range is reduced to even smaller */
5256 /* fractions. However, h (the power of 10 used to correct the */
5257 /* result at the end, see below) must be kept <=8 as otherwise */
5258 /* the final result cannot be computed. Hence the leverage is a */
5259 /* sliding value (8-h), where potentially the range is reduced */
5260 /* more for smaller values. */
5261 /* */
5262 /* The leverage that can be applied in this way is severely */
5263 /* limited by the cost of the raise-to-the power at the end, */
5264 /* which dominates when the number of iterations is small (less */
5265 /* than ten) or when rhs is short. As an example, the adjustment */
5266 /* x**10,000,000 needs 31 multiplications, all but one full-width. */
5267 /* */
5268 /* 3. The restrictions (especially precision) could be raised with */
5269 /* care, but the full decNumber range seems very hard within the */
5270 /* 32-bit limits. */
5271 /* */
5272 /* 4. The working precisions for the static buffers are twice the */
5273 /* obvious size to allow for calls from decNumberPower. */
5274 /* ------------------------------------------------------------------ */
5286 /* the argument is often copied to normalize it, so (unusually) it */
5287 /* is treated like other buffers, using DECBUFFER, +1 in case */
5288 /* DECBUFFER is 0 */
5292 /* the working precision will be no more than set->digits+8+1 */
5293 /* so for on-stack buffers DECBUFFER+9 is used, +1 in case DECBUFFER */
5294 /* is 0 (and twice that for the accumulator) */
5296 /* buffer for t, term (working precision plus) */
5300 /* buffer for a, accumulator (working precision * 2), at least 9 */
5304 /* decNumber for the divisor term; this needs at most 9 digits */
5305 /* and so can be fixed size [16 so can use standard context] */
5310 #if DECCHECK
5313 #endif
5321 }
5330 /* e**x when 0 < x < 0.66 is < 1+3x/2, hence can fast-path */
5331 /* positive and negative tiny cases which will result in inexact */
5332 /* 1. This also allows the later add-accumulate to always be */
5333 /* exact (because its length will never be more than twice the */
5334 /* working precision). */
5335 /* The comparator (tiny) needs just one digit, so use the */
5336 /* decNumber d for it (reused as the divisor, etc., below); its */
5337 /* exponent is such that if x is positive it will have */
5338 /* set->digits-1 zeros between the decimal point and the digit, */
5339 /* which is 4, and if x is negative one more zero there as the */
5340 /* more precise result will be of the form 0.9999999 rather than */
5341 /* 1.0000001. Hence, tiny will be 0.0000004 if digits=7 and x>0 */
5342 /* or 0.00000004 if digits=7 and x<0. If RHS not larger than */
5343 /* this then the result will be 1.000000 */
5361 /* set up the context to be used for calculating a, as this is */
5362 /* used on both paths below */
5364 /* accumulator bounds are as requested (could underflow) */
5369 /* calculate the adjusted (Hull & Abrham) exponent (where the */
5370 /* decimal point is just to the left of the coefficient msd) */
5372 /* if h>8 then 10**h cannot be calculated safely; however, when */
5373 /* h=8 then exp(|rhs|) will be at least exp(1E+7) which is at */
5374 /* least 6.59E+4342944, so (due to the restriction on Emax/Emin) */
5375 /* overflow (or underflow to 0) is guaranteed -- so this case can */
5376 /* be handled by simply forcing the appropriate excess */
5378 /* set up here so Power call below will over or underflow to */
5379 /* zero; set accumulator to either 2 or 0.02 */
5380 /* [stack buffer for a is always big enough for this] */
5386 }
5389 /* [could/should increase this for precisions >40 or so, too] */
5391 /* if h is 8, cannot normalize to a lower upper limit because */
5392 /* the final result will not be computable (see notes above), */
5393 /* but leverage can be applied whenever h is less than 8. */
5394 /* Apply as much as possible, up to a MAXLEVER digits, which */
5395 /* sets the tradeoff against the cost of the later a**(10**h). */
5396 /* As h is increased, the working precision below also */
5397 /* increases to compensate for the "constant digits at the */
5398 /* front" effect. */
5405 }
5406 /* Take a copy of RHS if it needs normalization (true whenever x>=1) */
5416 }
5420 /* decNumberShow(x); */
5421 }
5423 /* Now use the usual power series to evaluate exp(x). The */
5424 /* series starts as 1 + x + x^2/2 ... so prime ready for the */
5425 /* third term by setting the term variable t=x, the accumulator */
5426 /* a=1, and the divisor d=2. */
5428 /* First determine the working precision. From Hull & Abrham */
5429 /* this is set->digits+h+2. However, if x is 'over-precise' we */
5430 /* need to allow for all its digits to potentially participate */
5431 /* (consider an x where all the excess digits are 9s) so in */
5432 /* this case use x->digits+h+2 */
5435 /* a and t are variable precision, and depend on p, so space */
5436 /* must be allocated for them if necessary */
5438 /* the accumulator needs to be able to hold 2p digits so that */
5439 /* the additions on the second and subsequent iterations are */
5440 /* sufficiently exact. */
5448 }
5449 /* the term needs to be able to hold p digits (which is */
5450 /* guaranteed to be larger than x->digits, so the initial copy */
5451 /* is safe); it may also be used for the raise-to-power */
5452 /* calculation below, which needs an extra two digits */
5460 }
5467 /* set up the contexts for calculating a, t, and d */
5470 /* accumulator bounds are set above, set precision now */
5472 /* term bounds avoid any underflow or overflow */
5475 /* [dset.digits=16, etc., are sufficient] */
5477 /* finally ready to roll */
5479 #if DECCHECK
5480 iterations++;
5481 #endif
5482 /* only the status from the accumulation is interesting */
5483 /* [but it should remain unchanged after first add] */
5487 /* the iteration ends when the term cannot affect the result, */
5488 /* if rounded to p digits, which is when its value is smaller */
5489 /* than the accumulator by p+1 digits. There must also be */
5490 /* full precision in a. */
5496 #if DECCHECK
5497 /* just a sanity check; comment out test to show always */
5501 #endif
5504 /* apply postconditioning: a=a**(10**h) -- this is calculated */
5505 /* at a slightly higher precision than Hull & Abrham suggest */
5511 /* avoid the overhead and many extra digits of decNumberPower */
5512 /* as all that is needed is the short 'multipliers' loop; here */
5513 /* accumulate the answer into t */
5516 /* abandon if have had overflow or terminal underflow */
5523 }
5528 /* decNumberShow(t); */
5530 }
5532 /* Copy and round the result to res */
5543 /* [status is handled by caller] */
5547 /* ------------------------------------------------------------------ */
5548 /* Initial-estimate natural logarithm table */
5549 /* */
5550 /* LNnn -- 90-entry 16-bit table for values from .10 through .99. */
5551 /* The result is a 4-digit encode of the coefficient (c=the */
5552 /* top 14 bits encoding 0-9999) and a 2-digit encode of the */
5553 /* exponent (e=the bottom 2 bits encoding 0-3) */
5554 /* */
5555 /* The resulting value is given by: */
5556 /* */
5557 /* v = -c * 10**(-e-3) */
5558 /* */
5559 /* where e and c are extracted from entry k = LNnn[x-10] */
5560 /* where x is truncated (NB) into the range 10 through 99, */
5561 /* and then c = k>>2 and e = k&3. */
5562 /* ------------------------------------------------------------------ */
5575 /* ------------------------------------------------------------------ */
5576 /* decLnOp -- effect natural logarithm */
5577 /* */
5578 /* This computes C = ln(A) */
5579 /* */
5580 /* res is C, the result. C may be A */
5581 /* rhs is A */
5582 /* set is the context; note that rounding mode has no effect */
5583 /* */
5584 /* C must have space for set->digits digits. */
5585 /* */
5586 /* Notable cases: */
5587 /* A<0 -> Invalid */
5588 /* A=0 -> -Infinity (Exact) */
5589 /* A=+Infinity -> +Infinity (Exact) */
5590 /* A=1 exactly -> 0 (Exact) */
5591 /* */
5592 /* Restrictions (as for Exp): */
5593 /* */
5594 /* digits, emax, and -emin in the context must be less than */
5595 /* DEC_MAX_MATH+11 (1000010), and the rhs must be within these */
5596 /* bounds or a zero. This is an internal routine, so these */
5597 /* restrictions are contractual and not enforced. */
5598 /* */
5599 /* A finite result is rounded using DEC_ROUND_HALF_EVEN; it will */
5600 /* almost always be correctly rounded, but may be up to 1 ulp in */
5601 /* error in rare cases. */
5602 /* ------------------------------------------------------------------ */
5603 /* The result is calculated using Newton's method, with each */
5604 /* iteration calculating a' = a + x * exp(-a) - 1. See, for example, */
5605 /* Epperson 1989. */
5606 /* */
5607 /* The iteration ends when the adjustment x*exp(-a)-1 is tiny enough. */
5608 /* This has to be calculated at the sum of the precision of x and the */
5609 /* working precision. */
5610 /* */
5611 /* Implementation notes: */
5612 /* */
5613 /* 1. This is separated out as decLnOp so it can be called from */
5614 /* other Mathematical functions (e.g., Log 10) with a wider range */
5615 /* than normal. In particular, it can handle the slightly wider */
5616 /* (+9+2) range needed by a power function. */
5617 /* */
5618 /* 2. The speed of this function is about 10x slower than exp, as */
5619 /* it typically needs 4-6 iterations for short numbers, and the */
5620 /* extra precision needed adds a squaring effect, twice. */
5621 /* */
5622 /* 3. Fastpaths are included for ln(10) and ln(2), up to length 40, */
5623 /* as these are common requests. ln(10) is used by log10(x). */
5624 /* */
5625 /* 4. An iteration might be saved by widening the LNnn table, and */
5626 /* would certainly save at least one if it were made ten times */
5627 /* bigger, too (for truncated fractions 0.100 through 0.999). */
5628 /* However, for most practical evaluations, at least four or five */
5629 /* iterations will be needed -- so this would only speed up by */
5630 /* 20-25% and that probably does not justify increasing the table */
5631 /* size. */
5632 /* */
5633 /* 5. The static buffers are larger than might be expected to allow */
5634 /* for calls from decNumberPower. */
5635 /* ------------------------------------------------------------------ */
5646 /* buffers for a (accumulator, typically precision+2) and b */
5647 /* (adjustment calculator, same size) */
5659 #if DECCHECK
5662 #endif
5670 }
5679 /* Non-zero negatives are bad... */
5684 /* Here, rhs is positive, finite, and in range */
5686 /* lookaside fastpath code for ln(2) and ln(10) at common lengths */
5688 #if DECDPUN==1
5690 #else
5692 #endif
5706 /* Determine the working precision. This is normally the */
5707 /* requested precision + 2, with a minimum of 9. However, if */
5708 /* the rhs is 'over-precise' then allow for all its digits to */
5709 /* potentially participate (consider an rhs where all the excess */
5710 /* digits are 9s) so in this case use rhs->digits+2. */
5713 /* Allocate space for the accumulator and the high-precision */
5714 /* adjustment calculator, if necessary. The accumulator must */
5715 /* be able to hold p digits, and the adjustment up to */
5716 /* rhs->digits+p digits. They are also made big enough for 16 */
5717 /* digits so that they can be used for calculating the initial */
5718 /* estimate. */
5726 }
5735 }
5737 /* Prepare an initial estimate in acc. Calculate this by */
5738 /* considering the coefficient of x to be a normalized fraction, */
5739 /* f, with the decimal point at far left and multiplied by */
5740 /* 10**r. Then, rhs=f*10**r and 0.1<=f<1, and */
5741 /* ln(x) = ln(f) + ln(10)*r */
5742 /* Get the initial estimate for ln(f) from a small lookup */
5743 /* table (see above) indexed by the first two digits of f, */
5744 /* truncated. */
5752 /* now get top two digits of rhs into b by simple truncate and */
5753 /* force to integer */
5766 /* the initial estimate is now in a, with up to 4 digits correct. */
5767 /* When rhs is at or near Nmax the estimate will be low, so we */
5768 /* will approach it from below, avoiding overflow when calling exp. */
5772 /* accumulator bounds are as requested (could underflow, but */
5773 /* cannot overflow) */
5777 /* set up a context to be used for the multiply and subtract */
5781 /* [see decExpOp call below] */
5782 /* for each iteration double the number of digits to calculate, */
5783 /* up to a maximum of p */
5785 /* [initially 9 as then the sequence starts 7+2, 16+2, and */
5786 /* 34+2, which is ideal for standard-sized numbers] */
5790 #if DECCHECK
5791 iterations++;
5793 #endif
5794 /* calculate the adjustment (exp(-a)*x-1) into b. This is a */
5795 /* catastrophic subtraction but it really is the difference */
5796 /* from 1 that is of interest. */
5797 /* Use the internal entry point to Exp as it allows the double */
5798 /* range for calculating exp(-a) when a is the tiniest subnormal. */
5802 /* now multiply by rhs and subtract 1, at the wider precision */
5806 /* the iteration ends when the adjustment cannot affect the */
5807 /* result by >=0.5 ulp (at the requested digits), which */
5808 /* is when its value is smaller than the accumulator by */
5809 /* set->digits+1 digits (or it is zero) -- this is a looser */
5810 /* requirement than for Exp because all that happens to the */
5811 /* accumulator after this is the final rounding (but note that */
5812 /* there must also be full precision in a, or a=0). */
5822 }
5823 /* force padding if adjustment has gone to 0 before full length */
5825 }
5827 /* not done yet ... */
5830 /* lengthen the next calculation */
5837 #if DECCHECK
5838 /* just a sanity check; remove the test to show always */
5842 #endif
5844 /* Copy and round the result to res */
5854 /* [status is handled by caller] */
5858 /* ------------------------------------------------------------------ */
5859 /* decQuantizeOp -- force exponent to requested value */
5860 /* */
5861 /* This computes C = op(A, B), where op adjusts the coefficient */
5862 /* of C (by rounding or shifting) such that the exponent (-scale) */
5863 /* of C has the value B or matches the exponent of B. */
5864 /* The numerical value of C will equal A, except for the effects of */
5865 /* any rounding that occurred. */
5866 /* */
5867 /* res is C, the result. C may be A or B */
5868 /* lhs is A, the number to adjust */
5869 /* rhs is B, the requested exponent */
5870 /* set is the context */
5871 /* quant is 1 for quantize or 0 for rescale */
5872 /* status is the status accumulator (this can be called without */
5873 /* risk of control loss) */
5874 /* */
5875 /* C must have space for set->digits digits. */
5876 /* */
5877 /* Unless there is an error or the result is infinite, the exponent */
5878 /* after the operation is guaranteed to be that requested. */
5879 /* ------------------------------------------------------------------ */
5883 #if DECSUBSET
5886 #endif
5893 #if DECCHECK
5895 #endif
5898 #if DECSUBSET
5900 /* reduce operands and set lostDigits status, as needed */
5905 }
5910 }
5911 }
5912 #endif
5913 /* [following code does not require input rounding] */
5915 /* Handle special values */
5917 /* NaNs get usual processing */
5920 /* one infinity but not both is bad */
5923 /* both infinity: return lhs */
5926 }
5928 /* set requested exponent */
5931 /* Original rhs must be an integer that fits and is in range, */
5932 /* which could be from -1999999997 to +999999999, thanks to */
5933 /* subnormals */
5935 }
5937 #if DECSUBSET
5939 #endif
5948 /* the RHS has been processed, so it can be overwritten now if necessary */
5952 #if DECSUBSET
5954 #endif
5955 }
5958 /* if adjusted coefficient will definitely not fit, give up now */
5962 }
5965 /* this will decrease the length of the coefficient by adjust */
5966 /* digits, and must round as it does so */
5970 /* [note that the latter can be <1, here] */
5974 /* If just rounded a 999s case, exponent will be off by one; */
5975 /* adjust back (after checking space), if so. */
5977 /* re-check needed, e.g., for quantize(0.9999, 0.001) under */
5978 /* set->digits==3 */
5983 }
5986 }
5987 #if DECSUBSET
5989 #endif
5992 /* this will increase the length of the coefficient by -adjust */
5993 /* digits, by adding zero or more trailing zeros; this is */
5994 /* already checked for fit, above */
5996 /* if padding needed (adjust<0), add it now... */
6000 }
6004 /* Check for overflow [do not use Finalize in this case, as an */
6005 /* overflow here is a "don't fit" situation] */
6009 }
6013 }
6016 #if DECSUBSET
6019 #endif
6023 /* ------------------------------------------------------------------ */
6024 /* decCompareOp -- compare, min, or max two Numbers */
6025 /* */
6026 /* This computes C = A ? B and carries out one of four operations: */
6027 /* COMPARE -- returns the signum (as a number) giving the */
6028 /* result of a comparison unless one or both */
6029 /* operands is a NaN (in which case a NaN results) */
6030 /* COMPSIG -- as COMPARE except that a quiet NaN raises */
6031 /* Invalid operation. */
6032 /* COMPMAX -- returns the larger of the operands, using the */
6033 /* 754r maxnum operation */
6034 /* COMPMAXMAG -- ditto, comparing absolute values */
6035 /* COMPMIN -- the 754r minnum operation */
6036 /* COMPMINMAG -- ditto, comparing absolute values */
6037 /* COMTOTAL -- returns the signum (as a number) giving the */
6038 /* result of a comparison using 754r total ordering */
6039 /* */
6040 /* res is C, the result. C may be A and/or B (e.g., X=X?X) */
6041 /* lhs is A */
6042 /* rhs is B */
6043 /* set is the context */
6044 /* op is the operation flag */
6045 /* status is the usual accumulator */
6046 /* */
6047 /* C must have space for one digit for COMPARE or set->digits for */
6048 /* COMPMAX, COMPMIN, COMPMAXMAG, or COMPMINMAG. */
6049 /* ------------------------------------------------------------------ */
6050 /* The emphasis here is on speed for common cases, and avoiding */
6051 /* coefficient comparison if possible. */
6052 /* ------------------------------------------------------------------ */
6056 #if DECSUBSET
6059 #endif
6063 #if DECCHECK
6065 #endif
6068 #if DECSUBSET
6070 /* reduce operands and set lostDigits status, as needed */
6075 }
6080 }
6081 }
6082 #endif
6083 /* [following code does not require input rounding] */
6085 /* If total ordering then handle differing signs 'up front' */
6090 }
6094 }
6095 }
6097 /* handle NaNs specially; let infinities drop through */
6098 /* This assumes sNaN (even just one) leads to NaN. */
6105 /* signs are known to be the same; compute the ordering here */
6106 /* as if the signs are both positive, then invert for negatives */
6109 /* here if both NaNs */
6113 /* now it just depends on the payload */
6116 /* [Error not possible, as these are 'aligned'] */
6124 /* min or max -- 754r rules ignore single NaN */
6126 /* just one NaN; force choice to be the non-NaN operand */
6131 }
6136 }
6137 /* have numbers */
6146 /* operands are numerically equal or same NaN (and same sign, */
6147 /* tested first); if identical, leave result 0 */
6158 }
6159 }
6163 /* choose the operand for the result */
6166 /* choose according to sign then exponent (see 754r) */
6169 #if DECSUBSET
6173 }
6174 else
6175 #endif
6179 }
6183 /* [if equal, use lhs, technically identical] */
6184 }
6188 /* [ditto] */
6189 }
6191 /* here result will be non-0; reverse if looking for MIN */
6194 /* copy chosen to result, rounding if need be */
6197 }
6198 }
6199 #if DECSUBSET
6202 #endif
6206 /* ------------------------------------------------------------------ */
6207 /* decCompare -- compare two decNumbers by numerical value */
6208 /* */
6209 /* This routine compares A ? B without altering them. */
6210 /* */
6211 /* Arg1 is A, a decNumber which is not a NaN */
6212 /* Arg2 is B, a decNumber which is not a NaN */
6213 /* Arg3 is 1 for a sign-independent compare, 0 otherwise */
6214 /* */
6215 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
6216 /* (the only possible failure is an allocation error) */
6217 /* ------------------------------------------------------------------ */
6219 Flag abs) {
6228 /* RHS is non-zero */
6230 /* [here, both non-zero, result=1] */
6231 }
6240 }
6242 /* signums are the same; both are non-zero */
6247 }
6249 }
6250 /* must compare the coefficients, allowing for exponents */
6252 /* swap sides, and sign */
6257 }
6265 /* ------------------------------------------------------------------ */
6266 /* decUnitCompare -- compare two >=0 integers in Unit arrays */
6267 /* */
6268 /* This routine compares A ? B*10**E where A and B are unit arrays */
6269 /* A is a plain integer */
6270 /* B has an exponent of E (which must be non-negative) */
6271 /* */
6272 /* Arg1 is A first Unit (lsu) */
6273 /* Arg2 is A length in Units */
6274 /* Arg3 is B first Unit (lsu) */
6275 /* Arg4 is B length in Units */
6276 /* Arg5 is E (0 if the units are aligned) */
6277 /* */
6278 /* returns -1, 0, or 1 for A<B, A==B, or A>B, or BADINT if failure */
6279 /* (the only possible failure is an allocation error, which can */
6280 /* only occur if E!=0) */
6281 /* ------------------------------------------------------------------ */
6294 /* same number of units in both -- need unit-by-unit compare */
6300 }
6304 /* Unaligned. If one is >1 unit longer than the other, padded */
6305 /* approximately, then can return easily */
6309 /* Need to do a real subtract. For this, a result buffer is needed */
6310 /* even though only the sign is of interest. Its length needs */
6311 /* to be the larger of alength and padded blength, +2 */
6320 }
6321 /* Calculate units and remainder from exponent. */
6324 /* subtract [A+B*(-m)] */
6327 /* [UnitAddSub result may have leading zeros, even on zero] */
6330 /* check units of the result before freeing any storage */
6333 }
6334 /* clean up and return the result */
6339 /* ------------------------------------------------------------------ */
6340 /* decUnitAddSub -- add or subtract two >=0 integers in Unit arrays */
6341 /* */
6342 /* This routine performs the calculation: */
6343 /* */
6344 /* C=A+(B*M) */
6345 /* */
6346 /* Where M is in the range -DECDPUNMAX through +DECDPUNMAX. */
6347 /* */
6348 /* A may be shorter or longer than B. */
6349 /* */
6350 /* Leading zeros are not removed after a calculation. The result is */
6351 /* either the same length as the longer of A and B (adding any */
6352 /* shift), or one Unit longer than that (if a Unit carry occurred). */
6353 /* */
6354 /* A and B content are not altered unless C is also A or B. */
6355 /* C may be the same array as A or B, but only if no zero padding is */
6356 /* requested (that is, C may be B only if bshift==0). */
6357 /* C is filled from the lsu; only those units necessary to complete */
6358 /* the calculation are referenced. */
6359 /* */
6360 /* Arg1 is A first Unit (lsu) */
6361 /* Arg2 is A length in Units */
6362 /* Arg3 is B first Unit (lsu) */
6363 /* Arg4 is B length in Units */
6364 /* Arg5 is B shift in Units (>=0; pads with 0 units if positive) */
6365 /* Arg6 is C first Unit (lsu) */
6366 /* Arg7 is M, the multiplier */
6367 /* */
6368 /* returns the count of Units written to C, which will be non-zero */
6369 /* and negated if the result is negative. That is, the sign of the */
6370 /* returned Int is the sign of the result (positive for zero) and */
6371 /* the absolute value of the Int is the count of Units. */
6372 /* */
6373 /* It is the caller's responsibility to make sure that C size is */
6374 /* safe, allowing space if necessary for a one-Unit carry. */
6375 /* */
6376 /* This routine is severely performance-critical; *any* change here */
6377 /* must be measured (timed) to assure no performance degradation. */
6378 /* In particular, trickery here tends to be counter-productive, as */
6379 /* increased complexity of code hurts register optimizations on */
6380 /* register-poor architectures. Avoiding divisions is nearly */
6381 /* always a Good Idea, however. */
6382 /* */
6383 /* Special thanks to Rick McGuire (IBM Cambridge, MA) and Dave Clark */
6384 /* (IBM Warwick, UK) for some of the ideas used in this routine. */
6385 /* ------------------------------------------------------------------ */
6397 #endif
6399 #if DECTRACE
6402 #endif
6408 /* if in place [common], skip copy unless there's a gap [rare] */
6412 }
6416 }
6417 }
6422 }
6424 /* For speed, do the addition as two loops; the first where both A */
6425 /* and B contribute, and the second (if necessary) where only one or */
6426 /* other of the numbers contribute. */
6427 /* Carry handling is the same (i.e., duplicated) in each case. */
6430 a++;
6433 /* here carry is new Unit of digits; it could be +ve or -ve */
6438 }
6445 carry++;
6448 }
6449 /* negative case */
6455 carry++;
6457 #elif DECDPUN==3
6463 carry++;
6466 }
6467 /* negative case */
6473 carry++;
6475 #elif DECDPUN<=2
6476 /* Can use QUOT10 as carry <= 4 digits */
6482 }
6483 /* negative case */
6488 #else
6489 /* remainder operator is undefined if negative, so must test */
6494 }
6499 }
6500 /* negative case */
6504 #endif
6507 /* now may have one or other to complete */
6508 /* [pretest to avoid loop setup/shutdown] */
6512 a++;
6513 }
6516 b++;
6517 }
6518 /* here carry is new Unit of digits; it could be +ve or -ve and */
6519 /* magnitude up to DECDPUNMAX squared */
6524 }
6525 /* result for this unit is negative or >DECDPUNMAX */
6532 carry++;
6535 }
6536 /* negative case */
6542 carry++;
6544 #elif DECDPUN==3
6550 carry++;
6553 }
6554 /* negative case */
6560 carry++;
6562 #elif DECDPUN<=2
6568 }
6569 /* negative case */
6574 #else
6579 }
6580 /* remainder operator is undefined if negative, so must test */
6585 }
6586 /* negative case */
6590 #endif
6593 /* OK, all A and B processed; might still have carry or borrow */
6594 /* return number of Units in the result, negated if a borrow */
6600 }
6601 /* -ve carry: it's a borrow; complement needed */
6608 }
6612 }
6613 }
6614 /* add an extra unit iff it would be non-zero */
6615 #if DECTRACE
6617 #endif
6621 }
6625 /* ------------------------------------------------------------------ */
6626 /* decTrim -- trim trailing zeros or normalize */
6627 /* */
6628 /* dn is the number to trim or normalize */
6629 /* set is the context to use to check for clamp */
6630 /* all is 1 to remove all trailing zeros, 0 for just fraction ones */
6631 /* dropped returns the number of discarded trailing zeros */
6632 /* returns dn */
6633 /* */
6634 /* If clamp is set in the context then the number of zeros trimmed */
6635 /* may be limited if the exponent is high. */
6636 /* All fields are updated as required. This is a utility operation, */
6637 /* so special values are unchanged and no error is possible. */
6638 /* ------------------------------------------------------------------ */
6645 #if DECCHECK
6647 #endif
6655 }
6657 /* have a finite number which is even */
6662 /* slice by powers */
6663 #if DECDPUN<=4
6666 #else
6668 #endif
6669 /* have a trailing 0 */
6671 /* [if exp>0 then all trailing 0s are significant for trim] */
6675 }
6676 }
6679 up++;
6681 }
6685 /* may need to limit drop if clamping */
6690 }
6692 /* effect the drop */
6700 /* ------------------------------------------------------------------ */
6701 /* decReverse -- reverse a Unit array in place */
6702 /* */
6703 /* ulo is the start of the array */
6704 /* uhi is the end of the array (highest Unit to include) */
6705 /* */
6706 /* The units ulo through uhi are reversed in place (if the number */
6707 /* of units is odd, the middle one is untouched). Note that the */
6708 /* digit(s) in each unit are unaffected. */
6709 /* ------------------------------------------------------------------ */
6711 Unit temp;
6716 }
6720 /* ------------------------------------------------------------------ */
6721 /* decShiftToMost -- shift digits in array towards most significant */
6722 /* */
6723 /* uar is the array */
6724 /* digits is the count of digits in use in the array */
6725 /* shift is the number of zeros to pad with (least significant); */
6726 /* it must be zero or positive */
6727 /* */
6728 /* returns the new length of the integer in the array, in digits */
6729 /* */
6730 /* No overflow is permitted (that is, the uar array must be known to */
6731 /* be large enough to hold the result, after shifting). */
6732 /* ------------------------------------------------------------------ */
6742 }
6750 }
6754 /* split the source Unit and accumulate remainder for next */
6755 #if DECDPUN<=4
6759 #else
6762 #endif
6765 }
6768 /* propagate any partial unit to one below and clear the rest */
6772 }
6776 /* ------------------------------------------------------------------ */
6777 /* decShiftToLeast -- shift digits in array towards least significant */
6778 /* */
6779 /* uar is the array */
6780 /* units is length of the array, in units */
6781 /* shift is the number of digits to remove from the lsu end; it */
6782 /* must be zero or positive and <= than units*DECDPUN. */
6783 /* */
6784 /* returns the new length of the integer in the array, in units */
6785 /* */
6786 /* Removed digits are discarded (lost). Units not required to hold */
6787 /* the final result are unchanged. */
6788 /* ------------------------------------------------------------------ */
6798 }
6806 }
6808 /* messier */
6811 #if DECDPUN<=4
6813 #else
6815 #endif
6820 up++;
6822 #if DECDPUN<=4
6825 #else
6828 #endif
6832 }
6836 #if DECSUBSET
6837 /* ------------------------------------------------------------------ */
6838 /* decRoundOperand -- round an operand [used for subset only] */
6839 /* */
6840 /* dn is the number to round (dn->digits is > set->digits) */
6841 /* set is the relevant context */
6842 /* status is the status accumulator */
6843 /* */
6844 /* returns an allocated decNumber with the rounded result. */
6845 /* */
6846 /* lostDigits and other status may be set by this. */
6847 /* */
6848 /* Since the input is an operand, it must not be modified. */
6849 /* Instead, return an allocated decNumber, rounded as required. */
6850 /* It is the caller's responsibility to free the allocated storage. */
6851 /* */
6852 /* If no storage is available then the result cannot be used, so NULL */
6853 /* is returned. */
6854 /* ------------------------------------------------------------------ */
6861 /* Allocate storage for the returned decNumber, big enough for the */
6862 /* length specified by the context */
6868 }
6872 /* If that set Inexact then "lost digits" is raised... */
6877 #endif
6879 /* ------------------------------------------------------------------ */
6880 /* decCopyFit -- copy a number, truncating the coefficient if needed */
6881 /* */
6882 /* dest is the target decNumber */
6883 /* src is the source decNumber */
6884 /* set is the context [used for length (digits) and rounding mode] */
6885 /* residue is the residue accumulator */
6886 /* status contains the current status to be updated */
6887 /* */
6888 /* (dest==src is allowed and will be a no-op if fits) */
6889 /* All fields are updated as required. */
6890 /* ------------------------------------------------------------------ */
6898 /* ------------------------------------------------------------------ */
6899 /* decSetCoeff -- set the coefficient of a number */
6900 /* */
6901 /* dn is the number whose coefficient array is to be set. */
6902 /* It must have space for set->digits digits */
6903 /* set is the context [for size] */
6904 /* lsu -> lsu of the source coefficient [may be dn->lsu] */
6905 /* len is digits in the source coefficient [may be dn->digits] */
6906 /* residue is the residue accumulator. This has values as in */
6907 /* decApplyRound, and will be unchanged unless the */
6908 /* target size is less than len. In this case, the */
6909 /* coefficient is truncated and the residue is updated to */
6910 /* reflect the previous residue and the dropped digits. */
6911 /* status is the status accumulator, as usual */
6912 /* */
6913 /* The coefficient may already be in the number, or it can be an */
6914 /* external intermediate array. If it is in the number, lsu must == */
6915 /* dn->lsu and len must == dn->digits. */
6916 /* */
6917 /* Note that the coefficient length (len) may be < set->digits, and */
6918 /* in this case this merely copies the coefficient (or is a no-op */
6919 /* if dn->lsu==lsu). */
6920 /* */
6921 /* Note also that (only internally, from decQuantizeOp and */
6922 /* decSetSubnormal) the value of set->digits may be less than one, */
6923 /* indicating a round to left. This routine handles that case */
6924 /* correctly; caller ensures space. */
6925 /* */
6926 /* dn->digits, dn->lsu (and as required), and dn->exponent are */
6927 /* updated as necessary. dn->bits (sign) is unchanged. */
6928 /* */
6929 /* DEC_Rounded status is set if any digits are discarded. */
6930 /* DEC_Inexact status is set if any non-zero digits are discarded, or */
6931 /* incoming residue was non-0 (implies rounded) */
6932 /* ------------------------------------------------------------------ */
6933 /* mapping array: maps 0-9 to canonical residues, so that a residue */
6934 /* can be adjusted in the range [-1, +1] and achieve correct rounding */
6935 /* 0 1 2 3 4 5 6 7 8 9 */
6944 #if DECDPUN<=4
6946 #endif
6951 /* copy the coefficient array to the result number; no shift needed */
6957 }
6958 /* dn->exponent and residue are unchanged, record any inexactitude */
6961 }
6963 /* some digits must be discarded ... */
6969 /* guard digit is 0 */
6970 /* residue is all the number [NB could be all 0s] */
6976 }
6977 }
6984 /* partial discard [most common case] */
6985 /* here, at least the first (most significant) discarded digit exists */
6987 /* spin up the number, noting residue during the spin, until get to */
6988 /* the Unit with the first discarded digit. When reach it, extract */
6989 /* it and remember its position */
6997 /* here up -> Unit with first discarded digit */
7001 /* set residue directly */
7005 }
7008 }
7012 }
7017 /* on unit boundary, so shift-down copy loop is simple */
7020 }
7028 #if DECDPUN<=4
7031 #else
7034 #endif
7036 }
7037 /* discard digit is now at bottom of quot */
7038 #if DECDPUN<=4
7040 /* Vowels algorithm here not a win (9 instructions) */
7043 #else
7046 #endif
7047 /* here, discard1 is the guard digit, and residue is everything */
7048 /* else [use mapping array to accumulate residue safely] */
7051 /* here: up -> Unit of the array with bottom digit */
7052 /* cut is the division point for each Unit */
7053 /* quot holds the uncut high-order digits for the current unit */
7057 }
7061 /* shift-copy the coefficient array to the result number */
7066 up++;
7068 #if DECDPUN<=4
7071 #else
7074 #endif
7086 /* ------------------------------------------------------------------ */
7087 /* decApplyRound -- apply pending rounding to a number */
7088 /* */
7089 /* dn is the number, with space for set->digits digits */
7090 /* set is the context [for size and rounding mode] */
7091 /* residue indicates pending rounding, being any accumulated */
7092 /* guard and sticky information. It may be: */
7093 /* 6-9: rounding digit is >5 */
7094 /* 5: rounding digit is exactly half-way */
7095 /* 1-4: rounding digit is <5 and >0 */
7096 /* 0: the coefficient is exact */
7097 /* -1: as 1, but the hidden digits are subtractive, that */
7098 /* is, of the opposite sign to dn. In this case the */
7099 /* coefficient must be non-0. This case occurs when */
7100 /* subtracting a small number (which can be reduced to */
7101 /* a sticky bit); see decAddOp. */
7102 /* status is the status accumulator, as usual */
7103 /* */
7104 /* This routine applies rounding while keeping the length of the */
7105 /* coefficient constant. The exponent and status are unchanged */
7106 /* except if: */
7107 /* */
7108 /* -- the coefficient was increased and is all nines (in which */
7109 /* case Overflow could occur, and is handled directly here so */
7110 /* the caller does not need to re-test for overflow) */
7111 /* */
7112 /* -- the coefficient was decreased and becomes all nines (in which */
7113 /* case Underflow could occur, and is also handled directly). */
7114 /* */
7115 /* All fields in dn are updated as required. */
7116 /* */
7117 /* ------------------------------------------------------------------ */
7121 /* -1 if coefficient needs to be decremented */
7127 /* now decide whether, and how, to round, depending on mode */
7130 /* This is the same as DEC_ROUND_DOWN unless there is a */
7131 /* positive residue and the lsd of dn is 0 or 5, in which case */
7132 /* it is bumped; when residue is <0, the number is therefore */
7133 /* bumped down unless the final digit was 1 or 6 (in which */
7134 /* case it is bumped down and then up -- a no-op) */
7138 /* [bump==1 could be applied directly; use common path for clarity] */
7142 /* no change, except if negative residue */
7153 /* 0.5 goes up iff [new] lsd is odd */
7155 }
7167 /* same as _UP for positive numbers, and as _DOWN for negatives */
7168 /* [negative residue cannot occur on 0] */
7171 }
7174 }
7178 /* same as _UP for negative numbers, and as _DOWN for positive */
7179 /* [negative residue cannot occur on 0] */
7182 }
7185 }
7190 #if DECTRACE || (DECCHECK && DECVERB)
7192 #endif
7196 /* now bump the number, up or down, if need be */
7199 /* Simply use decUnitAddSub unless bumping up and the number is */
7200 /* all nines. In this special case set to 100... explicitly */
7201 /* and adjust the exponent by one (as otherwise could overflow */
7202 /* the array) */
7203 /* Similarly handle all-nines result if bumping down. */
7209 /* this is the last Unit (the msu) */
7211 /* here if it, too, is all nines */
7215 /* [which, very rarely, could cause Overflow...] */
7218 }
7220 }
7221 /* a full unit to check, with more to come */
7227 /* here checking for a pre-bump of 1000... (leading 1, all */
7228 /* other digits zero) */
7233 /* this is the last Unit (the msu) */
7235 /* here if have the 1000... case */
7238 /* others all to all-nines, too */
7242 /* iff the number was at the subnormal boundary (exponent=etiny) */
7243 /* then the exponent is now out of range, so it will in fact get */
7244 /* clamped to etiny and the final 9 dropped. */
7245 /* printf(">> emin=%d exp=%d sdig=%d\n", set->emin, */
7246 /* dn->exponent, set->digits); */
7252 }
7255 }
7257 }
7259 /* a full unit to check, with more to come */
7266 /* Actual bump needed. Do it. */
7270 #if DECSUBSET
7271 /* ------------------------------------------------------------------ */
7272 /* decFinish -- finish processing a number */
7273 /* */
7274 /* dn is the number */
7275 /* set is the context */
7276 /* residue is the rounding accumulator (as in decApplyRound) */
7277 /* status is the accumulator */
7278 /* */
7279 /* This finishes off the current number by: */
7280 /* 1. If not extended: */
7281 /* a. Converting a zero result to clean '0' */
7282 /* b. Reducing positive exponents to 0, if would fit in digits */
7283 /* 2. Checking for overflow and subnormals (always) */
7284 /* Note this is just Finalize when no subset arithmetic. */
7285 /* All fields are updated as required. */
7286 /* ------------------------------------------------------------------ */
7294 }
7296 /* >0; reduce to integer if possible */
7300 }
7301 }
7306 #endif
7308 /* ------------------------------------------------------------------ */
7309 /* decFinalize -- final check, clamp, and round of a number */
7310 /* */
7311 /* dn is the number */
7312 /* set is the context */
7313 /* residue is the rounding accumulator (as in decApplyRound) */
7314 /* status is the status accumulator */
7315 /* */
7316 /* This finishes off the current number by checking for subnormal */
7317 /* results, applying any pending rounding, checking for overflow, */
7318 /* and applying any clamping. */
7319 /* Underflow and overflow conditions are raised as appropriate. */
7320 /* All fields are updated as required. */
7321 /* ------------------------------------------------------------------ */
7327 /* Must be careful, here, when checking the exponent as the */
7328 /* adjusted exponent could overflow 31 bits [because it may already */
7329 /* be up to twice the expected]. */
7331 /* First test for subnormal. This must be done before any final */
7332 /* round as the result could be rounded to Nmin or 0. */
7334 Int comp;
7335 decNumber nmin;
7336 /* A very nasty case here is dn == Nmin and residue<0 */
7338 /* Go handle subnormals; this will apply round if needed. */
7341 }
7342 /* Equals case: only subnormal if dn=Nmin and negative residue */
7350 }
7355 }
7356 }
7358 /* now apply any pending round (this could raise overflow). */
7361 /* Check for overflow [redundant in the 'rare' case] or clamp */
7365 /* here when might have an overflow or clamp to do */
7369 }
7370 /* here when the result is normal but in clamp range */
7373 /* here when need to apply the IEEE exponent clamp (fold-down) */
7376 /* shift coefficient (if non-zero) */
7379 }
7385 /* ------------------------------------------------------------------ */
7386 /* decSetOverflow -- set number to proper overflow value */
7387 /* */
7388 /* dn is the number (used for sign [only] and result) */
7389 /* set is the context [used for the rounding mode, etc.] */
7390 /* status contains the current status to be updated */
7391 /* */
7392 /* This sets the sign of a number and sets its value to either */
7393 /* Infinity or the maximum finite value, depending on the sign of */
7394 /* dn and the rounding mode, following IEEE 854 rules. */
7395 /* ------------------------------------------------------------------ */
7406 }
7408 }
7425 }
7429 }
7434 /* ------------------------------------------------------------------ */
7435 /* decSetMaxValue -- set number to +Nmax (maximum normal value) */
7436 /* */
7437 /* dn is the number to set */
7438 /* set is the context [used for digits and emax] */
7439 /* */
7440 /* This sets the number to the maximum positive value. */
7441 /* ------------------------------------------------------------------ */
7446 /* fill in all nines to set maximum value */
7452 }
7459 /* ------------------------------------------------------------------ */
7460 /* decSetSubnormal -- process value whose exponent is <Emin */
7461 /* */
7462 /* dn is the number (used as input as well as output; it may have */
7463 /* an allowed subnormal value, which may need to be rounded) */
7464 /* set is the context [used for the rounding mode] */
7465 /* residue is any pending residue */
7466 /* status contains the current status to be updated */
7467 /* */
7468 /* If subset mode, set result to zero and set Underflow flags. */
7469 /* */
7470 /* Value may be zero with a low exponent; this does not set Subnormal */
7471 /* but the exponent will be clamped to Etiny. */
7472 /* */
7473 /* Otherwise ensure exponent is not out of range, and round as */
7474 /* necessary. Underflow is set if the result is Inexact. */
7475 /* ------------------------------------------------------------------ */
7481 #if DECSUBSET
7482 /* simple set to zero and 'hard underflow' for subset */
7485 /* always full overflow */
7488 }
7489 #endif
7491 /* Full arithmetic -- allow subnormals, rounded to minimum exponent */
7492 /* (Etiny) if needed */
7496 /* residue can never be non-zero here */
7497 #if DECCHECK
7501 }
7502 #endif
7506 }
7508 }
7513 /* residue can never be non-zero here, except in the Nmin-residue */
7514 /* case (which is a subnormal result), so can take fast-path here */
7515 /* it may already be inexact (from setting the coefficient) */
7518 }
7520 /* adjust>0, so need to rescale the result so exponent becomes Etiny */
7521 /* [this code is similar to that in rescale] */
7525 /* [note that the latter can be <1, here, similar to Rescale case] */
7529 /* Use 754R/854 default rule: Underflow is set iff Inexact */
7530 /* [independent of whether trapped] */
7533 /* if rounded up a 999s case, exponent will be off by one; adjust */
7534 /* back if so [it will fit, because it was shortened earlier] */
7538 }
7540 /* if rounded to zero, it is by definition clamped... */
7544 /* ------------------------------------------------------------------ */
7545 /* decCheckMath - check entry conditions for a math function */
7546 /* */
7547 /* This checks the context and the operand */
7548 /* */
7549 /* rhs is the operand to check */
7550 /* set is the context to check */
7551 /* status is unchanged if both are good */
7552 /* */
7553 /* returns non-zero if status is changed, 0 otherwise */
7554 /* */
7555 /* Restrictions enforced: */
7556 /* */
7557 /* digits, emax, and -emin in the context must be less than */
7558 /* DEC_MAX_MATH (999999), and A must be within these bounds if */
7559 /* non-zero. Invalid_operation is set in the status if a */
7560 /* restriction is violated. */
7561 /* ------------------------------------------------------------------ */
7575 /* ------------------------------------------------------------------ */
7576 /* decGetInt -- get integer from a number */
7577 /* */
7578 /* dn is the number [which will not be altered] */
7579 /* */
7580 /* returns one of: */
7581 /* BADINT if there is a non-zero fraction */
7582 /* the converted integer */
7583 /* BIGEVEN if the integer is even and magnitude > 2*10**9 */
7584 /* BIGODD if the integer is odd and magnitude > 2*10**9 */
7585 /* */
7586 /* This checks and gets a whole number from the input decNumber. */
7587 /* The sign can be determined from dn by the caller when BIGEVEN or */
7588 /* BIGODD is returned. */
7589 /* ------------------------------------------------------------------ */
7597 /* The number must be an integer that fits in 10 digits */
7598 /* Assert, here, that 10 is enough for any rescale Etiny */
7599 #if DEC_MAX_EMAX > 999999999
7600 #error GetInt may need updating [for Emax]
7601 #endif
7602 #if DEC_MIN_EMIN < -999999999
7603 #error GetInt may need updating [for Emin]
7604 #endif
7610 /* no fractional part [usual]; allow for positive exponent */
7612 }
7615 /* spin up whole units until reach the Unit with the unit digit */
7619 }
7623 /* slice off fraction digits and check for non-zero */
7624 #if DECDPUN<=4
7627 #else
7630 #endif
7632 /* it looks good */
7635 }
7636 }
7637 /* now it's known there's no fractional part */
7639 /* tricky code now, to accumulate up to 9.3 digits */
7644 /* collect any remaining unit(s) */
7648 }
7651 /* [that test also disallows the BADINT result case] */
7655 }
7656 }
7661 }
7667 /* ------------------------------------------------------------------ */
7668 /* decDecap -- decapitate the coefficient of a number */
7669 /* */
7670 /* dn is the number to be decapitated */
7671 /* drop is the number of digits to be removed from the left of dn; */
7672 /* this must be <= dn->digits (if equal, the coefficient is */
7673 /* set to 0) */
7674 /* */
7675 /* Returns dn; dn->digits will be <= the initial digits less drop */
7676 /* (after removing drop digits there may be leading zero digits */
7677 /* which will also be removed). Only dn->lsu and dn->digits change. */
7678 /* ------------------------------------------------------------------ */
7683 #if DECCHECK
7687 #endif
7691 }
7695 /* that may have left leading zero digits, so do a proper count... */
7700 /* ------------------------------------------------------------------ */
7701 /* decBiStr -- compare string with pairwise options */
7702 /* */
7703 /* targ is the string to compare */
7704 /* str1 is one of the strings to compare against (length may be 0) */
7705 /* str2 is the other; it must be the same length as str1 */
7706 /* */
7707 /* returns 1 if strings compare equal, (that is, it is the same */
7708 /* length as str1 and str2, and each character of targ is in either */
7709 /* str1 or str2 in the corresponding position), or 0 otherwise */
7710 /* */
7711 /* This is used for generic caseless compare, including the awkward */
7712 /* case of the Turkish dotted and dotless Is. Use as (for example): */
7713 /* if (decBiStr(test, "mike", "MIKE")) ... */
7714 /* ------------------------------------------------------------------ */
7718 /* *targ has a match in one (or both, if terminator) */
7724 /* ------------------------------------------------------------------ */
7725 /* decNaNs -- handle NaN operand or operands */
7726 /* */
7727 /* res is the result number */
7728 /* lhs is the first operand */
7729 /* rhs is the second operand, or NULL if none */
7730 /* context is used to limit payload length */
7731 /* status contains the current status */
7732 /* returns res in case convenient */
7733 /* */
7734 /* Called when one or both operands is a NaN, and propagates the */
7735 /* appropriate result to res. When an sNaN is found, it is changed */
7736 /* to a qNaN and Invalid operation is set. */
7737 /* ------------------------------------------------------------------ */
7741 /* This decision tree ends up with LHS being the source pointer, */
7742 /* and status updated if need be */
7749 }
7753 /* propagate the payload */
7758 /* copy safe number of units, then decapitate */
7763 /* maybe still too long */
7765 }
7770 /* [coefficient was copied/decapitated] */
7774 /* ------------------------------------------------------------------ */
7775 /* decStatus -- apply non-zero status */
7776 /* */
7777 /* dn is the number to set if error */
7778 /* status contains the current status (not yet in context) */
7779 /* set is the context */
7780 /* */
7781 /* If the status is an error status, the number is set to a NaN, */
7782 /* unless the error was an overflow, divide-by-zero, or underflow, */
7783 /* in which case the number will have already been set. */
7784 /* */
7785 /* The context status is then updated with the new status. Note that */
7786 /* this may raise a signal, so control may never return from this */
7787 /* routine (hence resources must be recovered before it is called). */
7788 /* ------------------------------------------------------------------ */
7791 /* if cause was an sNaN, clear and propagate [NaN is already set up] */
7796 }
7797 }
7802 /* ------------------------------------------------------------------ */
7803 /* decGetDigits -- count digits in a Units array */
7804 /* */
7805 /* uar is the Unit array holding the number (this is often an */
7806 /* accumulator of some sort) */
7807 /* len is the length of the array in units [>=1] */
7808 /* */
7809 /* returns the number of (significant) digits in the array */
7810 /* */
7811 /* All leading zeros are excluded, except the last if the array has */
7812 /* only zero Units. */
7813 /* ------------------------------------------------------------------ */
7814 /* This may be called twice during some operations. */
7818 #if DECDPUN>4
7820 #endif
7821 /* (at least 1 in final msu) */
7822 #if DECCHECK
7824 #endif
7831 /* found the first (most significant) non-zero Unit */
7834 digits++;
7837 digits++;
7840 digits++;
7843 #endif
7844 #endif
7845 #endif
7846 #endif
7852 #if DECTRACE | DECCHECK
7853 /* ------------------------------------------------------------------ */
7854 /* decNumberShow -- display a number [debug aid] */
7855 /* dn is the number to show */
7856 /* */
7857 /* Shows: sign, exponent, coefficient (msu first), digits */
7858 /* or: sign, special-value */
7859 /* ------------------------------------------------------------------ */
7860 /* this is public so other modules can use it */
7876 }
7877 /* if coefficient and exponent are 0, no more to do */
7881 /* drop through to report other information */
7883 }
7885 /* now carefully display the coefficient */
7901 }
7904 #endif
7906 #if DECTRACE || DECCHECK
7907 /* ------------------------------------------------------------------ */
7908 /* decDumpAr -- display a unit array [debug/check aid] */
7909 /* name is a single-character tag name */
7910 /* ar is the array to display */
7911 /* len is the length of the array in Units */
7912 /* ------------------------------------------------------------------ */
7914 Int i;
7916 #if DECDPUN==9
7918 #elif DECDPUN==8
7920 #elif DECDPUN==7
7922 #elif DECDPUN==6
7924 #elif DECDPUN==5
7926 #elif DECDPUN==4
7928 #elif DECDPUN==3
7930 #elif DECDPUN==2
7932 #else
7934 #endif
7939 }
7942 #endif
7944 #if DECCHECK
7945 /* ------------------------------------------------------------------ */
7946 /* decCheckOperands -- check operand(s) to a routine */
7947 /* res is the result structure (not checked; it will be set to */
7948 /* quiet NaN if error found (and it is not NULL)) */
7949 /* lhs is the first operand (may be DECUNRESU) */
7950 /* rhs is the second (may be DECUNUSED) */
7951 /* set is the context (may be DECUNCONT) */
7952 /* returns 0 if both operands, and the context are clean, or 1 */
7953 /* otherwise (in which case the context will show an error, */
7954 /* unless NULL). Note that res is not cleaned; caller should */
7955 /* handle this so res=NULL case is safe. */
7956 /* The caller is expected to abandon immediately if 1 is returned. */
7957 /* ------------------------------------------------------------------ */
7962 #if DECTRACE || DECVERB
7964 #endif
7970 #if DECTRACE || DECVERB
7973 #endif
7974 }
7978 #if DECTRACE
7979 /* this one not DECVERB as standard tests include NULL */
7981 #endif
7982 }
7985 }
7991 }
7992 }
7996 /* ------------------------------------------------------------------ */
7997 /* decCheckNumber -- check a number */
7998 /* dn is the number to check */
7999 /* returns 0 if the number is clean, or 1 otherwise */
8000 /* */
8001 /* The number is considered valid if it could be a result from some */
8002 /* operation in some valid context. */
8003 /* ------------------------------------------------------------------ */
8011 #if DECTRACE
8012 /* this one not DECVERB as standard tests include NULL */
8014 #endif
8017 /* check special values */
8020 #if DECTRACE || DECVERB
8023 #endif
8026 /* 2003.09.08: NaNs may now have coefficients, so next tests Inf only */
8029 #if DECTRACE || DECVERB
8031 #endif
8034 #if DECTRACE || DECVERB
8036 #endif
8040 /* 2002.12.26: negative NaNs can now appear through proposed IEEE */
8041 /* concrete formats (decimal64, etc.). */
8043 }
8045 /* check the coefficient */
8047 #if DECTRACE || DECVERB
8049 #endif
8059 #if DECTRACE || DECVERB
8062 #endif
8064 }
8066 #if DECTRACE || DECVERB
8069 #endif
8072 }
8074 /* check the exponent. Note that input operands can have exponents */
8075 /* which are out of the set->emin/set->emax and set->digits range */
8076 /* (just as they can have more digits than set->digits). */
8082 #if DECTRACE || DECVERB
8085 #endif
8088 #if DECTRACE || DECVERB
8091 #endif
8097 /* ------------------------------------------------------------------ */
8098 /* decCheckInexact -- check a normal finite inexact result has digits */
8099 /* dn is the number to check */
8100 /* set is the context (for status and precision) */
8101 /* sets Invalid operation, etc., if some digits are missing */
8102 /* [this check is not made for DECSUBSET compilation or when */
8103 /* subnormal is not set] */
8104 /* ------------------------------------------------------------------ */
8106 #if !DECSUBSET && DECEXTFLAG
8109 #if DECTRACE || DECVERB
8113 #endif
8115 }
8116 #else
8117 /* next is a noop for quiet compiler */
8119 #endif
8122 #endif
8124 #if DECALLOC
8125 #undef malloc
8126 #undef free
8127 /* ------------------------------------------------------------------ */
8128 /* decMalloc -- accountable allocation routine */
8129 /* n is the number of bytes to allocate */
8130 /* */
8131 /* Semantics is the same as the stdlib malloc routine, but bytes */
8132 /* allocated are accounted for globally, and corruption fences are */
8133 /* added before and after the 'actual' storage. */
8134 /* ------------------------------------------------------------------ */
8135 /* This routine allocates storage with an extra twelve bytes; 8 are */
8136 /* at the start and hold: */
8137 /* 0-3 the original length requested */
8138 /* 4-7 buffer corruption detection fence (DECFENCE, x4) */
8139 /* The 4 bytes at the end also hold a corruption fence (DECFENCE, x4) */
8140 /* ------------------------------------------------------------------ */
8153 /* printf(" alloc ++ dAB: %ld (%d)\n", decAllocBytes, n); */
8159 /* ------------------------------------------------------------------ */
8160 /* decFree -- accountable free routine */
8161 /* alloc is the storage to free */
8162 /* */
8163 /* Semantics is the same as the stdlib malloc routine, except that */
8164 /* the global storage accounting is updated and the fences are */
8165 /* checked to ensure that no routine has written 'out of bounds'. */
8166 /* ------------------------------------------------------------------ */
8167 /* This routine first checks that the fences have not been corrupted. */
8168 /* It then frees the storage using the 'truw' storage address (that */
8169 /* is, offset by 8). */
8170 /* ------------------------------------------------------------------ */
8188 /* printf(" free -- dAB: %d (%d)\n", decAllocBytes, -n); */
8190 #define malloc(a) decMalloc(a)
8191 #define free(a) decFree(a)
8192 #endif