github.com/matrixorigin/matrixone@v0.7.0/cgo/external/decNumber/decCommon.c (about)

     1  /* ------------------------------------------------------------------ */
     2  /* decCommon.c -- common code for all three fixed-size types          */
     3  /* ------------------------------------------------------------------ */
     4  /* Copyright (c) IBM Corporation, 2000, 2010.  All rights reserved.   */
     5  /*                                                                    */
     6  /* This software is made available under the terms of the             */
     7  /* ICU License -- ICU 1.8.1 and later.                                */
     8  /*                                                                    */
     9  /* The description and User's Guide ("The decNumber C Library") for   */
    10  /* this software is included in the package as decNumber.pdf.  This   */
    11  /* document is also available in HTML, together with specifications,  */
    12  /* testcases, and Web links, on the General Decimal Arithmetic page.  */
    13  /*                                                                    */
    14  /* Please send comments, suggestions, and corrections to the author:  */
    15  /*   mfc@uk.ibm.com                                                   */
    16  /*   Mike Cowlishaw, IBM Fellow                                       */
    17  /*   IBM UK, PO Box 31, Birmingham Road, Warwick CV34 5JL, UK         */
    18  /* ------------------------------------------------------------------ */
    19  /* This module comprises code that is shared between all the formats  */
    20  /* (decSingle, decDouble, and decQuad); it includes set and extract   */
    21  /* of format components, widening, narrowing, and string conversions. */
    22  /*                                                                    */
    23  /* Unlike decNumber, parameterization takes place at compile time     */
    24  /* rather than at runtime.  The parameters are set in the decDouble.c */
    25  /* (etc.) files, which then include this one to produce the compiled  */
    26  /* code.  The functions here, therefore, are code shared between      */
    27  /* multiple formats.                                                  */
    28  /* ------------------------------------------------------------------ */
    29  // Names here refer to decFloat rather than to decDouble, etc., and
    30  // the functions are in strict alphabetical order.
    31  // Constants, tables, and debug function(s) are included only for QUAD
    32  // (which will always be compiled if DOUBLE or SINGLE are used).
    33  //
    34  // Whenever a decContext is used, only the status may be set (using
    35  // OR) or the rounding mode read; all other fields are ignored and
    36  // untouched.
    37  
    38  // names for simpler testing and default context
    39  #if DECPMAX==7
    40    #define SINGLE     1
    41    #define DOUBLE     0
    42    #define QUAD       0
    43    #define DEFCONTEXT DEC_INIT_DECIMAL32
    44  #elif DECPMAX==16
    45    #define SINGLE     0
    46    #define DOUBLE     1
    47    #define QUAD       0
    48    #define DEFCONTEXT DEC_INIT_DECIMAL64
    49  #elif DECPMAX==34
    50    #define SINGLE     0
    51    #define DOUBLE     0
    52    #define QUAD       1
    53    #define DEFCONTEXT DEC_INIT_DECIMAL128
    54  #else
    55    #error Unexpected DECPMAX value
    56  #endif
    57  
    58  /* Assertions */
    59  
    60  #if DECPMAX!=7 && DECPMAX!=16 && DECPMAX!=34
    61    #error Unexpected Pmax (DECPMAX) value for this module
    62  #endif
    63  
    64  // Assert facts about digit characters, etc.
    65  #if ('9'&0x0f)!=9
    66    #error This module assumes characters are of the form 0b....nnnn
    67    // where .... are don't care 4 bits and nnnn is 0000 through 1001
    68  #endif
    69  #if ('9'&0xf0)==('.'&0xf0)
    70    #error This module assumes '.' has a different mask than a digit
    71  #endif
    72  
    73  // Assert ToString lay-out conditions
    74  #if DECSTRING<DECPMAX+9
    75    #error ToString needs at least 8 characters for lead-in and dot
    76  #endif
    77  #if DECPMAX+DECEMAXD+5 > DECSTRING
    78    #error Exponent form can be too long for ToString to lay out safely
    79  #endif
    80  #if DECEMAXD > 4
    81    #error Exponent form is too long for ToString to lay out
    82    // Note: code for up to 9 digits exists in archives [decOct]
    83  #endif
    84  
    85  /* Private functions used here and possibly in decBasic.c, etc. */
    86  static decFloat * decFinalize(decFloat *, bcdnum *, decContext *);
    87  static Flag decBiStr(const char *, const char *, const char *);
    88  
    89  /* Macros and private tables; those which are not format-dependent    */
    90  /* are only included if decQuad is being built.                       */
    91  
    92  /* ------------------------------------------------------------------ */
    93  /* Combination field lookup tables (uInts to save measurable work)    */
    94  /*                                                                    */
    95  /*   DECCOMBEXP  - 2 most-significant-bits of exponent (00, 01, or    */
    96  /*                 10), shifted left for format, or DECFLOAT_Inf/NaN  */
    97  /*   DECCOMBWEXP - The same, for the next-wider format (unless QUAD)  */
    98  /*   DECCOMBMSD  - 4-bit most-significant-digit                       */
    99  /*                 [0 if the index is a special (Infinity or NaN)]    */
   100  /*   DECCOMBFROM - 5-bit combination field from EXP top bits and MSD  */
   101  /*                 (placed in uInt so no shift is needed)             */
   102  /*                                                                    */
   103  /* DECCOMBEXP, DECCOMBWEXP, and DECCOMBMSD are indexed by the sign    */
   104  /*   and 5-bit combination field (0-63, the second half of the table  */
   105  /*   identical to the first half)                                     */
   106  /* DECCOMBFROM is indexed by expTopTwoBits*16 + msd                   */
   107  /*                                                                    */
   108  /* DECCOMBMSD and DECCOMBFROM are not format-dependent and so are     */
   109  /* only included once, when QUAD is being built                       */
   110  /* ------------------------------------------------------------------ */
   111  static const uInt DECCOMBEXP[64]={
   112    0, 0, 0, 0, 0, 0, 0, 0,
   113    1<<DECECONL, 1<<DECECONL, 1<<DECECONL, 1<<DECECONL,
   114    1<<DECECONL, 1<<DECECONL, 1<<DECECONL, 1<<DECECONL,
   115    2<<DECECONL, 2<<DECECONL, 2<<DECECONL, 2<<DECECONL,
   116    2<<DECECONL, 2<<DECECONL, 2<<DECECONL, 2<<DECECONL,
   117    0,           0,           1<<DECECONL, 1<<DECECONL,
   118    2<<DECECONL, 2<<DECECONL, DECFLOAT_Inf, DECFLOAT_NaN,
   119    0, 0, 0, 0, 0, 0, 0, 0,
   120    1<<DECECONL, 1<<DECECONL, 1<<DECECONL, 1<<DECECONL,
   121    1<<DECECONL, 1<<DECECONL, 1<<DECECONL, 1<<DECECONL,
   122    2<<DECECONL, 2<<DECECONL, 2<<DECECONL, 2<<DECECONL,
   123    2<<DECECONL, 2<<DECECONL, 2<<DECECONL, 2<<DECECONL,
   124    0,           0,           1<<DECECONL, 1<<DECECONL,
   125    2<<DECECONL, 2<<DECECONL, DECFLOAT_Inf, DECFLOAT_NaN};
   126  #if !QUAD
   127  static const uInt DECCOMBWEXP[64]={
   128    0, 0, 0, 0, 0, 0, 0, 0,
   129    1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL,
   130    1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL,
   131    2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL,
   132    2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL,
   133    0,            0,            1<<DECWECONL, 1<<DECWECONL,
   134    2<<DECWECONL, 2<<DECWECONL, DECFLOAT_Inf, DECFLOAT_NaN,
   135    0, 0, 0, 0, 0, 0, 0, 0,
   136    1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL,
   137    1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL, 1<<DECWECONL,
   138    2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL,
   139    2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL, 2<<DECWECONL,
   140    0,            0,            1<<DECWECONL, 1<<DECWECONL,
   141    2<<DECWECONL, 2<<DECWECONL, DECFLOAT_Inf, DECFLOAT_NaN};
   142  #endif
   143  
   144  #if QUAD
   145  const uInt DECCOMBMSD[64]={
   146    0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
   147    0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, 0, 0,
   148    0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
   149    0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, 0, 0};
   150  
   151  const uInt DECCOMBFROM[48]={
   152    0x00000000, 0x04000000, 0x08000000, 0x0C000000, 0x10000000, 0x14000000,
   153    0x18000000, 0x1C000000, 0x60000000, 0x64000000, 0x00000000, 0x00000000,
   154    0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x20000000, 0x24000000,
   155    0x28000000, 0x2C000000, 0x30000000, 0x34000000, 0x38000000, 0x3C000000,
   156    0x68000000, 0x6C000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000,
   157    0x00000000, 0x00000000, 0x40000000, 0x44000000, 0x48000000, 0x4C000000,
   158    0x50000000, 0x54000000, 0x58000000, 0x5C000000, 0x70000000, 0x74000000,
   159    0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000, 0x00000000};
   160  
   161  /* ------------------------------------------------------------------ */
   162  /* Request and include the tables to use for conversions              */
   163  /* ------------------------------------------------------------------ */
   164  #define DEC_BCD2DPD  1        // 0-0x999 -> DPD
   165  #define DEC_BIN2DPD  1        // 0-999 -> DPD
   166  #define DEC_BIN2BCD8 1        // 0-999 -> ddd, len
   167  #define DEC_DPD2BCD8 1        // DPD -> ddd, len
   168  #define DEC_DPD2BIN  1        // DPD -> 0-999
   169  #define DEC_DPD2BINK 1        // DPD -> 0-999000
   170  #define DEC_DPD2BINM 1        // DPD -> 0-999000000
   171  #include "decDPD.h"           // source of the lookup tables
   172  
   173  #endif
   174  
   175  /* ----------------------------------------------------------------- */
   176  /* decBiStr -- compare string with pairwise options                  */
   177  /*                                                                   */
   178  /*   targ is the string to compare                                   */
   179  /*   str1 is one of the strings to compare against (length may be 0) */
   180  /*   str2 is the other; it must be the same length as str1           */
   181  /*                                                                   */
   182  /*   returns 1 if strings compare equal, (that is, targ is the same  */
   183  /*   length as str1 and str2, and each character of targ is in one   */
   184  /*   of str1 or str2 in the corresponding position), or 0 otherwise  */
   185  /*                                                                   */
   186  /* This is used for generic caseless compare, including the awkward  */
   187  /* case of the Turkish dotted and dotless Is.  Use as (for example): */
   188  /*   if (decBiStr(test, "mike", "MIKE")) ...                         */
   189  /* ----------------------------------------------------------------- */
   190  static Flag decBiStr(const char *targ, const char *str1, const char *str2) {
   191    for (;;targ++, str1++, str2++) {
   192      if (*targ!=*str1 && *targ!=*str2) return 0;
   193      // *targ has a match in one (or both, if terminator)
   194      if (*targ=='\0') break;
   195      } // forever
   196    return 1;
   197    } // decBiStr
   198  
   199  /* ------------------------------------------------------------------ */
   200  /* decFinalize -- adjust and store a final result                     */
   201  /*                                                                    */
   202  /*  df  is the decFloat format number which gets the final result     */
   203  /*  num is the descriptor of the number to be checked and encoded     */
   204  /*         [its values, including the coefficient, may be modified]   */
   205  /*  set is the context to use                                         */
   206  /*  returns df                                                        */
   207  /*                                                                    */
   208  /* The num descriptor may point to a bcd8 string of any length; this  */
   209  /* string may have leading insignificant zeros.  If it has more than  */
   210  /* DECPMAX digits then the final digit can be a round-for-reround     */
   211  /* digit (i.e., it may include a sticky bit residue).                 */
   212  /*                                                                    */
   213  /* The exponent (q) may be one of the codes for a special value and   */
   214  /* can be up to 999999999 for conversion from string.                 */
   215  /*                                                                    */
   216  /* No error is possible, but Inexact, Underflow, and/or Overflow may  */
   217  /* be set.                                                            */
   218  /* ------------------------------------------------------------------ */
   219  // Constant whose size varies with format; also the check for surprises
   220  static uByte allnines[DECPMAX]=
   221  #if SINGLE
   222    {9, 9, 9, 9, 9, 9, 9};
   223  #elif DOUBLE
   224    {9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9};
   225  #elif QUAD
   226    {9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
   227     9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9};
   228  #endif
   229  
   230  static decFloat * decFinalize(decFloat *df, bcdnum *num,
   231                                decContext *set) {
   232    uByte *ub;                  // work
   233    uInt   dpd;                 // ..
   234    uInt   uiwork;              // for macros
   235    uByte *umsd=num->msd;       // local copy
   236    uByte *ulsd=num->lsd;       // ..
   237    uInt   encode;              // encoding accumulator
   238    Int    length;              // coefficient length
   239  
   240    #if DECCHECK
   241    Int clen=ulsd-umsd+1;
   242    #if QUAD
   243      #define COEXTRA 2                        // extra-long coefficent
   244    #else
   245      #define COEXTRA 0
   246    #endif
   247    if (clen<1 || clen>DECPMAX*3+2+COEXTRA)
   248      printf("decFinalize: suspect coefficient [length=%ld]\n", (LI)clen);
   249    if (num->sign!=0 && num->sign!=DECFLOAT_Sign)
   250      printf("decFinalize: bad sign [%08lx]\n", (LI)num->sign);
   251    if (!EXPISSPECIAL(num->exponent)
   252        && (num->exponent>1999999999 || num->exponent<-1999999999))
   253      printf("decFinalize: improbable exponent [%ld]\n", (LI)num->exponent);
   254    // decShowNum(num, "final");
   255    #endif
   256  
   257    // A special will have an 'exponent' which is very positive and a
   258    // coefficient < DECPMAX
   259    length=(uInt)(ulsd-umsd+1);                // coefficient length
   260  
   261    if (!NUMISSPECIAL(num)) {
   262      Int   drop;                              // digits to be dropped
   263      // skip leading insignificant zeros to calculate an exact length
   264      // [this is quite expensive]
   265      if (*umsd==0) {
   266        for (; umsd+3<ulsd && UBTOUI(umsd)==0;) umsd+=4;
   267        for (; *umsd==0 && umsd<ulsd;) umsd++;
   268        length=ulsd-umsd+1;                    // recalculate
   269        }
   270      drop=MAXI(length-DECPMAX, DECQTINY-num->exponent);
   271      // drop can now be > digits for bottom-clamp (subnormal) cases
   272      if (drop>0) {                            // rounding needed
   273        // (decFloatQuantize has very similar code to this, so any
   274        // changes may need to be made there, too)
   275        uByte *roundat;                        // -> re-round digit
   276        uByte reround;                         // reround value
   277        // printf("Rounding; drop=%ld\n", (LI)drop);
   278  
   279        num->exponent+=drop;                   // always update exponent
   280  
   281        // Three cases here:
   282        //   1. new LSD is in coefficient (almost always)
   283        //   2. new LSD is digit to left of coefficient (so MSD is
   284        //      round-for-reround digit)
   285        //   3. new LSD is to left of case 2 (whole coefficient is sticky)
   286        // [duplicate check-stickies code to save a test]
   287        // [by-digit check for stickies as runs of zeros are rare]
   288        if (drop<length) {                     // NB lengths not addresses
   289          roundat=umsd+length-drop;
   290          reround=*roundat;
   291          for (ub=roundat+1; ub<=ulsd; ub++) {
   292            if (*ub!=0) {                      // non-zero to be discarded
   293              reround=DECSTICKYTAB[reround];   // apply sticky bit
   294              break;                           // [remainder don't-care]
   295              }
   296            } // check stickies
   297          ulsd=roundat-1;                      // new LSD
   298          }
   299         else {                                // edge case
   300          if (drop==length) {
   301            roundat=umsd;
   302            reround=*roundat;
   303            }
   304           else {
   305            roundat=umsd-1;
   306            reround=0;
   307            }
   308          for (ub=roundat+1; ub<=ulsd; ub++) {
   309            if (*ub!=0) {                      // non-zero to be discarded
   310              reround=DECSTICKYTAB[reround];   // apply sticky bit
   311              break;                           // [remainder don't-care]
   312              }
   313            } // check stickies
   314          *umsd=0;                             // coefficient is a 0
   315          ulsd=umsd;                           // ..
   316          }
   317  
   318        if (reround!=0) {                      // discarding non-zero
   319          uInt bump=0;
   320          set->status|=DEC_Inexact;
   321          // if adjusted exponent [exp+digits-1] is < EMIN then num is
   322          // subnormal -- so raise Underflow
   323          if (num->exponent<DECEMIN && (num->exponent+(ulsd-umsd+1)-1)<DECEMIN)
   324            set->status|=DEC_Underflow;
   325  
   326          // next decide whether increment of the coefficient is needed
   327          if (set->round==DEC_ROUND_HALF_EVEN) {    // fastpath slowest case
   328            if (reround>5) bump=1;                  // >0.5 goes up
   329             else if (reround==5)                   // exactly 0.5000 ..
   330              bump=*ulsd & 0x01;                    // .. up iff [new] lsd is odd
   331            } // r-h-e
   332           else switch (set->round) {
   333            case DEC_ROUND_DOWN: {
   334              // no change
   335              break;} // r-d
   336            case DEC_ROUND_HALF_DOWN: {
   337              if (reround>5) bump=1;
   338              break;} // r-h-d
   339            case DEC_ROUND_HALF_UP: {
   340              if (reround>=5) bump=1;
   341              break;} // r-h-u
   342            case DEC_ROUND_UP: {
   343              if (reround>0) bump=1;
   344              break;} // r-u
   345            case DEC_ROUND_CEILING: {
   346              // same as _UP for positive numbers, and as _DOWN for negatives
   347              if (!num->sign && reround>0) bump=1;
   348              break;} // r-c
   349            case DEC_ROUND_FLOOR: {
   350              // same as _UP for negative numbers, and as _DOWN for positive
   351              // [negative reround cannot occur on 0]
   352              if (num->sign && reround>0) bump=1;
   353              break;} // r-f
   354            case DEC_ROUND_05UP: {
   355              if (reround>0) { // anything out there is 'sticky'
   356                // bump iff lsd=0 or 5; this cannot carry so it could be
   357                // effected immediately with no bump -- but the code
   358                // is clearer if this is done the same way as the others
   359                if (*ulsd==0 || *ulsd==5) bump=1;
   360                }
   361              break;} // r-r
   362            default: {      // e.g., DEC_ROUND_MAX
   363              set->status|=DEC_Invalid_context;
   364              #if DECCHECK
   365              printf("Unknown rounding mode: %ld\n", (LI)set->round);
   366              #endif
   367              break;}
   368            } // switch (not r-h-e)
   369          // printf("ReRound: %ld  bump: %ld\n", (LI)reround, (LI)bump);
   370  
   371          if (bump!=0) {                       // need increment
   372            // increment the coefficient; this might end up with 1000...
   373            // (after the all nines case)
   374            ub=ulsd;
   375            for(; ub-3>=umsd && UBTOUI(ub-3)==0x09090909; ub-=4)  {
   376              UBFROMUI(ub-3, 0);               // to 00000000
   377              }
   378            // [note ub could now be to left of msd, and it is not safe
   379            // to write to the the left of the msd]
   380            // now at most 3 digits left to non-9 (usually just the one)
   381            for (; ub>=umsd; *ub=0, ub--) {
   382              if (*ub==9) continue;            // carry
   383              *ub+=1;
   384              break;
   385              }
   386            if (ub<umsd) {                     // had all-nines
   387              *umsd=1;                         // coefficient to 1000...
   388              // usually the 1000... coefficient can be used as-is
   389              if ((ulsd-umsd+1)==DECPMAX) {
   390                num->exponent++;
   391                }
   392               else {
   393                // if coefficient is shorter than Pmax then num is
   394                // subnormal, so extend it; this is safe as drop>0
   395                // (or, if the coefficient was supplied above, it could
   396                // not be 9); this may make the result normal.
   397                ulsd++;
   398                *ulsd=0;
   399                // [exponent unchanged]
   400                #if DECCHECK
   401                if (num->exponent!=DECQTINY) // sanity check
   402                  printf("decFinalize: bad all-nines extend [^%ld, %ld]\n",
   403                         (LI)num->exponent, (LI)(ulsd-umsd+1));
   404                #endif
   405                } // subnormal extend
   406              } // had all-nines
   407            } // bump needed
   408          } // inexact rounding
   409  
   410        length=ulsd-umsd+1;               // recalculate (may be <DECPMAX)
   411        } // need round (drop>0)
   412  
   413      // The coefficient will now fit and has final length unless overflow
   414      // decShowNum(num, "rounded");
   415  
   416      // if exponent is >=emax may have to clamp, overflow, or fold-down
   417      if (num->exponent>DECEMAX-(DECPMAX-1)) { // is edge case
   418        // printf("overflow checks...\n");
   419        if (*ulsd==0 && ulsd==umsd) {     // have zero
   420          num->exponent=DECEMAX-(DECPMAX-1); // clamp to max
   421          }
   422         else if ((num->exponent+length-1)>DECEMAX) { // > Nmax
   423          // Overflow -- these could go straight to encoding, here, but
   424          // instead num is adjusted to keep the code cleaner
   425          Flag needmax=0;                 // 1 for finite result
   426          set->status|=(DEC_Overflow | DEC_Inexact);
   427          switch (set->round) {
   428            case DEC_ROUND_DOWN: {
   429              needmax=1;                  // never Infinity
   430              break;} // r-d
   431            case DEC_ROUND_05UP: {
   432              needmax=1;                  // never Infinity
   433              break;} // r-05
   434            case DEC_ROUND_CEILING: {
   435              if (num->sign) needmax=1;   // Infinity iff non-negative
   436              break;} // r-c
   437            case DEC_ROUND_FLOOR: {
   438              if (!num->sign) needmax=1;  // Infinity iff negative
   439              break;} // r-f
   440            default: break;               // Infinity in all other cases
   441            }
   442          if (!needmax) {                 // easy .. set Infinity
   443            num->exponent=DECFLOAT_Inf;
   444            *umsd=0;                      // be clean: coefficient to 0
   445            ulsd=umsd;                    // ..
   446            }
   447           else {                         // return Nmax
   448            umsd=allnines;                // use constant array
   449            ulsd=allnines+DECPMAX-1;
   450            num->exponent=DECEMAX-(DECPMAX-1);
   451            }
   452          }
   453         else { // no overflow but non-zero and may have to fold-down
   454          Int shift=num->exponent-(DECEMAX-(DECPMAX-1));
   455          if (shift>0) {                  // fold-down needed
   456            // fold down needed; must copy to buffer in order to pad
   457            // with zeros safely; fortunately this is not the worst case
   458            // path because cannot have had a round
   459            uByte buffer[ROUNDUP(DECPMAX+3, 4)]; // [+3 allows uInt padding]
   460            uByte *s=umsd;                // source
   461            uByte *t=buffer;              // safe target
   462            uByte *tlsd=buffer+(ulsd-umsd)+shift; // target LSD
   463            // printf("folddown shift=%ld\n", (LI)shift);
   464            for (; s<=ulsd; s+=4, t+=4) UBFROMUI(t, UBTOUI(s));
   465            for (t=tlsd-shift+1; t<=tlsd; t+=4) UBFROMUI(t, 0);  // pad 0s
   466            num->exponent-=shift;
   467            umsd=buffer;
   468            ulsd=tlsd;
   469            }
   470          } // fold-down?
   471        length=ulsd-umsd+1;               // recalculate length
   472        } // high-end edge case
   473      } // finite number
   474  
   475    /*------------------------------------------------------------------*/
   476    /* At this point the result will properly fit the decFloat          */
   477    /* encoding, and it can be encoded with no possibility of error     */
   478    /*------------------------------------------------------------------*/
   479    // Following code does not alter coefficient (could be allnines array)
   480  
   481    // fast path possible when DECPMAX digits
   482    if (length==DECPMAX) {
   483      return decFloatFromBCD(df, num->exponent, umsd, num->sign);
   484      } // full-length
   485  
   486    // slower path when not a full-length number; must care about length
   487    // [coefficient length here will be < DECPMAX]
   488    if (!NUMISSPECIAL(num)) {             // is still finite
   489      // encode the combination field and exponent continuation
   490      uInt uexp=(uInt)(num->exponent+DECBIAS); // biased exponent
   491      uInt code=(uexp>>DECECONL)<<4;      // top two bits of exp
   492      // [msd==0]
   493      // look up the combination field and make high word
   494      encode=DECCOMBFROM[code];           // indexed by (0-2)*16+msd
   495      encode|=(uexp<<(32-6-DECECONL)) & 0x03ffffff; // exponent continuation
   496      }
   497     else encode=num->exponent;           // special [already in word]
   498    encode|=num->sign;                    // add sign
   499  
   500    // private macro to extract a declet, n (where 0<=n<DECLETS and 0
   501    // refers to the declet from the least significant three digits)
   502    // and put the corresponding DPD code into dpd.  Access to umsd and
   503    // ulsd (pointers to the most and least significant digit of the
   504    // variable-length coefficient) is assumed, along with use of a
   505    // working pointer, uInt *ub.
   506    // As not full-length then chances are there are many leading zeros
   507    // [and there may be a partial triad]
   508    #define getDPDt(dpd, n) ub=ulsd-(3*(n))-2;                          \
   509      if (ub<umsd-2) dpd=0;                                             \
   510       else if (ub>=umsd) dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];  \
   511       else {dpd=*(ub+2); if (ub+1==umsd) dpd+=*(ub+1)*16; dpd=BCD2DPD[dpd];}
   512  
   513    // place the declets in the encoding words and copy to result (df),
   514    // according to endianness; in all cases complete the sign word
   515    // first
   516    #if DECPMAX==7
   517      getDPDt(dpd, 1);
   518      encode|=dpd<<10;
   519      getDPDt(dpd, 0);
   520      encode|=dpd;
   521      DFWORD(df, 0)=encode;     // just the one word
   522  
   523    #elif DECPMAX==16
   524      getDPDt(dpd, 4); encode|=dpd<<8;
   525      getDPDt(dpd, 3); encode|=dpd>>2;
   526      DFWORD(df, 0)=encode;
   527      encode=dpd<<30;
   528      getDPDt(dpd, 2); encode|=dpd<<20;
   529      getDPDt(dpd, 1); encode|=dpd<<10;
   530      getDPDt(dpd, 0); encode|=dpd;
   531      DFWORD(df, 1)=encode;
   532  
   533    #elif DECPMAX==34
   534      getDPDt(dpd,10); encode|=dpd<<4;
   535      getDPDt(dpd, 9); encode|=dpd>>6;
   536      DFWORD(df, 0)=encode;
   537  
   538      encode=dpd<<26;
   539      getDPDt(dpd, 8); encode|=dpd<<16;
   540      getDPDt(dpd, 7); encode|=dpd<<6;
   541      getDPDt(dpd, 6); encode|=dpd>>4;
   542      DFWORD(df, 1)=encode;
   543  
   544      encode=dpd<<28;
   545      getDPDt(dpd, 5); encode|=dpd<<18;
   546      getDPDt(dpd, 4); encode|=dpd<<8;
   547      getDPDt(dpd, 3); encode|=dpd>>2;
   548      DFWORD(df, 2)=encode;
   549  
   550      encode=dpd<<30;
   551      getDPDt(dpd, 2); encode|=dpd<<20;
   552      getDPDt(dpd, 1); encode|=dpd<<10;
   553      getDPDt(dpd, 0); encode|=dpd;
   554      DFWORD(df, 3)=encode;
   555    #endif
   556  
   557    // printf("Status: %08lx\n", (LI)set->status);
   558    // decFloatShow(df, "final2");
   559    return df;
   560    } // decFinalize
   561  
   562  /* ------------------------------------------------------------------ */
   563  /* decFloatFromBCD -- set decFloat from exponent, BCD8, and sign      */
   564  /*                                                                    */
   565  /*  df is the target decFloat                                         */
   566  /*  exp is the in-range unbiased exponent, q, or a special value in   */
   567  /*    the form returned by decFloatGetExponent                        */
   568  /*  bcdar holds DECPMAX digits to set the coefficient from, one       */
   569  /*    digit in each byte (BCD8 encoding); the first (MSD) is ignored  */
   570  /*    if df is a NaN; all are ignored if df is infinite.              */
   571  /*    All bytes must be in 0-9; results are undefined otherwise.      */
   572  /*  sig is DECFLOAT_Sign to set the sign bit, 0 otherwise             */
   573  /*  returns df, which will be canonical                               */
   574  /*                                                                    */
   575  /* No error is possible, and no status will be set.                   */
   576  /* ------------------------------------------------------------------ */
   577  decFloat * decFloatFromBCD(decFloat *df, Int exp, const uByte *bcdar,
   578                             Int sig) {
   579    uInt encode, dpd;                     // work
   580    const uByte *ub;                      // ..
   581  
   582    if (EXPISSPECIAL(exp)) encode=exp|sig;// specials already encoded
   583     else {                               // is finite
   584      // encode the combination field and exponent continuation
   585      uInt uexp=(uInt)(exp+DECBIAS);      // biased exponent
   586      uInt code=(uexp>>DECECONL)<<4;      // top two bits of exp
   587      code+=bcdar[0];                     // add msd
   588      // look up the combination field and make high word
   589      encode=DECCOMBFROM[code]|sig;       // indexed by (0-2)*16+msd
   590      encode|=(uexp<<(32-6-DECECONL)) & 0x03ffffff; // exponent continuation
   591      }
   592  
   593    // private macro to extract a declet, n (where 0<=n<DECLETS and 0
   594    // refers to the declet from the least significant three digits)
   595    // and put the corresponding DPD code into dpd.
   596    // Use of a working pointer, uInt *ub, is assumed.
   597  
   598    #define getDPDb(dpd, n) ub=bcdar+DECPMAX-1-(3*(n))-2;     \
   599      dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
   600  
   601    // place the declets in the encoding words and copy to result (df),
   602    // according to endianness; in all cases complete the sign word
   603    // first
   604    #if DECPMAX==7
   605      getDPDb(dpd, 1);
   606      encode|=dpd<<10;
   607      getDPDb(dpd, 0);
   608      encode|=dpd;
   609      DFWORD(df, 0)=encode;     // just the one word
   610  
   611    #elif DECPMAX==16
   612      getDPDb(dpd, 4); encode|=dpd<<8;
   613      getDPDb(dpd, 3); encode|=dpd>>2;
   614      DFWORD(df, 0)=encode;
   615      encode=dpd<<30;
   616      getDPDb(dpd, 2); encode|=dpd<<20;
   617      getDPDb(dpd, 1); encode|=dpd<<10;
   618      getDPDb(dpd, 0); encode|=dpd;
   619      DFWORD(df, 1)=encode;
   620  
   621    #elif DECPMAX==34
   622      getDPDb(dpd,10); encode|=dpd<<4;
   623      getDPDb(dpd, 9); encode|=dpd>>6;
   624      DFWORD(df, 0)=encode;
   625  
   626      encode=dpd<<26;
   627      getDPDb(dpd, 8); encode|=dpd<<16;
   628      getDPDb(dpd, 7); encode|=dpd<<6;
   629      getDPDb(dpd, 6); encode|=dpd>>4;
   630      DFWORD(df, 1)=encode;
   631  
   632      encode=dpd<<28;
   633      getDPDb(dpd, 5); encode|=dpd<<18;
   634      getDPDb(dpd, 4); encode|=dpd<<8;
   635      getDPDb(dpd, 3); encode|=dpd>>2;
   636      DFWORD(df, 2)=encode;
   637  
   638      encode=dpd<<30;
   639      getDPDb(dpd, 2); encode|=dpd<<20;
   640      getDPDb(dpd, 1); encode|=dpd<<10;
   641      getDPDb(dpd, 0); encode|=dpd;
   642      DFWORD(df, 3)=encode;
   643    #endif
   644    // decFloatShow(df, "fromB");
   645    return df;
   646    } // decFloatFromBCD
   647  
   648  /* ------------------------------------------------------------------ */
   649  /* decFloatFromPacked -- set decFloat from exponent and packed BCD    */
   650  /*                                                                    */
   651  /*  df is the target decFloat                                         */
   652  /*  exp is the in-range unbiased exponent, q, or a special value in   */
   653  /*    the form returned by decFloatGetExponent                        */
   654  /*  packed holds DECPMAX packed decimal digits plus a sign nibble     */
   655  /*    (all 6 codes are OK); the first (MSD) is ignored if df is a NaN */
   656  /*    and all except sign are ignored if df is infinite.  For DOUBLE  */
   657  /*    and QUAD the first (pad) nibble is also ignored in all cases.   */
   658  /*    All coefficient nibbles must be in 0-9 and sign in A-F; results */
   659  /*    are undefined otherwise.                                        */
   660  /*  returns df, which will be canonical                               */
   661  /*                                                                    */
   662  /* No error is possible, and no status will be set.                   */
   663  /* ------------------------------------------------------------------ */
   664  decFloat * decFloatFromPacked(decFloat *df, Int exp, const uByte *packed) {
   665    uByte bcdar[DECPMAX+2];               // work [+1 for pad, +1 for sign]
   666    const uByte *ip;                      // ..
   667    uByte *op;                            // ..
   668    Int   sig=0;                          // sign
   669  
   670    // expand coefficient and sign to BCDAR
   671    #if SINGLE
   672    op=bcdar+1;                           // no pad digit
   673    #else
   674    op=bcdar;                             // first (pad) digit ignored
   675    #endif
   676    for (ip=packed; ip<packed+((DECPMAX+2)/2); ip++) {
   677      *op++=*ip>>4;
   678      *op++=(uByte)(*ip&0x0f);            // [final nibble is sign]
   679      }
   680    op--;                                 // -> sign byte
   681    if (*op==DECPMINUS || *op==DECPMINUSALT) sig=DECFLOAT_Sign;
   682  
   683    if (EXPISSPECIAL(exp)) {              // Infinity or NaN
   684      if (!EXPISINF(exp)) bcdar[1]=0;     // a NaN: ignore MSD
   685       else memset(bcdar+1, 0, DECPMAX);  // Infinite: coefficient to 0
   686      }
   687    return decFloatFromBCD(df, exp, bcdar+1, sig);
   688    } // decFloatFromPacked
   689  
   690  /* ------------------------------------------------------------------ */
   691  /* decFloatFromPackedChecked -- set from exponent and packed; checked */
   692  /*                                                                    */
   693  /*  df is the target decFloat                                         */
   694  /*  exp is the in-range unbiased exponent, q, or a special value in   */
   695  /*    the form returned by decFloatGetExponent                        */
   696  /*  packed holds DECPMAX packed decimal digits plus a sign nibble     */
   697  /*    (all 6 codes are OK); the first (MSD) must be 0 if df is a NaN  */
   698  /*    and all digits must be 0 if df is infinite.  For DOUBLE and     */
   699  /*    QUAD the first (pad) nibble must be 0.                          */
   700  /*    All coefficient nibbles must be in 0-9 and sign in A-F.         */
   701  /*  returns df, which will be canonical or NULL if any of the         */
   702  /*    requirements are not met (if this case df is unchanged); that   */
   703  /*    is, the input data must be as returned by decFloatToPacked,     */
   704  /*    except that all six sign codes are acccepted.                   */
   705  /*                                                                    */
   706  /* No status will be set.                                             */
   707  /* ------------------------------------------------------------------ */
   708  decFloat * decFloatFromPackedChecked(decFloat *df, Int exp,
   709                                       const uByte *packed) {
   710    uByte bcdar[DECPMAX+2];               // work [+1 for pad, +1 for sign]
   711    const uByte *ip;                      // ..
   712    uByte *op;                            // ..
   713    Int   sig=0;                          // sign
   714  
   715    // expand coefficient and sign to BCDAR
   716    #if SINGLE
   717    op=bcdar+1;                           // no pad digit
   718    #else
   719    op=bcdar;                             // first (pad) digit here
   720    #endif
   721    for (ip=packed; ip<packed+((DECPMAX+2)/2); ip++) {
   722      *op=*ip>>4;
   723      if (*op>9) return NULL;
   724      op++;
   725      *op=(uByte)(*ip&0x0f);              // [final nibble is sign]
   726      if (*op>9 && ip<packed+((DECPMAX+2)/2)-1) return NULL;
   727      op++;
   728      }
   729    op--;                                 // -> sign byte
   730    if (*op<=9) return NULL;              // bad sign
   731    if (*op==DECPMINUS || *op==DECPMINUSALT) sig=DECFLOAT_Sign;
   732  
   733    #if !SINGLE
   734    if (bcdar[0]!=0) return NULL;         // bad pad nibble
   735    #endif
   736  
   737    if (EXPISNAN(exp)) {                  // a NaN
   738      if (bcdar[1]!=0) return NULL;       // bad msd
   739      } // NaN
   740     else if (EXPISINF(exp)) {            // is infinite
   741      Int i;
   742      for (i=0; i<DECPMAX; i++) {
   743        if (bcdar[i+1]!=0) return NULL;   // should be all zeros
   744        }
   745      } // infinity
   746     else {                               // finite
   747      // check the exponent is in range
   748      if (exp>DECEMAX-DECPMAX+1) return NULL;
   749      if (exp<DECEMIN-DECPMAX+1) return NULL;
   750      }
   751    return decFloatFromBCD(df, exp, bcdar+1, sig);
   752    } // decFloatFromPacked
   753  
   754  /* ------------------------------------------------------------------ */
   755  /* decFloatFromString -- conversion from numeric string               */
   756  /*                                                                    */
   757  /*  result  is the decFloat format number which gets the result of    */
   758  /*          the conversion                                            */
   759  /*  *string is the character string which should contain a valid      */
   760  /*          number (which may be a special value), \0-terminated      */
   761  /*          If there are too many significant digits in the           */
   762  /*          coefficient it will be rounded.                           */
   763  /*  set     is the context                                            */
   764  /*  returns result                                                    */
   765  /*                                                                    */
   766  /* The length of the coefficient and the size of the exponent are     */
   767  /* checked by this routine, so the correct error (Underflow or        */
   768  /* Overflow) can be reported or rounding applied, as necessary.       */
   769  /*                                                                    */
   770  /* There is no limit to the coefficient length for finite inputs;     */
   771  /* NaN payloads must be integers with no more than DECPMAX-1 digits.  */
   772  /* Exponents may have up to nine significant digits.                  */
   773  /*                                                                    */
   774  /* If bad syntax is detected, the result will be a quiet NaN.         */
   775  /* ------------------------------------------------------------------ */
   776  decFloat * decFloatFromString(decFloat *result, const char *string,
   777                                decContext *set) {
   778    Int    digits;                   // count of digits in coefficient
   779    const  char *dotchar=NULL;       // where dot was found [NULL if none]
   780    const  char *cfirst=string;      // -> first character of decimal part
   781    const  char *c;                  // work
   782    uByte *ub;                       // ..
   783    uInt   uiwork;                   // for macros
   784    bcdnum num;                      // collects data for finishing
   785    uInt   error=DEC_Conversion_syntax;      // assume the worst
   786    uByte  buffer[ROUNDUP(DECSTRING+11, 8)]; // room for most coefficents,
   787                                             // some common rounding, +3, & pad
   788    #if DECTRACE
   789    // printf("FromString %s ...\n", string);
   790    #endif
   791  
   792    for(;;) {                             // once-only 'loop'
   793      num.sign=0;                         // assume non-negative
   794      num.msd=buffer;                     // MSD is here always
   795  
   796      // detect and validate the coefficient, including any leading,
   797      // trailing, or embedded '.'
   798      // [could test four-at-a-time here (saving 10% for decQuads),
   799      // but that risks storage violation because the position of the
   800      // terminator is unknown]
   801      for (c=string;; c++) {              // -> input character
   802        if (((unsigned)(*c-'0'))<=9) continue; // '0' through '9' is good
   803        if (*c=='\0') break;              // most common non-digit
   804        if (*c=='.') {
   805          if (dotchar!=NULL) break;       // not first '.'
   806          dotchar=c;                      // record offset into decimal part
   807          continue;}
   808        if (c==string) {                  // first in string...
   809          if (*c=='-') {                  // valid - sign
   810            cfirst++;
   811            num.sign=DECFLOAT_Sign;
   812            continue;}
   813          if (*c=='+') {                  // valid + sign
   814            cfirst++;
   815            continue;}
   816          }
   817        // *c is not a digit, terminator, or a valid +, -, or '.'
   818        break;
   819        } // c loop
   820  
   821      digits=(uInt)(c-cfirst);            // digits (+1 if a dot)
   822  
   823      if (digits>0) {                     // had digits and/or dot
   824        const char *clast=c-1;            // note last coefficient char position
   825        Int exp=0;                        // exponent accumulator
   826        if (*c!='\0') {                   // something follows the coefficient
   827          uInt edig;                      // unsigned work
   828          // had some digits and more to come; expect E[+|-]nnn now
   829          const char *firstexp;           // exponent first non-zero
   830          if (*c!='E' && *c!='e') break;
   831          c++;                            // to (optional) sign
   832          if (*c=='-' || *c=='+') c++;    // step over sign (c=clast+2)
   833          if (*c=='\0') break;            // no digits!  (e.g., '1.2E')
   834          for (; *c=='0';) c++;           // skip leading zeros [even last]
   835          firstexp=c;                     // remember start [maybe '\0']
   836          // gather exponent digits
   837          edig=(uInt)*c-(uInt)'0';
   838          if (edig<=9) {                  // [check not bad or terminator]
   839            exp+=edig;                    // avoid initial X10
   840            c++;
   841            for (;; c++) {
   842              edig=(uInt)*c-(uInt)'0';
   843              if (edig>9) break;
   844              exp=exp*10+edig;
   845              }
   846            }
   847          // if not now on the '\0', *c must not be a digit
   848          if (*c!='\0') break;
   849  
   850          // (this next test must be after the syntax checks)
   851          // if definitely more than the possible digits for format then
   852          // the exponent may have wrapped, so simply set it to a certain
   853          // over/underflow value
   854          if (c>firstexp+DECEMAXD) exp=DECEMAX*2;
   855          if (*(clast+2)=='-') exp=-exp;  // was negative
   856          } // exponent part
   857  
   858        if (dotchar!=NULL) {              // had a '.'
   859          digits--;                       // remove from digits count
   860          if (digits==0) break;           // was dot alone: bad syntax
   861          exp-=(Int)(clast-dotchar);      // adjust exponent
   862          // [the '.' can now be ignored]
   863          }
   864        num.exponent=exp;                 // exponent is good; store it
   865  
   866        // Here when whole string has been inspected and syntax is good
   867        // cfirst->first digit or dot, clast->last digit or dot
   868        error=0;                          // no error possible now
   869  
   870        // if the number of digits in the coefficient will fit in buffer
   871        // then it can simply be converted to bcd8 and copied -- decFinalize
   872        // will take care of leading zeros and rounding; the buffer is big
   873        // enough for all canonical coefficients, including 0.00000nn...
   874        ub=buffer;
   875        if (digits<=(Int)(sizeof(buffer)-3)) { // [-3 allows by-4s copy]
   876          c=cfirst;
   877          if (dotchar!=NULL) {                 // a dot to worry about
   878            if (*(c+1)=='.') {                 // common canonical case
   879              *ub++=(uByte)(*c-'0');           // copy leading digit
   880              c+=2;                            // prepare to handle rest
   881              }
   882             else for (; c<=clast;) {          // '.' could be anywhere
   883              // as usual, go by fours when safe; NB it has been asserted
   884              // that a '.' does not have the same mask as a digit
   885              if (c<=clast-3                             // safe for four
   886               && (UBTOUI(c)&0xf0f0f0f0)==CHARMASK) {    // test four
   887                UBFROMUI(ub, UBTOUI(c)&0x0f0f0f0f);      // to BCD8
   888                ub+=4;
   889                c+=4;
   890                continue;
   891                }
   892              if (*c=='.') {                   // found the dot
   893                c++;                           // step over it ..
   894                break;                         // .. and handle the rest
   895                }
   896              *ub++=(uByte)(*c++-'0');
   897              }
   898            } // had dot
   899          // Now no dot; do this by fours (where safe)
   900          for (; c<=clast-3; c+=4, ub+=4) UBFROMUI(ub, UBTOUI(c)&0x0f0f0f0f);
   901          for (; c<=clast; c++, ub++) *ub=(uByte)(*c-'0');
   902          num.lsd=buffer+digits-1;             // record new LSD
   903          } // fits
   904  
   905         else {                                // too long for buffer
   906          // [This is a rare and unusual case; arbitrary-length input]
   907          // strip leading zeros [but leave final 0 if all 0's]
   908          if (*cfirst=='.') cfirst++;          // step past dot at start
   909          if (*cfirst=='0') {                  // [cfirst always -> digit]
   910            for (; cfirst<clast; cfirst++) {
   911              if (*cfirst!='0') {              // non-zero found
   912                if (*cfirst=='.') continue;    // [ignore]
   913                break;                         // done
   914                }
   915              digits--;                        // 0 stripped
   916              } // cfirst
   917            } // at least one leading 0
   918  
   919          // the coefficient is now as short as possible, but may still
   920          // be too long; copy up to Pmax+1 digits to the buffer, then
   921          // just record any non-zeros (set round-for-reround digit)
   922          for (c=cfirst; c<=clast && ub<=buffer+DECPMAX; c++) {
   923            // (see commentary just above)
   924            if (c<=clast-3                          // safe for four
   925             && (UBTOUI(c)&0xf0f0f0f0)==CHARMASK) { // four digits
   926              UBFROMUI(ub, UBTOUI(c)&0x0f0f0f0f);   // to BCD8
   927              ub+=4;
   928              c+=3;                            // [will become 4]
   929              continue;
   930              }
   931            if (*c=='.') continue;             // [ignore]
   932            *ub++=(uByte)(*c-'0');
   933            }
   934          ub--;                                // -> LSD
   935          for (; c<=clast; c++) {              // inspect remaining chars
   936            if (*c!='0') {                     // sticky bit needed
   937              if (*c=='.') continue;           // [ignore]
   938              *ub=DECSTICKYTAB[*ub];           // update round-for-reround
   939              break;                           // no need to look at more
   940              }
   941            }
   942          num.lsd=ub;                          // record LSD
   943          // adjust exponent for dropped digits
   944          num.exponent+=digits-(Int)(ub-buffer+1);
   945          } // too long for buffer
   946        } // digits and/or dot
   947  
   948       else {                             // no digits or dot were found
   949        // only Infinities and NaNs are allowed, here
   950        if (*c=='\0') break;              // nothing there is bad
   951        buffer[0]=0;                      // default a coefficient of 0
   952        num.lsd=buffer;                   // ..
   953        if (decBiStr(c, "infinity", "INFINITY")
   954         || decBiStr(c, "inf", "INF")) num.exponent=DECFLOAT_Inf;
   955         else {                           // should be a NaN
   956          num.exponent=DECFLOAT_qNaN;     // assume quiet NaN
   957          if (*c=='s' || *c=='S') {       // probably an sNaN
   958            num.exponent=DECFLOAT_sNaN;   // effect the 's'
   959            c++;                          // and step over it
   960            }
   961          if (*c!='N' && *c!='n') break;  // check caseless "NaN"
   962          c++;
   963          if (*c!='a' && *c!='A') break;  // ..
   964          c++;
   965          if (*c!='N' && *c!='n') break;  // ..
   966          c++;
   967          // now either nothing, or nnnn payload (no dots), expected
   968          // -> start of integer, and skip leading 0s [including plain 0]
   969          for (cfirst=c; *cfirst=='0';) cfirst++;
   970          if (*cfirst!='\0') {            // not empty or all-0, payload
   971            // payload found; check all valid digits and copy to buffer as bcd8
   972            ub=buffer;
   973            for (c=cfirst;; c++, ub++) {
   974              if ((unsigned)(*c-'0')>9) break; // quit if not 0-9
   975              if (c-cfirst==DECPMAX-1) break;  // too many digits
   976              *ub=(uByte)(*c-'0');        // good bcd8
   977              }
   978            if (*c!='\0') break;          // not all digits, or too many
   979            num.lsd=ub-1;                 // record new LSD
   980            }
   981          } // NaN or sNaN
   982        error=0;                          // syntax is OK
   983        } // digits=0 (special expected)
   984      break;                              // drop out
   985      }                                   // [for(;;) once-loop]
   986  
   987    // decShowNum(&num, "fromStr");
   988  
   989    if (error!=0) {
   990      set->status|=error;
   991      num.exponent=DECFLOAT_qNaN;         // set up quiet NaN
   992      num.sign=0;                         // .. with 0 sign
   993      buffer[0]=0;                        // .. and coefficient
   994      num.lsd=buffer;                     // ..
   995      // decShowNum(&num, "oops");
   996      }
   997  
   998    // decShowNum(&num, "dffs");
   999    decFinalize(result, &num, set);       // round, check, and lay out
  1000    // decFloatShow(result, "fromString");
  1001    return result;
  1002    } // decFloatFromString
  1003  
  1004  /* ------------------------------------------------------------------ */
  1005  /* decFloatFromWider -- conversion from next-wider format             */
  1006  /*                                                                    */
  1007  /*  result  is the decFloat format number which gets the result of    */
  1008  /*          the conversion                                            */
  1009  /*  wider   is the decFloatWider format number which will be narrowed */
  1010  /*  set     is the context                                            */
  1011  /*  returns result                                                    */
  1012  /*                                                                    */
  1013  /* Narrowing can cause rounding, overflow, etc., but not Invalid      */
  1014  /* operation (sNaNs are copied and do not signal).                    */
  1015  /* ------------------------------------------------------------------ */
  1016  // narrow-to is not possible for decQuad format numbers; simply omit
  1017  #if !QUAD
  1018  decFloat * decFloatFromWider(decFloat *result, const decFloatWider *wider,
  1019                               decContext *set) {
  1020    bcdnum num;                           // collects data for finishing
  1021    uByte  bcdar[DECWPMAX];               // room for wider coefficient
  1022    uInt   widerhi=DFWWORD(wider, 0);     // top word
  1023    Int    exp;
  1024  
  1025    GETWCOEFF(wider, bcdar);
  1026  
  1027    num.msd=bcdar;                        // MSD is here always
  1028    num.lsd=bcdar+DECWPMAX-1;             // LSD is here always
  1029    num.sign=widerhi&0x80000000;          // extract sign [DECFLOAT_Sign=Neg]
  1030  
  1031    // decode the wider combination field to exponent
  1032    exp=DECCOMBWEXP[widerhi>>26];         // decode from wider combination field
  1033    // if it is a special there's nothing to do unless sNaN; if it's
  1034    // finite then add the (wider) exponent continuation and unbias
  1035    if (EXPISSPECIAL(exp)) exp=widerhi&0x7e000000; // include sNaN selector
  1036     else exp+=GETWECON(wider)-DECWBIAS;
  1037    num.exponent=exp;
  1038  
  1039    // decShowNum(&num, "dffw");
  1040    return decFinalize(result, &num, set);// round, check, and lay out
  1041    } // decFloatFromWider
  1042  #endif
  1043  
  1044  /* ------------------------------------------------------------------ */
  1045  /* decFloatGetCoefficient -- get coefficient as BCD8                  */
  1046  /*                                                                    */
  1047  /*  df is the decFloat from which to extract the coefficient          */
  1048  /*  bcdar is where DECPMAX bytes will be written, one BCD digit in    */
  1049  /*    each byte (BCD8 encoding); if df is a NaN the first byte will   */
  1050  /*    be zero, and if it is infinite they will all be zero            */
  1051  /*  returns the sign of the coefficient (DECFLOAT_Sign if negative,   */
  1052  /*    0 otherwise)                                                    */
  1053  /*                                                                    */
  1054  /* No error is possible, and no status will be set.  If df is a       */
  1055  /* special value the array is set to zeros (for Infinity) or to the   */
  1056  /* payload of a qNaN or sNaN.                                         */
  1057  /* ------------------------------------------------------------------ */
  1058  Int decFloatGetCoefficient(const decFloat *df, uByte *bcdar) {
  1059    if (DFISINF(df)) memset(bcdar, 0, DECPMAX);
  1060     else {
  1061      GETCOEFF(df, bcdar);           // use macro
  1062      if (DFISNAN(df)) bcdar[0]=0;   // MSD needs correcting
  1063      }
  1064    return GETSIGN(df);
  1065    } // decFloatGetCoefficient
  1066  
  1067  /* ------------------------------------------------------------------ */
  1068  /* decFloatGetExponent -- get unbiased exponent                       */
  1069  /*                                                                    */
  1070  /*  df is the decFloat from which to extract the exponent             */
  1071  /*  returns the exponent, q.                                          */
  1072  /*                                                                    */
  1073  /* No error is possible, and no status will be set.  If df is a       */
  1074  /* special value the first seven bits of the decFloat are returned,   */
  1075  /* left adjusted and with the first (sign) bit set to 0 (followed by  */
  1076  /* 25 0 bits).  e.g., -sNaN would return 0x7e000000 (DECFLOAT_sNaN).  */
  1077  /* ------------------------------------------------------------------ */
  1078  Int decFloatGetExponent(const decFloat *df) {
  1079    if (DFISSPECIAL(df)) return DFWORD(df, 0)&0x7e000000;
  1080    return GETEXPUN(df);
  1081    } // decFloatGetExponent
  1082  
  1083  /* ------------------------------------------------------------------ */
  1084  /* decFloatSetCoefficient -- set coefficient from BCD8                */
  1085  /*                                                                    */
  1086  /*  df is the target decFloat (and source of exponent/special value)  */
  1087  /*  bcdar holds DECPMAX digits to set the coefficient from, one       */
  1088  /*    digit in each byte (BCD8 encoding); the first (MSD) is ignored  */
  1089  /*    if df is a NaN; all are ignored if df is infinite.              */
  1090  /*  sig is DECFLOAT_Sign to set the sign bit, 0 otherwise             */
  1091  /*  returns df, which will be canonical                               */
  1092  /*                                                                    */
  1093  /* No error is possible, and no status will be set.                   */
  1094  /* ------------------------------------------------------------------ */
  1095  decFloat * decFloatSetCoefficient(decFloat *df, const uByte *bcdar,
  1096                                    Int sig) {
  1097    uInt exp;                        // for exponent
  1098    uByte bcdzero[DECPMAX];          // for infinities
  1099  
  1100    // Exponent/special code is extracted from df
  1101    if (DFISSPECIAL(df)) {
  1102      exp=DFWORD(df, 0)&0x7e000000;
  1103      if (DFISINF(df)) {
  1104        memset(bcdzero, 0, DECPMAX);
  1105        return decFloatFromBCD(df, exp, bcdzero, sig);
  1106        }
  1107      }
  1108     else exp=GETEXPUN(df);
  1109    return decFloatFromBCD(df, exp, bcdar, sig);
  1110    } // decFloatSetCoefficient
  1111  
  1112  /* ------------------------------------------------------------------ */
  1113  /* decFloatSetExponent -- set exponent or special value               */
  1114  /*                                                                    */
  1115  /*  df  is the target decFloat (and source of coefficient/payload)    */
  1116  /*  set is the context for reporting status                           */
  1117  /*  exp is the unbiased exponent, q, or a special value in the form   */
  1118  /*    returned by decFloatGetExponent                                 */
  1119  /*  returns df, which will be canonical                               */
  1120  /*                                                                    */
  1121  /* No error is possible, but Overflow or Underflow might occur.       */
  1122  /* ------------------------------------------------------------------ */
  1123  decFloat * decFloatSetExponent(decFloat *df, decContext *set, Int exp) {
  1124    uByte  bcdcopy[DECPMAX];         // for coefficient
  1125    bcdnum num;                      // work
  1126    num.exponent=exp;
  1127    num.sign=decFloatGetCoefficient(df, bcdcopy); // extract coefficient
  1128    if (DFISSPECIAL(df)) {           // MSD or more needs correcting
  1129      if (DFISINF(df)) memset(bcdcopy, 0, DECPMAX);
  1130      bcdcopy[0]=0;
  1131      }
  1132    num.msd=bcdcopy;
  1133    num.lsd=bcdcopy+DECPMAX-1;
  1134    return decFinalize(df, &num, set);
  1135    } // decFloatSetExponent
  1136  
  1137  /* ------------------------------------------------------------------ */
  1138  /* decFloatRadix -- returns the base (10)                             */
  1139  /*                                                                    */
  1140  /*   df is any decFloat of this format                                */
  1141  /* ------------------------------------------------------------------ */
  1142  uInt decFloatRadix(const decFloat *df) {
  1143    if (df) return 10;                         // to placate compiler
  1144    return 10;
  1145    } // decFloatRadix
  1146  
  1147  /* The following function is not available if DECPRINT=0              */
  1148  #if DECPRINT
  1149  /* ------------------------------------------------------------------ */
  1150  /* decFloatShow -- printf a decFloat in hexadecimal and decimal       */
  1151  /*   df  is the decFloat to show                                      */
  1152  /*   tag is a tag string displayed with the number                    */
  1153  /*                                                                    */
  1154  /* This is a debug aid; the precise format of the string may change.  */
  1155  /* ------------------------------------------------------------------ */
  1156  void decFloatShow(const decFloat *df, const char *tag) {
  1157    char hexbuf[DECBYTES*2+DECBYTES/4+1]; // NB blank after every fourth
  1158    char buff[DECSTRING];                 // for value in decimal
  1159    Int i, j=0;
  1160  
  1161    for (i=0; i<DECBYTES; i++) {
  1162      #if DECLITEND
  1163        sprintf(&hexbuf[j], "%02x", df->bytes[DECBYTES-1-i]);
  1164      #else
  1165        sprintf(&hexbuf[j], "%02x", df->bytes[i]);
  1166      #endif
  1167      j+=2;
  1168      // the next line adds blank (and terminator) after final pair, too
  1169      if ((i+1)%4==0) {strcpy(&hexbuf[j], " "); j++;}
  1170      }
  1171    decFloatToString(df, buff);
  1172    printf(">%s> %s [big-endian]  %s\n", tag, hexbuf, buff);
  1173    return;
  1174    } // decFloatShow
  1175  #endif
  1176  
  1177  /* ------------------------------------------------------------------ */
  1178  /* decFloatToBCD -- get sign, exponent, and BCD8 from a decFloat      */
  1179  /*                                                                    */
  1180  /*  df is the source decFloat                                         */
  1181  /*  exp will be set to the unbiased exponent, q, or to a special      */
  1182  /*    value in the form returned by decFloatGetExponent               */
  1183  /*  bcdar is where DECPMAX bytes will be written, one BCD digit in    */
  1184  /*    each byte (BCD8 encoding); if df is a NaN the first byte will   */
  1185  /*    be zero, and if it is infinite they will all be zero            */
  1186  /*  returns the sign of the coefficient (DECFLOAT_Sign if negative,   */
  1187  /*    0 otherwise)                                                    */
  1188  /*                                                                    */
  1189  /* No error is possible, and no status will be set.                   */
  1190  /* ------------------------------------------------------------------ */
  1191  Int decFloatToBCD(const decFloat *df, Int *exp, uByte *bcdar) {
  1192    if (DFISINF(df)) {
  1193      memset(bcdar, 0, DECPMAX);
  1194      *exp=DFWORD(df, 0)&0x7e000000;
  1195      }
  1196     else {
  1197      GETCOEFF(df, bcdar);           // use macro
  1198      if (DFISNAN(df)) {
  1199        bcdar[0]=0;                  // MSD needs correcting
  1200        *exp=DFWORD(df, 0)&0x7e000000;
  1201        }
  1202       else {                        // finite
  1203        *exp=GETEXPUN(df);
  1204        }
  1205      }
  1206    return GETSIGN(df);
  1207    } // decFloatToBCD
  1208  
  1209  /* ------------------------------------------------------------------ */
  1210  /* decFloatToEngString -- conversion to numeric string, engineering   */
  1211  /*                                                                    */
  1212  /*  df is the decFloat format number to convert                       */
  1213  /*  string is the string where the result will be laid out            */
  1214  /*                                                                    */
  1215  /* string must be at least DECPMAX+9 characters (the worst case is    */
  1216  /* "-0.00000nnn...nnn\0", which is as long as the exponent form when  */
  1217  /* DECEMAXD<=4); this condition is asserted above                     */
  1218  /*                                                                    */
  1219  /* No error is possible, and no status will be set                    */
  1220  /* ------------------------------------------------------------------ */
  1221  char * decFloatToEngString(const decFloat *df, char *string){
  1222    uInt msd;                        // coefficient MSD
  1223    Int  exp;                        // exponent top two bits or full
  1224    uInt comb;                       // combination field
  1225    char *cstart;                    // coefficient start
  1226    char *c;                         // output pointer in string
  1227    char *s, *t;                     // .. (source, target)
  1228    Int  pre, e;                     // work
  1229    const uByte *u;                  // ..
  1230    uInt  uiwork;                    // for macros [one compiler needs
  1231                                     // volatile here to avoid bug, but
  1232                                     // that doubles execution time]
  1233  
  1234    // Source words; macro handles endianness
  1235    uInt sourhi=DFWORD(df, 0);       // word with sign
  1236    #if DECPMAX==16
  1237    uInt sourlo=DFWORD(df, 1);
  1238    #elif DECPMAX==34
  1239    uInt sourmh=DFWORD(df, 1);
  1240    uInt sourml=DFWORD(df, 2);
  1241    uInt sourlo=DFWORD(df, 3);
  1242    #endif
  1243  
  1244    c=string;                        // where result will go
  1245    if (((Int)sourhi)<0) *c++='-';   // handle sign
  1246    comb=sourhi>>26;                 // sign+combination field
  1247    msd=DECCOMBMSD[comb];            // decode the combination field
  1248    exp=DECCOMBEXP[comb];            // ..
  1249  
  1250    if (EXPISSPECIAL(exp)) {         // special
  1251      if (exp==DECFLOAT_Inf) {       // infinity
  1252        strcpy(c,   "Inf");
  1253        strcpy(c+3, "inity");
  1254        return string;               // easy
  1255        }
  1256      if (sourhi&0x02000000) *c++='s'; // sNaN
  1257      strcpy(c, "NaN");              // complete word
  1258      c+=3;                          // step past
  1259      // quick exit if the payload is zero
  1260      #if DECPMAX==7
  1261      if ((sourhi&0x000fffff)==0) return string;
  1262      #elif DECPMAX==16
  1263      if (sourlo==0 && (sourhi&0x0003ffff)==0) return string;
  1264      #elif DECPMAX==34
  1265      if (sourlo==0 && sourml==0 && sourmh==0
  1266       && (sourhi&0x00003fff)==0) return string;
  1267      #endif
  1268      // otherwise drop through to add integer; set correct exp etc.
  1269      exp=0; msd=0;                  // setup for following code
  1270      }
  1271     else { // complete exponent; top two bits are in place
  1272      exp+=GETECON(df)-DECBIAS;      // .. + continuation and unbias
  1273      }
  1274  
  1275    /* convert the digits of the significand to characters */
  1276    cstart=c;                        // save start of coefficient
  1277    if (msd) *c++=(char)('0'+(char)msd);  // non-zero most significant digit
  1278  
  1279    // Decode the declets.  After extracting each declet, it is
  1280    // decoded to a 4-uByte sequence by table lookup; the four uBytes
  1281    // are the three encoded BCD8 digits followed by a 1-byte length
  1282    // (significant digits, except that 000 has length 0).  This allows
  1283    // us to left-align the first declet with non-zero content, then
  1284    // the remaining ones are full 3-char length.  Fixed-length copies
  1285    // are used because variable-length memcpy causes a subroutine call
  1286    // in at least two compilers.  (The copies are length 4 for speed
  1287    // and are safe because the last item in the array is of length
  1288    // three and has the length byte following.)
  1289    #define dpd2char(dpdin) u=&DPD2BCD8[((dpdin)&0x3ff)*4];        \
  1290           if (c!=cstart) {UBFROMUI(c, UBTOUI(u)|CHARMASK); c+=3;} \
  1291            else if (*(u+3)) {                                     \
  1292             UBFROMUI(c, UBTOUI(u+3-*(u+3))|CHARMASK); c+=*(u+3);}
  1293  
  1294    #if DECPMAX==7
  1295    dpd2char(sourhi>>10);                 // declet 1
  1296    dpd2char(sourhi);                     // declet 2
  1297  
  1298    #elif DECPMAX==16
  1299    dpd2char(sourhi>>8);                  // declet 1
  1300    dpd2char((sourhi<<2) | (sourlo>>30)); // declet 2
  1301    dpd2char(sourlo>>20);                 // declet 3
  1302    dpd2char(sourlo>>10);                 // declet 4
  1303    dpd2char(sourlo);                     // declet 5
  1304  
  1305    #elif DECPMAX==34
  1306    dpd2char(sourhi>>4);                  // declet 1
  1307    dpd2char((sourhi<<6) | (sourmh>>26)); // declet 2
  1308    dpd2char(sourmh>>16);                 // declet 3
  1309    dpd2char(sourmh>>6);                  // declet 4
  1310    dpd2char((sourmh<<4) | (sourml>>28)); // declet 5
  1311    dpd2char(sourml>>18);                 // declet 6
  1312    dpd2char(sourml>>8);                  // declet 7
  1313    dpd2char((sourml<<2) | (sourlo>>30)); // declet 8
  1314    dpd2char(sourlo>>20);                 // declet 9
  1315    dpd2char(sourlo>>10);                 // declet 10
  1316    dpd2char(sourlo);                     // declet 11
  1317    #endif
  1318  
  1319    if (c==cstart) *c++='0';         // all zeros, empty -- make "0"
  1320  
  1321    if (exp==0) {                    // integer or NaN case -- easy
  1322      *c='\0';                       // terminate
  1323      return string;
  1324      }
  1325    /* non-0 exponent */
  1326  
  1327    e=0;                             // assume no E
  1328    pre=(Int)(c-cstart)+exp;         // length+exp  [c->LSD+1]
  1329    // [here, pre-exp is the digits count (==1 for zero)]
  1330  
  1331    if (exp>0 || pre<-5) {           // need exponential form
  1332      e=pre-1;                       // calculate E value
  1333      pre=1;                         // assume one digit before '.'
  1334      if (e!=0) {                    // engineering: may need to adjust
  1335        Int adj;                     // adjustment
  1336        // The C remainder operator is undefined for negative numbers, so
  1337        // a positive remainder calculation must be used here
  1338        if (e<0) {
  1339          adj=(-e)%3;
  1340          if (adj!=0) adj=3-adj;
  1341          }
  1342         else { // e>0
  1343          adj=e%3;
  1344          }
  1345        e=e-adj;
  1346        // if dealing with zero still produce an exponent which is a
  1347        // multiple of three, as expected, but there will only be the
  1348        // one zero before the E, still.  Otherwise note the padding.
  1349        if (!DFISZERO(df)) pre+=adj;
  1350         else {  // is zero
  1351          if (adj!=0) {              // 0.00Esnn needed
  1352            e=e+3;
  1353            pre=-(2-adj);
  1354            }
  1355          } // zero
  1356        } // engineering adjustment
  1357      } // exponential form
  1358    // printf("e=%ld pre=%ld exp=%ld\n", (LI)e, (LI)pre, (LI)exp);
  1359  
  1360    /* modify the coefficient, adding 0s, '.', and E+nn as needed */
  1361    if (pre>0) {                     // ddd.ddd (plain), perhaps with E
  1362                                     // or dd00 padding for engineering
  1363      char *dotat=cstart+pre;
  1364      if (dotat<c) {                      // if embedded dot needed...
  1365        // move by fours; there must be space for junk at the end
  1366        // because there is still space for exponent
  1367        s=dotat+ROUNDDOWN4(c-dotat);      // source
  1368        t=s+1;                            // target
  1369        // open the gap [cannot use memcpy]
  1370        for (; s>=dotat; s-=4, t-=4) UBFROMUI(t, UBTOUI(s));
  1371        *dotat='.';
  1372        c++;                              // length increased by one
  1373        } // need dot?
  1374       else for (; c<dotat; c++) *c='0';  // pad for engineering
  1375      } // pre>0
  1376     else {
  1377      /* -5<=pre<=0: here for plain 0.ddd or 0.000ddd forms (may have
  1378         E, but only for 0.00E+3 kind of case -- with plenty of spare
  1379         space in this case */
  1380      pre=-pre+2;                         // gap width, including "0."
  1381      t=cstart+ROUNDDOWN4(c-cstart)+pre;  // preferred first target point
  1382      // backoff if too far to the right
  1383      if (t>string+DECSTRING-5) t=string+DECSTRING-5; // adjust to fit
  1384      // now shift the entire coefficient to the right, being careful not
  1385      // to access to the left of string [cannot use memcpy]
  1386      for (s=t-pre; s>=string; s-=4, t-=4) UBFROMUI(t, UBTOUI(s));
  1387      // for Quads and Singles there may be a character or two left...
  1388      s+=3;                               // where next would come from
  1389      for(; s>=cstart; s--, t--) *(t+3)=*(s);
  1390      // now have fill 0. through 0.00000; use overlaps to avoid tests
  1391      if (pre>=4) {
  1392        memcpy(cstart+pre-4, "0000", 4);
  1393        memcpy(cstart, "0.00", 4);
  1394        }
  1395       else { // 2 or 3
  1396        *(cstart+pre-1)='0';
  1397        memcpy(cstart, "0.", 2);
  1398        }
  1399      c+=pre;                             // to end
  1400      }
  1401  
  1402    // finally add the E-part, if needed; it will never be 0, and has
  1403    // a maximum length of 3 or 4 digits (asserted above)
  1404    if (e!=0) {
  1405      memcpy(c, "E+", 2);                 // starts with E, assume +
  1406      c++;
  1407      if (e<0) {
  1408        *c='-';                           // oops, need '-'
  1409        e=-e;                             // uInt, please
  1410        }
  1411      c++;
  1412      // Three-character exponents are easy; 4-character a little trickier
  1413      #if DECEMAXD<=3
  1414        u=&BIN2BCD8[e*4];                 // -> 3 digits + length byte
  1415        // copy fixed 4 characters [is safe], starting at non-zero
  1416        // and with character mask to convert BCD to char
  1417        UBFROMUI(c, UBTOUI(u+3-*(u+3))|CHARMASK);
  1418        c+=*(u+3);                        // bump pointer appropriately
  1419      #elif DECEMAXD==4
  1420        if (e<1000) {                     // 3 (or fewer) digits case
  1421          u=&BIN2BCD8[e*4];               // -> 3 digits + length byte
  1422          UBFROMUI(c, UBTOUI(u+3-*(u+3))|CHARMASK); // [as above]
  1423          c+=*(u+3);                      // bump pointer appropriately
  1424          }
  1425         else {                           // 4-digits
  1426          Int thou=((e>>3)*1049)>>17;     // e/1000
  1427          Int rem=e-(1000*thou);          // e%1000
  1428          *c++=(char)('0'+(char)thou);    // the thousands digit
  1429          u=&BIN2BCD8[rem*4];             // -> 3 digits + length byte
  1430          UBFROMUI(c, UBTOUI(u)|CHARMASK);// copy fixed 3+1 characters [is safe]
  1431          c+=3;                           // bump pointer, always 3 digits
  1432          }
  1433      #endif
  1434      }
  1435    *c='\0';                              // terminate
  1436    //printf("res %s\n", string);
  1437    return string;
  1438    } // decFloatToEngString
  1439  
  1440  /* ------------------------------------------------------------------ */
  1441  /* decFloatToPacked -- convert decFloat to Packed decimal + exponent  */
  1442  /*                                                                    */
  1443  /*  df is the source decFloat                                         */
  1444  /*  exp will be set to the unbiased exponent, q, or to a special      */
  1445  /*    value in the form returned by decFloatGetExponent               */
  1446  /*  packed is where DECPMAX nibbles will be written with the sign as  */
  1447  /*    final nibble (0x0c for +, 0x0d for -); a NaN has a first nibble */
  1448  /*    of zero, and an infinity is all zeros. decDouble and decQuad    */
  1449  /*    have a additional leading zero nibble, leading to result        */
  1450  /*    lengths of 4, 9, and 18 bytes.                                  */
  1451  /*  returns the sign of the coefficient (DECFLOAT_Sign if negative,   */
  1452  /*    0 otherwise)                                                    */
  1453  /*                                                                    */
  1454  /* No error is possible, and no status will be set.                   */
  1455  /* ------------------------------------------------------------------ */
  1456  Int decFloatToPacked(const decFloat *df, Int *exp, uByte *packed) {
  1457    uByte bcdar[DECPMAX+2];          // work buffer
  1458    uByte *ip=bcdar, *op=packed;     // work pointers
  1459    if (DFISINF(df)) {
  1460      memset(bcdar, 0, DECPMAX+2);
  1461      *exp=DECFLOAT_Inf;
  1462      }
  1463     else {
  1464      GETCOEFF(df, bcdar+1);         // use macro
  1465      if (DFISNAN(df)) {
  1466        bcdar[1]=0;                  // MSD needs clearing
  1467        *exp=DFWORD(df, 0)&0x7e000000;
  1468        }
  1469       else {                        // finite
  1470        *exp=GETEXPUN(df);
  1471        }
  1472      }
  1473    // now pack; coefficient currently at bcdar+1
  1474    #if SINGLE
  1475      ip++;                          // ignore first byte
  1476    #else
  1477      *ip=0;                         // need leading zero
  1478    #endif
  1479    // set final byte to Packed BCD sign value
  1480    bcdar[DECPMAX+1]=(DFISSIGNED(df) ? DECPMINUS : DECPPLUS);
  1481    // pack an even number of bytes...
  1482    for (; op<packed+((DECPMAX+2)/2); op++, ip+=2) {
  1483      *op=(uByte)((*ip<<4)+*(ip+1));
  1484      }
  1485    return (bcdar[DECPMAX+1]==DECPMINUS ? DECFLOAT_Sign : 0);
  1486    } // decFloatToPacked
  1487  
  1488  /* ------------------------------------------------------------------ */
  1489  /* decFloatToString -- conversion to numeric string                   */
  1490  /*                                                                    */
  1491  /*  df is the decFloat format number to convert                       */
  1492  /*  string is the string where the result will be laid out            */
  1493  /*                                                                    */
  1494  /* string must be at least DECPMAX+9 characters (the worst case is    */
  1495  /* "-0.00000nnn...nnn\0", which is as long as the exponent form when  */
  1496  /* DECEMAXD<=4); this condition is asserted above                     */
  1497  /*                                                                    */
  1498  /* No error is possible, and no status will be set                    */
  1499  /* ------------------------------------------------------------------ */
  1500  char * decFloatToString(const decFloat *df, char *string){
  1501    uInt msd;                        // coefficient MSD
  1502    Int  exp;                        // exponent top two bits or full
  1503    uInt comb;                       // combination field
  1504    char *cstart;                    // coefficient start
  1505    char *c;                         // output pointer in string
  1506    char *s, *t;                     // .. (source, target)
  1507    Int  pre, e;                     // work
  1508    const uByte *u;                  // ..
  1509    uInt  uiwork;                    // for macros [one compiler needs
  1510                                     // volatile here to avoid bug, but
  1511                                     // that doubles execution time]
  1512  
  1513    // Source words; macro handles endianness
  1514    uInt sourhi=DFWORD(df, 0);       // word with sign
  1515    #if DECPMAX==16
  1516    uInt sourlo=DFWORD(df, 1);
  1517    #elif DECPMAX==34
  1518    uInt sourmh=DFWORD(df, 1);
  1519    uInt sourml=DFWORD(df, 2);
  1520    uInt sourlo=DFWORD(df, 3);
  1521    #endif
  1522  
  1523    c=string;                        // where result will go
  1524    if (((Int)sourhi)<0) *c++='-';   // handle sign
  1525    comb=sourhi>>26;                 // sign+combination field
  1526    msd=DECCOMBMSD[comb];            // decode the combination field
  1527    exp=DECCOMBEXP[comb];            // ..
  1528  
  1529    if (!EXPISSPECIAL(exp)) {        // finite
  1530      // complete exponent; top two bits are in place
  1531      exp+=GETECON(df)-DECBIAS;      // .. + continuation and unbias
  1532      }
  1533     else {                          // IS special
  1534      if (exp==DECFLOAT_Inf) {       // infinity
  1535        strcpy(c, "Infinity");
  1536        return string;               // easy
  1537        }
  1538      if (sourhi&0x02000000) *c++='s'; // sNaN
  1539      strcpy(c, "NaN");              // complete word
  1540      c+=3;                          // step past
  1541      // quick exit if the payload is zero
  1542      #if DECPMAX==7
  1543      if ((sourhi&0x000fffff)==0) return string;
  1544      #elif DECPMAX==16
  1545      if (sourlo==0 && (sourhi&0x0003ffff)==0) return string;
  1546      #elif DECPMAX==34
  1547      if (sourlo==0 && sourml==0 && sourmh==0
  1548       && (sourhi&0x00003fff)==0) return string;
  1549      #endif
  1550      // otherwise drop through to add integer; set correct exp etc.
  1551      exp=0; msd=0;                  // setup for following code
  1552      }
  1553  
  1554    /* convert the digits of the significand to characters */
  1555    cstart=c;                        // save start of coefficient
  1556    if (msd) *c++=(char)('0'+(char)msd);  // non-zero most significant digit
  1557  
  1558    // Decode the declets.  After extracting each declet, it is
  1559    // decoded to a 4-uByte sequence by table lookup; the four uBytes
  1560    // are the three encoded BCD8 digits followed by a 1-byte length
  1561    // (significant digits, except that 000 has length 0).  This allows
  1562    // us to left-align the first declet with non-zero content, then
  1563    // the remaining ones are full 3-char length.  Fixed-length copies
  1564    // are used because variable-length memcpy causes a subroutine call
  1565    // in at least two compilers.  (The copies are length 4 for speed
  1566    // and are safe because the last item in the array is of length
  1567    // three and has the length byte following.)
  1568    #define dpd2char(dpdin) u=&DPD2BCD8[((dpdin)&0x3ff)*4];        \
  1569           if (c!=cstart) {UBFROMUI(c, UBTOUI(u)|CHARMASK); c+=3;} \
  1570            else if (*(u+3)) {                                     \
  1571             UBFROMUI(c, UBTOUI(u+3-*(u+3))|CHARMASK); c+=*(u+3);}
  1572  
  1573    #if DECPMAX==7
  1574    dpd2char(sourhi>>10);                 // declet 1
  1575    dpd2char(sourhi);                     // declet 2
  1576  
  1577    #elif DECPMAX==16
  1578    dpd2char(sourhi>>8);                  // declet 1
  1579    dpd2char((sourhi<<2) | (sourlo>>30)); // declet 2
  1580    dpd2char(sourlo>>20);                 // declet 3
  1581    dpd2char(sourlo>>10);                 // declet 4
  1582    dpd2char(sourlo);                     // declet 5
  1583  
  1584    #elif DECPMAX==34
  1585    dpd2char(sourhi>>4);                  // declet 1
  1586    dpd2char((sourhi<<6) | (sourmh>>26)); // declet 2
  1587    dpd2char(sourmh>>16);                 // declet 3
  1588    dpd2char(sourmh>>6);                  // declet 4
  1589    dpd2char((sourmh<<4) | (sourml>>28)); // declet 5
  1590    dpd2char(sourml>>18);                 // declet 6
  1591    dpd2char(sourml>>8);                  // declet 7
  1592    dpd2char((sourml<<2) | (sourlo>>30)); // declet 8
  1593    dpd2char(sourlo>>20);                 // declet 9
  1594    dpd2char(sourlo>>10);                 // declet 10
  1595    dpd2char(sourlo);                     // declet 11
  1596    #endif
  1597  
  1598    if (c==cstart) *c++='0';         // all zeros, empty -- make "0"
  1599  
  1600    //[This fast path is valid but adds 3-5 cycles to worst case length]
  1601    //if (exp==0) {                  // integer or NaN case -- easy
  1602    //  *c='\0';                     // terminate
  1603    //  return string;
  1604    //  }
  1605  
  1606    e=0;                             // assume no E
  1607    pre=(Int)(c-cstart)+exp;         // length+exp  [c->LSD+1]
  1608    // [here, pre-exp is the digits count (==1 for zero)]
  1609  
  1610    if (exp>0 || pre<-5) {           // need exponential form
  1611      e=pre-1;                       // calculate E value
  1612      pre=1;                         // assume one digit before '.'
  1613      } // exponential form
  1614  
  1615    /* modify the coefficient, adding 0s, '.', and E+nn as needed */
  1616    if (pre>0) {                     // ddd.ddd (plain), perhaps with E
  1617      char *dotat=cstart+pre;
  1618      if (dotat<c) {                      // if embedded dot needed...
  1619        // [memmove is a disaster, here]
  1620        // move by fours; there must be space for junk at the end
  1621        // because exponent is still possible
  1622        s=dotat+ROUNDDOWN4(c-dotat);      // source
  1623        t=s+1;                            // target
  1624        // open the gap [cannot use memcpy]
  1625        for (; s>=dotat; s-=4, t-=4) UBFROMUI(t, UBTOUI(s));
  1626        *dotat='.';
  1627        c++;                              // length increased by one
  1628        } // need dot?
  1629  
  1630      // finally add the E-part, if needed; it will never be 0, and has
  1631      // a maximum length of 3 or 4 digits (asserted above)
  1632      if (e!=0) {
  1633        memcpy(c, "E+", 2);               // starts with E, assume +
  1634        c++;
  1635        if (e<0) {
  1636          *c='-';                         // oops, need '-'
  1637          e=-e;                           // uInt, please
  1638          }
  1639        c++;
  1640        // Three-character exponents are easy; 4-character a little trickier
  1641        #if DECEMAXD<=3
  1642          u=&BIN2BCD8[e*4];               // -> 3 digits + length byte
  1643          // copy fixed 4 characters [is safe], starting at non-zero
  1644          // and with character mask to convert BCD to char
  1645          UBFROMUI(c, UBTOUI(u+3-*(u+3))|CHARMASK);
  1646          c+=*(u+3);                      // bump pointer appropriately
  1647        #elif DECEMAXD==4
  1648          if (e<1000) {                   // 3 (or fewer) digits case
  1649            u=&BIN2BCD8[e*4];             // -> 3 digits + length byte
  1650            UBFROMUI(c, UBTOUI(u+3-*(u+3))|CHARMASK); // [as above]
  1651            c+=*(u+3);                    // bump pointer appropriately
  1652            }
  1653           else {                         // 4-digits
  1654            Int thou=((e>>3)*1049)>>17;   // e/1000
  1655            Int rem=e-(1000*thou);        // e%1000
  1656            *c++=(char)('0'+(char)thou);  // the thousands digit
  1657            u=&BIN2BCD8[rem*4];           // -> 3 digits + length byte
  1658            UBFROMUI(c, UBTOUI(u)|CHARMASK); // copy fixed 3+1 characters [is safe]
  1659            c+=3;                         // bump pointer, always 3 digits
  1660            }
  1661        #endif
  1662        }
  1663      *c='\0';                            // add terminator
  1664      //printf("res %s\n", string);
  1665      return string;
  1666      } // pre>0
  1667  
  1668    /* -5<=pre<=0: here for plain 0.ddd or 0.000ddd forms (can never have E) */
  1669    // Surprisingly, this is close to being the worst-case path, so the
  1670    // shift is done by fours; this is a little tricky because the
  1671    // rightmost character to be written must not be beyond where the
  1672    // rightmost terminator could be -- so backoff to not touch
  1673    // terminator position if need be (this can make exact alignments
  1674    // for full Doubles, but in some cases needs care not to access too
  1675    // far to the left)
  1676  
  1677    pre=-pre+2;                           // gap width, including "0."
  1678    t=cstart+ROUNDDOWN4(c-cstart)+pre;    // preferred first target point
  1679    // backoff if too far to the right
  1680    if (t>string+DECSTRING-5) t=string+DECSTRING-5; // adjust to fit
  1681    // now shift the entire coefficient to the right, being careful not
  1682    // to access to the left of string [cannot use memcpy]
  1683    for (s=t-pre; s>=string; s-=4, t-=4) UBFROMUI(t, UBTOUI(s));
  1684    // for Quads and Singles there may be a character or two left...
  1685    s+=3;                                 // where next would come from
  1686    for(; s>=cstart; s--, t--) *(t+3)=*(s);
  1687    // now have fill 0. through 0.00000; use overlaps to avoid tests
  1688    if (pre>=4) {
  1689      memcpy(cstart+pre-4, "0000", 4);
  1690      memcpy(cstart, "0.00", 4);
  1691      }
  1692     else { // 2 or 3
  1693      *(cstart+pre-1)='0';
  1694      memcpy(cstart, "0.", 2);
  1695      }
  1696    *(c+pre)='\0';                        // terminate
  1697    return string;
  1698    } // decFloatToString
  1699  
  1700  /* ------------------------------------------------------------------ */
  1701  /* decFloatToWider -- conversion to next-wider format                 */
  1702  /*                                                                    */
  1703  /*  source  is the decFloat format number which gets the result of    */
  1704  /*          the conversion                                            */
  1705  /*  wider   is the decFloatWider format number which will be narrowed */
  1706  /*  returns wider                                                     */
  1707  /*                                                                    */
  1708  /* Widening is always exact; no status is set (sNaNs are copied and   */
  1709  /* do not signal).  The result will be canonical if the source is,    */
  1710  /* and may or may not be if the source is not.                        */
  1711  /* ------------------------------------------------------------------ */
  1712  // widening is not possible for decQuad format numbers; simply omit
  1713  #if !QUAD
  1714  decFloatWider * decFloatToWider(const decFloat *source, decFloatWider *wider) {
  1715    uInt msd;
  1716  
  1717    /* Construct and copy the sign word */
  1718    if (DFISSPECIAL(source)) {
  1719      // copy sign, combination, and first bit of exponent (sNaN selector)
  1720      DFWWORD(wider, 0)=DFWORD(source, 0)&0xfe000000;
  1721      msd=0;
  1722      }
  1723     else { // is finite number
  1724      uInt exp=GETEXPUN(source)+DECWBIAS; // get unbiased exponent and rebias
  1725      uInt code=(exp>>DECWECONL)<<29;     // set two bits of exp [msd=0]
  1726      code|=(exp<<(32-6-DECWECONL)) & 0x03ffffff; // add exponent continuation
  1727      code|=DFWORD(source, 0)&0x80000000; // add sign
  1728      DFWWORD(wider, 0)=code;             // .. and place top word in wider
  1729      msd=GETMSD(source);                 // get source coefficient MSD [0-9]
  1730      }
  1731    /* Copy the coefficient and clear any 'unused' words to left */
  1732    #if SINGLE
  1733      DFWWORD(wider, 1)=(DFWORD(source, 0)&0x000fffff)|(msd<<20);
  1734    #elif DOUBLE
  1735      DFWWORD(wider, 2)=(DFWORD(source, 0)&0x0003ffff)|(msd<<18);
  1736      DFWWORD(wider, 3)=DFWORD(source, 1);
  1737      DFWWORD(wider, 1)=0;
  1738    #endif
  1739    return wider;
  1740    } // decFloatToWider
  1741  #endif
  1742  
  1743  /* ------------------------------------------------------------------ */
  1744  /* decFloatVersion -- return package version string                   */
  1745  /*                                                                    */
  1746  /*  returns a constant string describing this package                 */
  1747  /* ------------------------------------------------------------------ */
  1748  const char *decFloatVersion(void) {
  1749    return DECVERSION;
  1750    } // decFloatVersion
  1751  
  1752  /* ------------------------------------------------------------------ */
  1753  /* decFloatZero -- set to canonical (integer) zero                    */
  1754  /*                                                                    */
  1755  /*  df is the decFloat format number to integer +0 (q=0, c=+0)        */
  1756  /*  returns df                                                        */
  1757  /*                                                                    */
  1758  /* No error is possible, and no status can be set.                    */
  1759  /* ------------------------------------------------------------------ */
  1760  decFloat * decFloatZero(decFloat *df){
  1761    DFWORD(df, 0)=ZEROWORD;     // set appropriate top word
  1762    #if DOUBLE || QUAD
  1763      DFWORD(df, 1)=0;
  1764      #if QUAD
  1765        DFWORD(df, 2)=0;
  1766        DFWORD(df, 3)=0;
  1767      #endif
  1768    #endif
  1769    // decFloatShow(df, "zero");
  1770    return df;
  1771    } // decFloatZero
  1772  
  1773  /* ------------------------------------------------------------------ */
  1774  /* Private generic function (not format-specific) for development use */
  1775  /* ------------------------------------------------------------------ */
  1776  // This is included once only, for all to use
  1777  #if QUAD && (DECCHECK || DECTRACE)
  1778    /* ---------------------------------------------------------------- */
  1779    /* decShowNum -- display bcd8 number in debug form                  */
  1780    /*                                                                  */
  1781    /*   num is the bcdnum to display                                   */
  1782    /*   tag is a string to label the display                           */
  1783    /* ---------------------------------------------------------------- */
  1784    void decShowNum(const bcdnum *num, const char *tag) {
  1785      const char *csign="+";              // sign character
  1786      uByte *ub;                          // work
  1787      uInt  uiwork;                       // for macros
  1788      if (num->sign==DECFLOAT_Sign) csign="-";
  1789  
  1790      printf(">%s> ", tag);
  1791      if (num->exponent==DECFLOAT_Inf) printf("%sInfinity", csign);
  1792      else if (num->exponent==DECFLOAT_qNaN) printf("%sqNaN", csign);
  1793      else if (num->exponent==DECFLOAT_sNaN) printf("%ssNaN", csign);
  1794      else {                              // finite
  1795       char qbuf[10];                     // for right-aligned q
  1796       char *c;                           // work
  1797       const uByte *u;                    // ..
  1798       Int e=num->exponent;               // .. exponent
  1799       strcpy(qbuf, "q=");
  1800       c=&qbuf[2];                        // where exponent will go
  1801       // lay out the exponent
  1802       if (e<0) {
  1803         *c++='-';                        // add '-'
  1804         e=-e;                            // uInt, please
  1805         }
  1806       #if DECEMAXD>4
  1807         #error Exponent form is too long for ShowNum to lay out
  1808       #endif
  1809       if (e==0) *c++='0';                // 0-length case
  1810        else if (e<1000) {                // 3 (or fewer) digits case
  1811         u=&BIN2BCD8[e*4];                // -> 3 digits + length byte
  1812         UBFROMUI(c, UBTOUI(u+3-*(u+3))|CHARMASK); // [as above]
  1813         c+=*(u+3);                       // bump pointer appropriately
  1814         }
  1815        else {                            // 4-digits
  1816         Int thou=((e>>3)*1049)>>17;      // e/1000
  1817         Int rem=e-(1000*thou);           // e%1000
  1818         *c++=(char)('0'+(char)thou);     // the thousands digit
  1819         u=&BIN2BCD8[rem*4];              // -> 3 digits + length byte
  1820         UBFROMUI(c, UBTOUI(u)|CHARMASK); // copy fixed 3+1 characters [is safe]
  1821         c+=3;                            // bump pointer, always 3 digits
  1822         }
  1823       *c='\0';                           // add terminator
  1824       printf("%7s c=%s", qbuf, csign);
  1825       }
  1826  
  1827      if (!EXPISSPECIAL(num->exponent) || num->msd!=num->lsd || *num->lsd!=0) {
  1828        for (ub=num->msd; ub<=num->lsd; ub++) { // coefficient...
  1829          printf("%1x", *ub);
  1830          if ((num->lsd-ub)%3==0 && ub!=num->lsd) printf(" "); // 4-space
  1831          }
  1832        }
  1833      printf("\n");
  1834      } // decShowNum
  1835  #endif