C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
realio.cpp
00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
00007 **                          Universitaet Wuppertal, Germany   
00008 **
00009 **  This library is free software; you can redistribute it and/or
00010 **  modify it under the terms of the GNU Library General Public
00011 **  License as published by the Free Software Foundation; either
00012 **  version 2 of the License, or (at your option) any later version.
00013 **
00014 **  This library is distributed in the hope that it will be useful,
00015 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 **  Library General Public License for more details.
00018 **
00019 **  You should have received a copy of the GNU Library General Public
00020 **  License along with this library; if not, write to the Free
00021 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 /* CVS $Id: realio.cpp,v 1.31 2014/01/30 17:23:48 cxsc Exp $ */
00025 
00026 #include <cstring>
00027 #include "real.hpp"
00028 #include "ioflags.hpp"
00029 #include "RtsFunc.h"
00030 #include "dot.hpp"
00031 
00032 namespace cxsc {
00033 
00034 int d_init_dm (void);
00035 void d_outp(char *buffer, Dotprecision c,
00036               int FormatFlag, int FracDigits, int rnd,
00037               int *length);
00038 
00039 #if _WIN32
00040 extern __declspec(thread) char *dm;
00041 #elif __APPLE__ && !CXSC_FORCE_TLS
00042 extern char *dm; 
00043 #else
00044 extern __thread char *dm;
00045 #endif
00046 
00047 
00048 char* GetHexDigits (char* s, a_btyp& x, int count);
00049 
00050 int r_outpx (char *buffer, a_real s,a_intg FormatFlag, a_intg FracDigits, a_intg rnd, a_intg *length);
00051 
00052 #define WASGOOD          0 
00053 #define ISINFINITYPLUS  -1
00054 #define ISINFINITYMINUS -2
00055 #define ISQUIETNAN      -3
00056 #define ISSIGNALINGNAN -4
00057 
00058 std::string realToHex(const real& a) 
00059 {
00060    a_btyp* b = (a_btyp*) &a;
00061    int i;
00062    char cs[256];
00063 
00064    for (i=0; i < realwidth-19; i++) cs[i] = ' ';
00065    cs[i] = 0;
00066    sprintf (&cs[strlen(cs)], "%s", ((b[HIGHREAL] & 0x80000000L) ? "-" : "+"));
00067    sprintf (&cs[strlen(cs)],"%c",   '1');
00068    sprintf (&cs[strlen(cs)],"%05lX", (b[HIGHREAL] & 0xFFFFFL));
00069    sprintf (&cs[strlen(cs)],"%08lX",(long unsigned int)(b[LOWREAL]));
00070    sprintf (&cs[strlen(cs)],"e%03X",(unsigned int)((b[HIGHREAL] >> 20) & 0x7FF));
00071    return string(cs);
00072 }
00073 
00074 std::string & operator<<(std::string &s, const real& a) throw()
00075 {
00076 #if _WIN32
00077    static __declspec(thread) char cs[256];
00078 #elif __APPLE__ && !CXSC_FORCE_TLS
00079     static char cs[256];
00080 #else
00081    static __thread char cs[256];
00082 #endif
00083    
00084    if (ioflags.isset(IOFlags::hex))
00085    {
00086       a_btyp* b = (a_btyp*) &a;
00087       int i;
00088       char cs[256];
00089 
00090       for (i=0; i < realwidth-19; i++) cs[i] = ' ';
00091       cs[i] = 0;
00092       sprintf (&cs[strlen(cs)], "%s", ((b[HIGHREAL] & 0x80000000L) ? "-" : "+"));
00093       sprintf (&cs[strlen(cs)],"%c",   '1');
00094       sprintf (&cs[strlen(cs)],"%05lX", (b[HIGHREAL] & 0xFFFFFL));
00095       sprintf (&cs[strlen(cs)],"%08lX",(long unsigned int)(b[LOWREAL]));
00096       sprintf (&cs[strlen(cs)],"e%03X",(unsigned int)((b[HIGHREAL] >> 20) & 0x7FF));
00097       s+=cs;
00098    } else if (ioflags.isset(IOFlags::rndnone))
00099    {
00100       if (IsSignalingNaN(a)) s+="<SignallingNaN>";
00101       else if (IsQuietNaN(a)) s+="<QuietNaN>";
00102       else if (IsInfinity(a)) s+="<Infinity>";
00103       else 
00104       {
00105          if (realdigits && realwidth)
00106             sprintf (cs, "%*.*g", realwidth, realdigits, a.w); // no need for "lg"
00107          else if (realwidth)
00108             sprintf (cs, "%*g", realwidth, a.w);
00109          else
00110             sprintf (cs, "%g", a.w);
00111          s+=cs;
00112       }
00113    } else
00114    {
00115       rndtype rnd;
00116       a_intg length, formatflag, addblanks;
00117       a_intg digits = realdigits;
00118       char *str;
00119 
00120       if (d_init_dm () == -1) 
00121       { 
00122          // throw
00123          // errmon (ERR_ALL(NOMOREMEMORY));
00124          // errmon (ERR_ALL(NOCONTINUEPOSSIBLE));
00125       }
00126 
00127       if (ioflags.isset(IOFlags::rndup)) 
00128          rnd = RND_UP;
00129       else if (ioflags.isset(IOFlags::rnddown)) 
00130          rnd = RND_DOWN;
00131       else 
00132          rnd = RND_NEXT;
00133 
00134       if (ioflags.isset(IOFlags::variable)) 
00135          formatflag = realwidth;
00136       else if (ioflags.isset(IOFlags::varfixwidth)) 
00137          formatflag = realwidth, digits = -digits;
00138       else 
00139          formatflag = (ioflags.isset(IOFlags::fixed)) ? 0 : -1;
00140 
00141       switch (r_outpx (dm, a.w, formatflag, digits, rnd, &length)) 
00142       { 
00143          case WASGOOD:
00144             dm[length] = 0;
00145             str = dm;
00146             if (*str == '+') 
00147             {
00148                if (ioflags.isset(IOFlags::blank))
00149                   *str = ' ';
00150                else if (ioflags.isset(IOFlags::noblank)) 
00151                   str++;
00152             }
00153             break;
00154          case ISINFINITYPLUS:  
00155             str = (char*)"<+Infinity>";     
00156             break;
00157          case ISINFINITYMINUS: 
00158             str = (char*)"<-Infinity>";     
00159             break;
00160          case ISQUIETNAN:      
00161             str = (char*)"<QuietNaN>";      
00162             break;
00163          case ISSIGNALINGNAN:  
00164             str = (char*)"<SignalingNaN>";  
00165             break;
00166          default:              
00167             str = (char*)"<ERROR>";
00168             break;
00169       }
00170       length = strlen(str);
00171       addblanks = (length < realwidth) ? realwidth - length : 0;
00172 
00173       if (ioflags.isset(IOFlags::rightjust)) 
00174          for (;addblanks; addblanks--) 
00175             s+= ' ';
00176 
00177       s+=str;
00178       
00179       for (;addblanks; addblanks--) 
00180          s+= ' ';
00181    }
00182    return s;
00183 }
00184 
00185 std::ostream & operator <<(std::ostream &o,const real &a) throw()
00186 {
00187    std::string s="";
00188    s << a;
00189    o << s;
00190    return o;
00191 }
00192 
00193 std::string & operator>> (std::string & str, real& a) throw()
00194 {
00195    char *s=new char[str.size()+1];
00196    char *orgs=s;
00197    strcpy(s,str.c_str());  
00198 
00199    if (ioflags.isset(IOFlags::hex))
00200    {
00201       a_btyp* b = (a_btyp*) &a;
00202       a_btyp x;
00203 
00204       b[0] = b[1] = 0;
00205       s = cskipwhitespaces (s);
00206       if (*s == '-') 
00207       {
00208          b[HIGHREAL] |= 0x80000000L; 
00209          s++; 
00210       } else if (*s == '+') 
00211          s++;
00212 
00213     // ----------------------------------------------------
00214     // es wird ohne Pruefung folgendes Format vorausgestzt :
00215     //    1xxxxxxxxxxxxxeXXX
00216     // wobei x = HexDigits der Mantisse
00217     //       X = HexDigits des Exponents
00218                       // saemmtliche HexDigits muessen gross geschrieben sein !
00219 
00220       if (*s) 
00221          s++;                                      // skip '1'
00222       s = GetHexDigits (s, x, 5);  
00223       b[HIGHREAL] |= x;    // get mantissa
00224       s = GetHexDigits (s, x, 8);  
00225       b[LOWREAL]   = x;
00226       if (*s) 
00227          s++;                                      // skip 'e'     
00228       s = GetHexDigits (s, x, 3);
00229       b[HIGHREAL] |= x << 20;                 // get exponent
00230       if (*s) 
00231          s++;                            // skip at least one more char
00232    } else
00233    {
00234       rndtype rndfl;
00235 
00236       if (ioflags.isset(IOFlags::rndup)) 
00237          rndfl = RND_UP;
00238       else if (ioflags.isset(IOFlags::rnddown)) 
00239          rndfl = RND_DOWN;
00240       else 
00241          rndfl = RND_NEXT;
00242 
00243       str=s;
00244       dotprecision dot;
00245       str >> dot;
00246       strcpy(s,str.c_str()); // Ooooooooohhhh... :((
00247       a = rnd (dot, rndfl);
00248    }
00249    str=s;
00250    delete [] orgs;
00251    return str;
00252 }   
00253 void operator >>(const char *a,real &b) throw()
00254 {
00255    std::string c(a);
00256    c>>b;
00257 }
00258 void operator >>(const string &a,real &b) throw()
00259 {
00260    std::string c(a);
00261    c>>b;
00262 }
00263 std::istream& operator>> (std::istream& s, real& a) throw()
00264 {
00265    if (ioflags.isset(IOFlags::hex))
00266    {
00267       char inp[20];
00268       int i;
00269       char c;
00270 
00271       // !! There are no checks about the right input-format, it is assumed
00272       //    the fixed format : S1xxxxxxxxxxxxxeXXX
00273       //    whereby S   means the sign
00274       //            '1' is the leading implicit binary one of the mantissa
00275       //            x   are the mantissa hexdigits
00276       //            'e' signals the beginning of the exponent
00277       //            X   are the exponent hexdigits
00278 
00279       // skip white spaces and
00280       // get the sign, the 14 Mantissadigits and the exponent
00281       // (total 19 chars)
00282                                                
00283       c = skipwhitespaces (s);
00284       for (i=0; i < 19; i++) 
00285       {
00286          inp[i] = c;
00287          if (s.good()) s.get(c); else c = '\0';
00288       }
00289       inp[i] = '\0';
00290 
00291       inp >> a;
00292    } else
00293    {
00294       rndtype rndfl;
00295 
00296       if (ioflags.isset(IOFlags::rndup)) 
00297          rndfl = RND_UP;
00298       else if (ioflags.isset(IOFlags::rnddown)) 
00299          rndfl = RND_DOWN;
00300       else 
00301          rndfl = RND_NEXT;
00302 
00303       dotprecision dot;
00304       s >> dot;
00305       a = rnd (dot, rndfl);
00306    }
00307    return s;
00308 }                                                                                                                    
00309 
00310 //----------------------------------------------------------------------------
00311 // GetHexDigits
00312 //
00313 //  Interpretiert die naechsten "count" Zeichen aus dem String "s"
00314 //  als HexDigits, der binaere Wert wird dann in "x" zurueckgegeben.
00315 //  Als Returnwert wird ein Zeiger auf die Stringposition nach den
00316 //  HexDigits zurueckgegeben.
00317 
00318 char* GetHexDigits (char* s, a_btyp& x, int count)
00319 {
00320   int i, c;
00321 
00322   for (x=0,i=0; i < count && *s; s++,i++) {
00323     if ((c = *s) >= 'A') c -= 'A' - 10; else c -= '0';
00324     if (c < 0 || c > 0xF) c = 0;
00325     x = (x << 4) | c;
00326   }
00327   return s;
00328 }        
00329 
00330 } // namespace cxsc
00331 
00332 /****************************************************************/
00333 /*                                                              */
00334 /*      Filename        : r_outpx.c                             */
00335 /*                                                              */
00336 /*      Entries         : void r_outp                           */
00337 /*                         (buffer,s,FormatFlag,                */
00338 /*                          FracDigits,rnd,length)              */
00339 /*                        char *buffer;                         */
00340 /*                        a_real s;                             */
00341 /*                        a_intg rnd,FormatFlag,FracDigits;     */
00342 /*                        a_intg *length;                       */
00343 /*                                                              */
00344 /*      Arguments       : buffer - output buffer holding string */
00345 /*                        s - IEEE value                        */
00346 /*                        FormatFlag - format selection         */
00347 /*                              -1 = scientific format          */
00348 /*                               0 = fixed format               */
00349 /*                             > 0 = variable format            */
00350 /*                                   value is the total field   */
00351 /*                                   width                      */
00352 /*                        FracDigits - number of fraction digits*/
00353 /*                        rnd - rounding mode                   */
00354 /*                              -1 = round downwards            */
00355 /*                               0 = round to nearest           */
00356 /*                               1 = round upwards              */
00357 /*                        length - size of buffer string        */
00358 /*                                                              */
00359 /*      Description     : Convert an IEEE double format number  */
00360 /*                        to a character string.                */
00361 /*                                                              */
00362 /****************************************************************/
00363 
00364 /*#ifndef ALL_IN_ONE
00365 #ifdef AIX
00366 #include "/u/p88c/runtime/o_defs.h"
00367 #else
00368 #include "o_defs.h"
00369 #endif
00370 #define local
00371 #endif*/
00372 
00373 namespace cxsc {
00374 
00375 #define WASGOOD          0 
00376 #define ISINFINITYPLUS  -1
00377 #define ISINFINITYMINUS -2
00378 #define ISQUIETNAN      -3
00379 #define ISSIGNALINGNAN -4
00380 
00381 #define MANT_INFINITY(a)        ((a)[0]==HIDDEN_BIT && (a)[1]==ZERO)
00382 #define SIGNALING(a)            ((a) & SIGNAL_BIT)
00383 #define SIGNAL_BIT              ((a_btyp)0x00080000L)  
00384 
00385 int r_outpx(char *buffer, a_real s, a_intg FormatFlag,
00386             a_intg FracDigits, a_intg rnd, a_intg *length)
00387 {
00388    a_intg ActWidth,DecPlaces,expo,IntDigits,MinNumChars;
00389    a_intg dexpo,digits,bdp,k,l,addpoint;
00390    a_intg HoldWidth = (FracDigits < 0);
00391    a_bool vz,zero;
00392    a_btyp mant[BSIZE];
00393 
00394    if (HoldWidth) 
00395       FracDigits = -FracDigits;
00396    *length = 0;
00397 
00398    zero = b_deko(s,&expo,mant,&vz);
00399 
00400 /*                                                                   */
00401 /* infinity or NaN                                                   */
00402 /*                                                                   */
00403    if (expo>EXPO_MAX) 
00404    {
00405       if (MANT_INFINITY(mant)) 
00406       {                         /* Infinity   */
00407          k = (vz) ? ISINFINITYMINUS : ISINFINITYPLUS;
00408       } else 
00409       {                                             /* NaN        */
00410          k = (SIGNALING(mant[0])) ? ISSIGNALINGNAN : ISQUIETNAN;
00411       }
00412       return k;
00413    }
00414 
00415    for (k=D_U_RATIO;k<BSIZE;k++) 
00416       mant[k] = ZERO;
00417 
00418    if (vz) 
00419       rnd = -rnd;
00420 
00421    if (FormatFlag > 0 && FormatFlag-FracDigits <= 2) 
00422    {
00423       FormatFlag = FracDigits+3;
00424    }
00425    if (!zero && FormatFlag > 0) 
00426    {
00427       dexpo = (a_intg) ((expo*30103L)/100000L);
00428       if (dexpo < -((FracDigits+1)/2) ||
00429          (dexpo > 0 && dexpo >= FormatFlag-FracDigits-2)) 
00430       {
00431          FormatFlag = -1;
00432          if (HoldWidth) 
00433          {
00434             FracDigits = (FracDigits < 5) ? 0 : FracDigits-5;
00435          }
00436       }
00437    }
00438 
00439    if (expo>800) 
00440       bdp = (BUFFERSIZE-8)-(FormatFlag == -1 ? 0 : FracDigits);
00441    else if (expo<-800) 
00442       bdp = 8;
00443    else 
00444       bdp = B_D_P;
00445 
00446    addpoint = (FracDigits == 0) ? 0 : 1;
00447 /*                                                                   */
00448 /* floating-point representation                                     */
00449 /*                                                                   */
00450    if (FormatFlag == -1) 
00451    {
00452       DecPlaces = FracDigits;
00453       ActWidth =  DecPlaces+ExpDigits+4+addpoint;
00454 /*                                                                   */
00455 /* number is zero                                                    */
00456 /*                                                                   */
00457       if (zero) 
00458       {
00459          *buffer++ = (vz) ? '-' : '+';
00460          *buffer++ = '0';
00461          if (FracDigits) 
00462          {
00463             *buffer++ = XSC_DECIMAL_POINT;
00464             for (k=0;k<DecPlaces;k++) *buffer++ = '0';
00465          }
00466          *buffer++ = EXPONENT_E;
00467          *buffer++ = PLUS_SIGN;
00468          for (k=0;k<ExpDigits;k++) 
00469             *buffer++ = '0';
00470 
00471          *length = ActWidth;
00472 
00473          return WASGOOD;
00474       }
00475 /*                                                                   */
00476 /* determine output string                                           */
00477 /*                                                                   */
00478       digits = DecPlaces+(1+2);
00479       dexpo = -1;
00480       b_out(mant,expo,digits,buffer,&bdp,&dexpo);
00481 
00482       if (dexpo>0 && dexpo>DecPlaces+2) 
00483          digits = dexpo+1;
00484 
00485       b_rnd(rnd,buffer,digits,DecPlaces+1,&bdp,&dexpo);
00486 
00487       if (FracDigits) 
00488       {
00489          buffer[bdp-dexpo-1] = buffer[bdp-dexpo];
00490          buffer[bdp-dexpo] = XSC_DECIMAL_POINT;
00491       }
00492       buffer[bdp-dexpo+DecPlaces+1] = EXPONENT_E;
00493       buffer[bdp-dexpo+DecPlaces+2] = (dexpo<0) ? MINUS_SIGN : PLUS_SIGN;
00494       expo = (dexpo<0) ? -dexpo : dexpo;
00495       for (k=ExpDigits;k>0;k--) 
00496       {
00497          buffer[bdp-dexpo+DecPlaces+2+k] = expo%10+'0';
00498          expo /= 10;
00499       }
00500       *length = 2+addpoint+DecPlaces+2+ExpDigits;
00501 
00502       l = bdp - 1 - addpoint - dexpo;
00503       buffer[l] = (vz) ? '-' : '+';
00504       for (k=0;k<*length;k++,l++) 
00505          buffer[k] = buffer[l];
00506    } else
00507 /*                                                                   */
00508 /* fixed-point representation                                        */
00509 /*                                                                   */
00510    {
00511 /*                                                                   */
00512 /* number is zero                                                    */
00513 /*                                                                   */
00514       if (zero) 
00515       {
00516          *length = FracDigits+2+addpoint;
00517          *buffer++ = (vz) ? '-' : '+';
00518          *buffer++ = '0';
00519          if (FracDigits) 
00520          {
00521             *buffer++ = XSC_DECIMAL_POINT;
00522             for (k=0;k<FracDigits;k++) 
00523                *buffer++ = '0';
00524          }
00525          return WASGOOD;
00526       }
00527 
00528       /* estimate number of decimal digits */
00529       if (expo>=0) 
00530       {
00531          IntDigits = ((expo+1)*61)/200+1;
00532       } else 
00533       {
00534          IntDigits = 0;
00535          dexpo = 0;
00536       }
00537 
00538       /* fill fractional part with zeros                   */
00539       for (k=0;k<=FracDigits+2;k++) 
00540          buffer[bdp+k] = '0';
00541 
00542       digits = IntDigits+FracDigits+2;
00543       b_out(mant,expo,digits,buffer,&bdp,&dexpo);
00544 
00545       /* correct setting of IntDigits */
00546       if (expo>=0) 
00547       {
00548          IntDigits = dexpo+1;
00549       } else 
00550       {
00551          IntDigits = 1;
00552          digits++;
00553       }
00554 
00555       b_rnd(rnd,buffer,digits,IntDigits+FracDigits,&bdp,&dexpo);
00556 
00557       /* correct setting of IntDigits after rounding */
00558       IntDigits = dexpo+1;
00559 
00560       /* value is zero after rounding;                     */
00561       if (vz) 
00562       {
00563          for (k=(bdp+1)-IntDigits;k<(bdp+1)+FracDigits;k++) 
00564          {
00565             if (buffer[k]!='0') 
00566                break;
00567          }
00568          if (k==(bdp+1)+FracDigits) 
00569             vz = FALSE;
00570       }
00571 
00572       MinNumChars = IntDigits+FracDigits+1+addpoint;
00573 
00574       if (FracDigits) 
00575       {
00576          /* generate decimal point                            */
00577          for (k=bdp-IntDigits;k<bdp;k++) 
00578             buffer[k] = buffer[k+1];
00579          buffer[bdp] = XSC_DECIMAL_POINT;
00580       }
00581 
00582       *length = MinNumChars;
00583 /*                                                                   */
00584 /* sign of number                                                    */
00585 /*                                                                   */
00586       l = bdp - 1 - addpoint - dexpo;
00587       buffer[l] = (vz) ? '-' : '+';
00588       for (k=0;k<*length;k++,l++) 
00589          buffer[k] = buffer[l];
00590    }
00591 
00592    return WASGOOD;
00593 }
00594 
00595 } // namespace cxsc
00596