Open64: osprey/common/util/c_q_to

00001 /*
00002  * Copyright (C) 2007. PathScale, LLC. All Rights Reserved.
00003  */
00004 /*
00005  * Copyright 2004, 2005, 2006 PathScale, Inc.  All Rights Reserved.
00006  */
00007 
00008 /*
00009 
00010   Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
00011 
00012   This program is free software; you can redistribute it and/or modify it
00013   under the terms of version 2 of the GNU General Public License as
00014   published by the Free Software Foundation.
00015 
00016   This program is distributed in the hope that it would be useful, but
00017   WITHOUT ANY WARRANTY; without even the implied warranty of
00018   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
00019 
00020   Further, this software is distributed without any warranty that it is
00021   free of the rightful claim of any third person regarding infringement 
00022   or the like.  Any license provided herein, whether implied or 
00023   otherwise, applies only to this software file.  Patent licenses, if 
00024   any, provided herein do not apply to combinations of this program with 
00025   other software, or any other product whatsoever.  
00026 
00027   You should have received a copy of the GNU General Public License along
00028   with this program; if not, write the Free Software Foundation, Inc., 59
00029   Temple Place - Suite 330, Boston MA 02111-1307, USA.
00030 
00031   Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
00032   Mountain View, CA 94043, or:
00033 
00034   http://www.sgi.com
00035 
00036   For further information regarding this notice, see:
00037 
00038   http://oss.sgi.com/projects/GenInfo/NoticeExplan
00039 
00040 */
00041 
00042 
00043 /* ====================================================================
00044  * ====================================================================
00045  *
00046  * Module: c_q_to_a
00047  * $Revision: 1.6 $
00048  * $Date: 04/12/21 14:57:27-08:00 $
00049  * $Author: bos@eng-25.internal.keyresearch.com $
00050  * $Source: /home/bos/bk/kpro64-pending/common/util/SCCS/s.c_q_to_a.c $
00051  *
00052  * Revision history:
00053  *  26-jul-93 - Original Version
00054  *
00055  * Description: Convert a quad precision binary number to decimal.
00056  *
00057  * ====================================================================
00058  * ====================================================================
00059  */
00060 
00061 
00062 #if QUAD_DEBUG
00063 #include <stdio.h>
00064 #endif
00065 
00066 /* Macros to pull apart parts of single and  double precision
00067  * floating point numbers in IEEE format
00068  */
00069 
00070 /* double precision */
00071 
00072 typedef  union {
00073   struct {
00074 #if defined(_LITTLE_ENDIAN)
00075     unsigned int  lo  :32;
00076     unsigned int  hi  :20;
00077     unsigned int  exp :11;
00078     unsigned int  sign  :1;
00079 #else
00080     unsigned int  sign  :1;
00081     unsigned int  exp :11;
00082     unsigned int  hi  :20;
00083     unsigned int  lo  :32;
00084 #endif
00085   } fparts;
00086   struct {
00087 #if defined(_LITTLE_ENDIAN)
00088     unsigned int  lo  :32;
00089     unsigned int  hi  :19;
00090     unsigned int  qnan_bit  :1;
00091     unsigned int  exp :11;
00092     unsigned int  sign  :1;
00093 #else
00094     unsigned int  sign  :1;
00095     unsigned int  exp :11;
00096     unsigned int  qnan_bit  :1;
00097     unsigned int  hi  :19;
00098     unsigned int  lo  :32;
00099 #endif
00100   } nparts;
00101   struct {
00102 #if defined(_LITTLE_ENDIAN)
00103     unsigned int lo;
00104     unsigned int hi;
00105 #else
00106     unsigned int hi;
00107     unsigned int lo;
00108 #endif
00109   } fwords;
00110   double  d;
00111 } _dval;
00112 
00113 /* quad precision */
00114 
00115 typedef  union {
00116   struct {
00117     _dval hi;
00118     _dval lo;
00119   } qparts;
00120   struct {
00121     unsigned long long hi;
00122     unsigned long long lo;
00123   } fwords;
00124 } _ldblval;
00125 
00126 /* parts of a double precision floating point number */
00127 
00128 #define SIGNBIT(X)  (((_dval *)&(X))->fparts.sign)
00129 #define EXPONENT(X) (((_dval *)&(X))->fparts.exp)
00130 
00131 #include "defs.h"
00132 #include "quad.h"
00133 
00134 
00135 #define min(x,y) ((x)<(y)? (x): (y))
00136 #define max(x,y) ((x)>(y)? (x): (y))
00137 
00138 extern void c_qtenscale( _ldblval *, INT32, INT32 *);
00139 static INT32  c_qtoa(char *buffer, INT32 ndigit, QUAD x, INT32 fflag);
00140 extern INT  c_q_to_a(char *str, QUAD x, INT *p_err );
00141 
00142 #if defined(BUILD_OS_DARWIN)
00143 /* Can't use "pragma weak" to create aliases in Mach-O */
00144 INT c_q_to_a(char *str, QUAD x, INT *p_err );
00145 INT __c_q_to_a(char *str, QUAD x, INT *p_err )
00146   { return c_q_to_a(str, x, p_err); }
00147 #else /* defined(BUILD_OS_DARWIN) */
00148 #pragma weak c_q_to_a = __c_q_to_a
00149 #define c_q_to_a __c_q_to_a
00150 #endif /* defined(BUILD_OS_DARWIN) */
00151 
00152 
00153 /* ====================================================================
00154  *
00155  * FunctionName: c_q_to_a
00156  *
00157  * Description: Convert a quad precision binary number to decimal.
00158  *
00159  * ====================================================================
00160  */
00161 
00162 
00163 INT
00164 c_q_to_a(char *str, QUAD x, INT *p_err )
00165 {
00166 char  buffer[42];
00167 char  *pbuff;
00168 INT32 ndigit;
00169 INT32 fflag;
00170 INT32 m, n;
00171 
00172   *p_err = 0;
00173 
00174   ndigit = 34;
00175   fflag = 0;
00176 
00177   n = c_qtoa(buffer, ndigit, x, fflag);
00178 
00179   if ( strncmp(buffer, "inf", 3) == 0 )
00180   {
00181     strncpy( str, buffer, 3 );
00182 
00183     *(str + 3) = '\0';
00184 
00185     return ( *p_err );
00186   }
00187 
00188   if ( strncmp(buffer, "nan", 3) == 0 )
00189   {
00190     strncpy( str, buffer, 3 );
00191 
00192     *(str + 3) = '\0';
00193 
00194     return ( *p_err );
00195   }
00196 
00197   if ( strncmp(buffer+1, "inf", 3) == 0 )
00198   {
00199     strncpy( str, buffer, 4 );
00200 
00201     *(str + 4) = '\0';
00202 
00203     return ( *p_err );
00204   }
00205 
00206   if ( strncmp(buffer+1, "nan", 3) == 0 )
00207   {
00208     strncpy( str, buffer+1, 3 );
00209 
00210     *(str + 3) = '\0';
00211 
00212     return ( *p_err );
00213   }
00214 
00215   if ( (n <= -1000) || (n >= 1000) )
00216   {
00217     *p_err = 1;
00218 
00219     *str = '\0';
00220 
00221     return ( *p_err );
00222   }
00223 
00224   pbuff = buffer;
00225 
00226   if ( (buffer[0] == ' ') || (buffer[0] == '-') )
00227     *str++ = *pbuff++;
00228 
00229   *str++ = *pbuff++;  /* copy first digit to output string */
00230 
00231   *str++ = '.';
00232 
00233   if ( ndigit > 1 )
00234   {
00235     strncpy(str, pbuff, ndigit-1);
00236 
00237     str += ndigit - 1;  /* skip over digits copied */
00238   }
00239 
00240   *str++ = 'E';
00241 
00242   if ( n < 0 )
00243   {
00244     n = -n;
00245     *str++ = '-';
00246   }
00247 
00248   m = n/100;
00249 
00250   if ( m!= 0 )
00251     *str++ = m + '0';
00252 
00253   n = n%100;
00254 
00255   m = n/10;
00256 
00257   *str++ = m + '0';
00258 
00259   m = n%10;
00260 
00261   *str++ = m + '0';
00262 
00263   *str++ = '\0';
00264 
00265   return ( *p_err );
00266 }
00267 
00268 
00269 /* ====================================================================
00270  *
00271  * FunctionName: c_qtoa
00272  *
00273  * Description: Convert a quad precision binary number to decimal.
00274  *
00275  * Except for minor changes, this routine is the same as _qtoa()
00276  * in libc.  Included here, so we can do cross compilations on
00277  * machines not supporting quad precision.
00278  *
00279  * ====================================================================
00280  */
00281 
00282 static
00283 INT32       /* returns decimal exponent */
00284 c_qtoa(buffer,ndigit,x,fflag) 
00285      char *buffer;    /* Destination for ascii string */ 
00286      INT32 ndigit;    /* number of digits to convert.
00287            * 1 <= ndigits <= 34
00288          * Log2(1e34) = 112.9 bits, so the
00289          * 107 significant bits can be con-
00290          * verted to decimal with room left
00291          * over.
00292          */ 
00293      QUAD x;      /* number to convert to ascii */ 
00294      INT32 fflag;   /* 0 for e format, 1 for f format */ 
00295 {
00296   _ldblval value;   /* the value manipulated */ 
00297 #define HI (value.fwords.hi)  /* top 64b of 128b fraction */ 
00298 #define LO (value.fwords.lo)  /* bottom 64b of 128b fraction */
00299 
00300   INT32 signhi, signlo;   /* sign of HI = sign of number */ 
00301   INT32 exphi, explo;   /* exponents of HI and LO */ 
00302   INT32 dexp;     /* decimal exponent */ 
00303   INT32 bexp;     /* binary exponent */ 
00304   UINT32 lohi;      /* top 4 bits of lo */ 
00305   UINT64 lolo;      /* bottom 60 bits of lo */ 
00306   INT32 nd;     /* number of digits output */ 
00307   INT32 guard;      /* guard digit */
00308 
00309 #if QUAD_DEBUG 
00310   char *buffer0=buffer;   /* initial value of buffer */
00311 #endif
00312 
00313 
00314   value.qparts.hi.d = x.hi;
00315   value.qparts.lo.d = x.lo;
00316 
00317 #if QUAD_DEBUG 
00318   printf("initially: HI=0x%016llx\n",HI); 
00319   printf("LO=0x%016llx\n",LO); 
00320 #endif
00321 
00322   signhi = SIGNBIT(HI); 
00323   HI = (HI<<1)>>1;    /* discard sign bit */ 
00324   *buffer++ = signhi? '-': ' ';
00325 
00326 #if QUAD_DEBUG 
00327   printf("after de-signing hi: HI=0x%016llx\n",HI); 
00328   printf("signhi=%d\n",signhi); 
00329   *buffer = '\0'; 
00330   printf(" buffer=%s\n",buffer0);
00331 #endif
00332 
00333   if( HI == 0 ) {   /* true zero */ 
00334     char *b;      /* fencepost for buffer */
00335 
00336     b = buffer + ndigit; 
00337     while(buffer < b) 
00338       *buffer++ = '0';
00339 
00340 #if QUAD_DEBUG 
00341     printf("HI was zero: HI=0x%016llx\n",HI); 
00342     printf("signhi=%d\n",signhi); 
00343     *buffer = '\0'; 
00344     printf(" buffer=%s\n",buffer0);
00345 #endif
00346 
00347     return(0);
00348   }
00349 
00350   exphi = EXPONENT(HI); 
00351   EXPONENT(HI) = 0;
00352 
00353 #if QUAD_DEBUG 
00354   printf("HI non-zero: HI=0x%016llx\n",HI); 
00355   printf("exphi=%d\n",exphi); 
00356 #endif
00357 
00358   if(exphi == 2047) {   /* infinity or NaN */ 
00359     char *s;      /* infinity or nan string */ 
00360     char *infinity = "inf"; 
00361     char *nan = "nan";
00362 
00363     s = HI==0? infinity: nan;
00364 
00365     do 
00366       *buffer++ = *s; 
00367     while(*s++);
00368 
00369 #if QUAD_DEBUG 
00370     printf("NaN or Infinity: buffer=%s\n",buffer0); 
00371 #endif
00372 
00373     return(0); 
00374   }
00375 
00376   /* processing a number. create 128b fraction */ 
00377   /* hi is non-zero, or the following while loop wouldn't terminate */ 
00378   if(exphi == 0) {    /* HI denorm */ 
00379     exphi++;
00380     while( !(HI & (1ULL<<52))) { /* normalize - need not be fast */ 
00381       HI <<=1; 
00382       exphi--;
00383 
00384 #if QUAD_DEBUG 
00385       printf("HI denorm: HI=0x%016llx\n",HI); 
00386       printf("exphi=%d\n",exphi); 
00387 #endif
00388 
00389     }
00390   } 
00391   else {      /* HI is normal number */ 
00392     HI |= 1ULL<<52;   /* expose hidden bit */
00393 
00394 #if QUAD_DEBUG 
00395     printf("HI normal: HI=0x%016llx\n",HI); 
00396     printf("exphi=%d\n",exphi); 
00397 #endif
00398 
00399   }
00400 
00401   HI <<=11;     /* left align */
00402 
00403 #if QUAD_DEBUG 
00404   printf("HI aligned: HI=0x%016llx\n",HI); 
00405 #endif
00406 
00407   /* handle LO */
00408 
00409   signlo = SIGNBIT(LO) ^ signhi; /* HI,LO signs differ? */
00410   SIGNBIT(LO) = 0;    /* discard sign */
00411 
00412 #if QUAD_DEBUG
00413   printf("handle LO: LO=0x%016llx\n",LO);
00414   printf("           signlo=%d\n",signlo);
00415 #endif
00416 
00417   if(LO) {      /* lo is non-zero */
00418 
00419     explo = EXPONENT(LO);
00420     EXPONENT(LO) = 0;
00421 
00422 #if QUAD_DEBUG
00423   printf("non-zero LO: LO=0x%016llx\n",LO);
00424   printf("             explo=%d\n",explo);
00425 #endif
00426 
00427     LO <<=1;      /* shift LO to expose 107th bit */
00428 
00429 #if QUAD_DEBUG
00430     printf("left-shifted LO: LO=0x%016llx\n",LO);
00431 #endif
00432 
00433     if(explo == 0) {    /* denorm */
00434       INT32 shift_n;    /* shift count */
00435 
00436       LO <<= 1;     /* shift up one more bit */
00437 #if QUAD_DEBUG
00438       printf("shifted up LO: LO=0x%016llx\n",LO);
00439 #endif
00440       shift_n = 53 - (exphi-explo);
00441       /* normalize to HI */
00442       if(shift_n > 0) {
00443   LO <<= shift_n; 
00444       }
00445       else if(shift_n < 0) {
00446   LO >>= -shift_n;
00447       }
00448     }
00449     else {      /* LO normal number */
00450       LO |= 1ULL<<53;   /* expose hidden bit */
00451       
00452 #if QUAD_DEBUG
00453       printf("expose hidden bit: LO=0x%016llx\n",LO);
00454 #endif
00455       
00456       LO >>= (exphi-explo) - 53; /* normalize to HI */
00457     }
00458 
00459 #if QUAD_DEBUG
00460     printf("normalize LO to HI: LO=0x%016llx\n",LO);
00461 #endif
00462     
00463     if( signlo ) {    /* LO is negative */
00464       HI-= 1<<11;     /* unround HI */
00465       LO = (1ULL<<54) - LO; /* complement LO */
00466       
00467 #if QUAD_DEBUG
00468       printf("LO negative: HI=0x%016llx\n",HI);
00469       printf("             LO=0x%016llx\n",LO);
00470 #endif
00471       
00472     }
00473 
00474     HI |= LO>>43;   /* combine HI and LO into 107b fract. */
00475     LO <<=21;
00476 #if QUAD_DEBUG
00477   printf("LO&HI combined: HI=0x%016llx\n",HI);
00478   printf("                LO=0x%016llx\n",LO);
00479 #endif
00480 
00481   }
00482 
00483   /* At this point we have a 128b binary fraction HILO, the low 21
00484    * bits of which are zero, and a biased binary exponent exphi.
00485    *
00486    * Generate the decimal exponent by multiplying the binary 
00487    * exponent by log10(2). 
00488    */
00489 
00490   /* exphi is at most 9 bits plus sign. Shifts avoid overflow
00491    * and recover integer part.
00492    */
00493 
00494   dexp = ( ((INT32) (exphi-1023)) * (0x4D104D42L >> 21) );
00495   /* truncate to -infinity */
00496   dexp >>= 11;
00497   dexp++;
00498 
00499 #if QUAD_DEBUG
00500   printf("before scaling: dexp=%d\n",dexp);
00501 #endif
00502 
00503   /* Scale fraction by decimal exponent. */
00504 
00505   if(dexp)      /* value *= 10 ** -dexp */
00506     c_qtenscale(&value, -dexp, &bexp); 
00507   else
00508     bexp=0;
00509 
00510 #if QUAD_DEBUG
00511   printf("after scaling: HI=0x%016llx\n",HI);
00512   printf("               LO=0x%016llx\n",LO);
00513   printf("               bexp=%d\n",bexp);
00514 #endif
00515 
00516   exphi += bexp;
00517 
00518 #if QUAD_DEBUG
00519   printf("adjusted exphi: exphi=%d\n",exphi);
00520 #endif
00521 
00522   /* unbias exphi */
00523 
00524   exphi -= 1022;
00525 
00526   /* Remove binary exponent by denormalizing */
00527   /* Shift an extra 4 bits to allow for later digit generation */
00528 
00529   LO = (LO >> 4-exphi) | (HI << (60+exphi));
00530   HI >>= 4-exphi;
00531 
00532 #if QUAD_DEBUG
00533   printf("denormalized: HI=0x%016llx\n",HI);
00534   printf("              LO=0x%016llx\n",LO);
00535 #endif
00536 
00537   /* Split LO inpreparation to producing digits */
00538 
00539   lohi = LO >> 60;
00540   lolo = LO & ((1ULL<<60)-1);
00541 
00542 #if QUAD_DEBUG
00543   printf("split LO: lohi=0x%02x\n",lohi);
00544   printf("          lolo=0x%016llx\n",lolo);
00545 #endif
00546 
00547   /* find first non-zero digit */
00548   while( !(HI >> 60) ) {
00549     lolo *= 10;
00550     lohi = lohi * 10 + (lolo>>60);
00551     lolo &= (1ULL<<60)-1;
00552     HI = HI * 10 + (lohi>>4);
00553     lohi &= 0xF;
00554     dexp--;
00555 
00556 #if QUAD_DEBUG
00557   printf("finding first digit: HI=0x%016llx\n",HI);
00558   printf("                     lohi=0x%02x\n",lohi);
00559   printf("                     lolo=0x%016llx\n",lolo);
00560   printf("                     dexp=%d\n",dexp);
00561 #endif
00562 
00563   }
00564   /* produce digits */
00565 
00566   if(fflag) {     /* f format */
00567     nd = min(34, ndigit + dexp + 1); /* need 1 integer digit */
00568   }
00569   else {
00570     nd = ndigit;    /* e format */
00571   }
00572 
00573 #if QUAD_DEBUG
00574   printf("number of digits: nd=%d\n",nd);
00575   printf("                  ndigit=%d\n",ndigit);
00576   printf("                  dexp=%d\n",dexp);
00577 #endif
00578 
00579   /* produce digits - always at least one */
00580   do {
00581     *buffer++ = '0' + (HI>>60);
00582     HI &= (1ULL<<60) - 1;
00583     lolo *= 10;
00584     lohi = lohi * 10 + (lolo>>60);
00585     lolo &= (1ULL<<60)-1;
00586     HI = HI * 10 + (lohi>>4);
00587     lohi &= 0xF;
00588 
00589 #if QUAD_DEBUG
00590   printf("another digit: HI=0x%016llx\n",HI);
00591   printf("               lohi=0x%02x\n",lohi);
00592   printf("               lolo=0x%016llx\n",lolo);
00593   *buffer = '\0';
00594   printf("               buffer=%s\n",buffer0);
00595 #endif
00596 
00597   }
00598   while(--nd > 0);
00599 
00600   /* decimal round */
00601 
00602   /*  *buffer guard rest  Action
00603    *  
00604    *  dc  0-4 dc  none
00605    *  dc  6-9 dc  round
00606    *  odd 5 dc  round
00607    *  even  5 0 none
00608    *  even  5 !=0 round
00609    */
00610 
00611   guard = HI>>60;
00612   if(guard >= 5) {
00613     if(guard > 5 ||
00614        (*buffer-1)&1 ||   /* buffer odd */
00615        HI & ((1ULL<<60) - 1) ||
00616        lohi || lolo) {    /* rest non-zero */
00617       char *b=buffer;   /* buffer pointer */
00618 
00619       /* must round */
00620 
00621 #if QUAD_DEBUG
00622   printf("must round\n");
00623 #endif
00624       for(b=buffer-1; ; b--) {
00625   if(*b < '0' || '9' < *b) { /* attempt to carry into sign */
00626     *(b+1) = '1';
00627     dexp++;
00628     if(fflag && nd != min(34, ndigit + dexp + 1)) { 
00629       *buffer++ = '0';  /* need one more digit */
00630 #if QUAD_DEBUG
00631       printf("added a zero: fflag=%d\n",fflag);
00632       printf("              ndigit=%d\n",ndigit);
00633       printf("              nd=%d\n",nd);
00634 #endif
00635     }
00636 
00637 #if QUAD_DEBUG
00638   printf("carry into sign: dexp=%d\n",dexp);
00639   printf("                 buffer=%s\n",buffer0);
00640 #endif
00641 
00642     break;
00643   }
00644   (*b)++;
00645 
00646 #if QUAD_DEBUG
00647   *buffer = '\0';
00648   printf("round: buffer=%s\n",buffer0);
00649 #endif
00650 
00651   if(*b > '9') {    /* carry */
00652     *b = '0';
00653 
00654 #if QUAD_DEBUG
00655   printf("carry: buffer=%s\n",buffer0);
00656 #endif
00657 
00658   }
00659   else
00660     break;    /* no carry */
00661       }
00662       
00663     }
00664   }
00665 
00666 #if REALLYNEEDED
00667   if(fflag && ndigit+dexp < 0) { /* no significant digits */
00668 
00669     dexp--;     /* This is a lie. However, that's the
00670          * way dtoa behaves (sometimes!), so 
00671          * make qtoa work the same way.
00672          */
00673   }
00674 #endif
00675 
00676   *buffer = '\0';
00677 
00678 #if QUAD_DEBUG
00679   printf("done: buffer=%s\n",buffer0);
00680   printf("      dexp=%d\n", dexp);
00681 #endif
00682 
00683   return dexp;
00684 }
00685
osprey/common/util/c_q_to_a.c
