nsz Git - libfirm/blob - ir/tv/IeeeCC754/src/IeeeCC754.cc

   1  /*
   2 ##########################################################################
   3 #                                                                        #
   4 # Program: IeeeCC754                                                     #
   5 #                                                                        #
   6 # Description:                                                           #
   7 #   IeeeCC754 or IEEE 754 Compliance Checker is a precision and range    #
   8 #   independent tool to test whether an implementation of                #
   9 #   floating-point arithmetic (in hardware or software) is compliant     #
  10 #   with the principles of the IEEE 754-854 floating-point standards.    #
  11 #   You can find out more about the testing tool IeeeCC754 at            #
  12 #                                                                        #
  13 #         http://win-www.uia.ac.be/u/cant/ieeecc754.html                 #
  14 #                                                                        #
  15 #   This tool is in parts based on and greatly benefited from the        #
  16 #   the program FPTEST developed by Jerome Coonen. For a full            #
  17 #   description of the extensions to FPTEST and a reference to           #
  18 #   the original Coonen program, please refer to the URL given above.    #
  19 #   For the options available with the program IeeeCC754 and its         #
  20 #   compatibility with David Hough's hexadecimal UCB format, we          #
  21 #   also refer to the file readme.usage.                                 #
  22 #                                                                        #
  23 #  Usage: see readme.usage                                               #
  24 #                                                                        #
  25 #  Responsible authors:                                                  #
  26 #         Brigitte Verdonk                                               #
  27 #         Annie Cuyt                                                     #
  28 #                                                                        #
  29 #  Contributors:                                                         #
  30 #         Tarun Agarwal(05-07/2002)                                      #
  31 #         Johan Bogo (1998-1999)                                         #
  32 #         Tim Gevers (10-12/2000)                                        #
  33 #         Debby Ooms (1996-1997)                                         #
  34 #         Geert Vermuyten (1996-1997)                                    #
  35 #         Dennis Verschaeren (09/1996-06/2000)                           #
  36 #                                                                        #
  37 #  Copyright (C) 2000  University of Antwerp                             #
  38 #                                                                        #
  39 #  This program can be obtained from the authors, free, but WITHOUT ANY  #
  40 #  WARRANTY; without even the implied warranty of MERCHANTABILITY or     #
  41 #  FITNESS FOR A PARTICULAR PURPOSE.                                     #
  42 #                                                                        #
  43 #  Contact:                                                              #
  44 #       Brigitte.Verdonk@uia.ua.ac.be                                    #
  45 #       Department of Mathematics and Computer Science                   #
  46 #       University of Antwerp (UIA)                                      #
  47 #       Universiteitsplein 1                                             #
  48 #       B2610 Antwerp, BELGIUM                                           #
  49 #                                                                        #
  50 ##########################################################################
  51
  52 Filename:
  53        $RCSfile$
  54
  55 Last updated:
  56        $Date$
  57
  58 */
  59
  60 #include <stdlib.h>
  61 #include <stdio.h>
  62 #include <ctype.h>
  63 #include <string.h>
  64 #include <iostream.h>
  65 #include <fstream.h>
  66 #include <strstream.h>
  67 #include <iomanip.h>
  68 #include <math.h>
  69 #include <Bitstring.h>
  70 #include <UCB.h>
  71 #include <limits.h>
  72 #include <string>
  73 extern "C" {
  74 typedef double LLDBL;
  75 #include <fltcalc.h>
  76 }
  77
  78 using namespace std;
  79
  80 #include <DriverFloatRepr.h>
  81 template<class M> class UCB;
  82 UCB<DriverFloatRepr> ucb;
  83
  84 #define TD_OVERFLOW             0x08
  85 #define TD_UNDERFLOW            0x10
  86 #define TD_UNDERFLOWV           0x40  // inexactness, no denormalization loss
  87 #define TD_UNDERFLOWW           0x80  // tininess before rounding
  88
  89 #define TD_INVALID              0x01
  90 #define TD_INEXACT              0x20
  91 #define TD_DIVBYZERO            0x04
  92 #define IEEEDEFAULTENV          0
  93
  94 #define TD_TONEAREST    0  // deze numberlen worden gebruikt       // FP_RN
  95 #define TD_UPWARD       2  // ipv de originele roundingnumberlen   // FP_RP
  96 #define TD_DOWNWARD     3  // zoals die voorkomen in "ieeefp.h"   // FP_RM
  97 #define TD_TOWARDZERO   1  // voor compatibiliteit met de UCB     // FP_RZ
  98 // outputfile qua volgorde van
  99 // rounding
 100 #define ALL_FORMATS 3
 101 #define TESTEVEN 1
 102 #define TESTODD 2
 103 #define NO_FLAGS 31
 104 #define NO_FLAGS_INEXACT      16
 105 #define NO_FLAGS_OVERFLOW     8
 106 #define NO_FLAGS_UNDERFLOW    4
 107 #define NO_FLAGS_INVALID      2
 108 #define NO_FLAGS_DIV_BY_ZERO  1
 109 #define ROUND_ALL  31
 110 #define ROUND_NEAREST  8
 111 #define ROUND_UP       4
 112 #define ROUND_DOWN     2
 113 #define ROUND_ZERO     1
 114 #define ALL_MODES ((1<<TD_TONEAREST) | (1<<TD_DOWNWARD) | (1<<TD_UPWARD) | (1<<TD_TOWARDZERO))
 115 #define NO_OF_BITS   32
 116 #define USAGE "Error in calling " << argv[0] <<endl<<\
 117 "For instructions on calling the program, see the readme.usage file."<<endl\
 118
 119
 120 int cant_infinity=0 , small_norm=0;
 121 int HexNoIsGiven=0;
 122 int posDec;
 123 int posUpDec;
 124 int incexp=0;
 125 int NaN1,NaN2,NaN3;
 126 int err1,err2,err3;
 127 int isIeeeVector=0;
 128 unsigned long expon;
 129 unsigned long hidden;
 130 unsigned long mant;
 131 unsigned long length;
 132 unsigned long rest;
 133 unsigned long sign_exp_length;
 134 unsigned long kl_biased_exp;
 135 unsigned long sign_exp_rest ;
 136 FILE *input=NULL, *output=NULL;
 137 unsigned long readings;
 138 unsigned long  sexpon;
 139 unsigned long dexpon;
 140 unsigned long  smant=0;
 141 unsigned long dmant=0;
 142 unsigned long  shidden;
 143 unsigned long dhidden;
 144 unsigned long  slength;
 145 unsigned long  srest;
 146 unsigned long dlength;
 147 unsigned long drest;
 148 unsigned long ssign_exp_length; // sign+exponent length(in multiples of NO_OF_BITS)
 149 unsigned long ssign_exp_rest;
 150 unsigned long dsign_exp_length;
 151 unsigned long dsign_exp_rest;
 152 unsigned long skl_biased_exp;
 153 unsigned long dkl_biased_exp;
 154 unsigned long testFormats;
 155 unsigned long testModes;
 156 unsigned long flushVector;
 157 unsigned long killerVector;
 158 unsigned long *operand1, *operand2, *result;
 159 long xcp, dots, wrong_input, lines_in, err_form;
 160 char dot, opcode, op2code, op3code, op4code;
 161 int comp_sort;
 162 int CmndLinePos=1, def=0, defe = 0,deft=0,falt =0;
 163 int ucbInput=0;
 164 int check =-1;
 165 int noFlags=0;
 166 int Round =ROUND_ALL;
 167 int dectobin=0;
 168 int tiny=1;
 169 int inf=1;
 170 int NaN=1;
 171 int signedZero=1;
 172 int ieee=1;
 173 char *intstr = NULL;
 174 int skipVector=0,jumpOverflow=0,jumpUnderflow=0,jumpInvalid=0,jumpDivZero=0;
 175 int pre = 0, post = 0; // to accomodate '?' before or after a result
 176 int NaNresult = 0;
 177 char *binary;
 178 char *upbinary;
 179 int nextdig = 0;
 180 int mysignbit = 0;
 181 char *decimal = new char[ maxstr ];
 182 char *updecimal = new char[ maxstr ];
 183 char *logfilename = "ieee.log";
 184
 185
 186 void ParseArgs( int, char** );
 187 void ParsecaseO( int argc, char** argv);
 188 void ParsecaseS();
 189 void ParsecaseD( int argc, char** argv);
 190 void ParsecaseL();
 191 void ParsecaseQ();
 192 void ParsecaseM();
 193 void ParsecaseE( int argc, char** argv);
 194 void ParsecaseT( int argc, char** argv);
 195 void ParsecaseR( int argc, char** argv);
 196 void ParsecaseN( int argc, char** argv);
 197 void ParsecaseJ( int argc, char** argv);
 198 void ParsecaseH();
 199 void PostParseAction( int argc, char** argv );
 200 void readParam( unsigned long & );
 201 int readAline( );
 202 void Calc_General_Stuff( );
 203 void putDot( );
 204 void getModes( );
 205 void getExceptions( );
 206 unsigned long getinteger();
 207 unsigned long getHex();
 208 void getNumber( unsigned long *number);
 209 void getbinary( unsigned long *number);
 210 void getdecimal( );
 211 void getFPsig(unsigned long *number);
 212 void getFPexp(unsigned long *number);
 213 void getBbias(unsigned long *number);
 214 void getBsignificant(unsigned long *number);
 215 long int getBexpon();
 216 void getDsignificand();
 217 void getDexpon();
 218 void exitmem( );
 219 int FindNextSpace(  string pstring , int beginpos, int occurance );
 220 void writeUCB();
 221 void writeUCBopcode(ostrstream &stream);
 222 void writeUCBprecision(ostrstream &stream);
 223 void writeUCBround(int m,ostrstream &stream);
 224 void writeUCBexcep(ostrstream &stream);
 225 void writeUCBoperand(ostrstream &stream);
 226 void writeUCBresult(int m,ostrstream &stream);
 227 void writeUCBcomutative(ostrstream &stream);
 228 void writeUCBresultBin (int m, ostrstream &stream);
 229 void writeUCBresultDec (int m, ostrstream &stream);
 230 void writeUCBresultNum (ostrstream &stream);
 231 void getFPInteger( unsigned long *number);
 232 void Infinity( unsigned long *number);
 233 void Quiet_NaN( unsigned long *number);
 234 void Signaling_NaN( unsigned long *number);
 235 void Smallest_Norm( unsigned long *number);
 236 void ModifierPorM(unsigned long *number,int McasePorM);
 237 void getModifier(unsigned long *number);
 238 void ModifierIorD(unsigned long *number,int McaseIorD);
 239 void getFPRootNumber(  unsigned long *number);
 240 void getFPHexNo();
 241 int getSign();
 242 void initialize_values(int source);
 243 unsigned long null_exp( const unsigned long *number);
 244 unsigned long one_exp( const unsigned long *number);
 245 void Plus( unsigned long *number, unsigned long k);
 246 void PlusBias( unsigned long *number);
 247 void Minus( unsigned long *number, unsigned long k);
 248 void MinusBias( unsigned long *number);
 249 void PlusBk( unsigned long *number, unsigned long k );
 250 void MinusBk( unsigned long *number, unsigned long k);
 251 void Incr_at_pos( unsigned long *number, unsigned long pos, unsigned long k);
 252 void Incr( unsigned long *number, unsigned long k);
 253 void Decr_at_pos (long unsigned int *, long unsigned int, long unsigned int);
 254 void Decr( unsigned long *number, unsigned long k);
 255 void ModifierU( unsigned long *number, unsigned long ulps);
 256 char* WriteNumber( const unsigned long *, unsigned long, unsigned long,unsigned long, unsigned long, const int );
 257 void skipSpace( );
 258 int mygetc( FILE* );
 259 void myungetc( int, FILE* );
 260
 261
 262
 263 /*-----------------------------------------------------------------------------
 264 name        :main()
 265 description :this is the Ieee compliance checker for floating point
 266                 implementation according to IEEE754-854 specifications.
 267 called from :/
 268 call fns    :ParseArgs(),Calc_General_Stuff(),readAline(),putdot(),ucb.close()
 269 exceptions  :/
 270 algorithm   :
 271 Global var. used :
 272 -----------------------------------------------------------------------------*/
 273
 274
 275 int main( int argc,char * argv[ ] )
 276 {
 277   init_fltcalc(0);
 278   readings=dots=wrong_input=lines_in=err_form=0;
 279   if ( argc >1 )
 280
 281     ParseArgs( argc, argv ); /* to parse arguments in the executing statement*/
 282
 283   else
 284
 285     {
 286       cout << USAGE;
 287       exit( 1 );
 288     }
 289
 290   Calc_General_Stuff( );
 291
 292   operand1 = new unsigned long[ slength ];
 293   if ( operand1 == NULL )
 294     {
 295       printf( "\nNot enough memory!\n" );
 296       printf( "Program aborted.\n" );
 297       fclose( input );
 298       fclose( output );
 299       exit( EXIT_SUCCESS );
 300     }
 301
 302   operand2 = new unsigned long[ slength ];
 303   if ( operand2 == NULL )
 304     {
 305       printf( "\nNot enough memory\n" );
 306       printf( "Program aborted.\n" );
 307       fclose( input );
 308       fclose( output );
 309       delete[ ]  operand1;
 310       exit( EXIT_SUCCESS );
 311     }
 312
 313   result = new unsigned long[ dlength ];
 314
 315   if ( result == NULL )
 316     {
 317       printf( "\nNot enough memory\n" );
 318       printf( "Program aborted.\n" );
 319       fclose( input );
 320       fclose( output );
 321       delete[ ]  operand1;
 322       delete[ ]  operand2;
 323       exit( EXIT_SUCCESS );
 324     }
 325
 326   dot='.';
 327
 328   while ( readAline( ) )
 329     {
 330       skipVector = 0;
 331       readings += 1;
 332       putDot( );
 333     }
 334
 335   printf( "\n\nCounting %ld lines.\n", readings );
 336
 337   delete[ ]  operand1;
 338   delete[ ]  operand2;
 339   delete[ ]  result;
 340   ucb.Close( tiny );
 341   fclose( input );
 342   if ( output )
 343     fclose( output );
 344   return 0;
 345 }
 346
 347
 348
 349
 350 /*-----------------------------------------------------------------------------
 351 name        :ParseArgs()
 352 description :parse the command line arguments to get the precision and
 353                 rounding modes.
 354 called from :main()
 355 call fns    :ParsecaseO/S/D/L/Q/M/E/T/R/N/J/H , PostParseAction
 356 Global var. used :check, def,defe,deft,flat,ucbInput
 357 -----------------------------------------------------------------------------*/
 358 void ParseArgs( int argc, char** argv )
 359 {
 360   int j;
 361   dmant = 0;
 362   dexpon = 0;
 363   shidden=0;
 364
 365   while ( CmndLinePos < argc )
 366     {
 367       if ( ( strlen( argv[ CmndLinePos ] )  <= 3 )  && ( argv[CmndLinePos][ 0 ] == '-' ) )
 368         {
 369           switch ( argv[CmndLinePos][ 1 ] )
 370             {
 371             case 'c':                   //coonen format
 372               if ( check == -1 )
 373                 check = 1;
 374               else
 375                 falt =1;
 376               break;
 377
 378             case 'u':                   //ucb format (precisely one of
 379                                         //      -c | -u |-o should be there)
 380               if ( check !=-1 )
 381                 falt =1;
 382               def =1;
 383               ucbInput =1;
 384               break;
 385
 386             case 'o':                   //change coonen to ucb format
 387               ParsecaseO(argc,argv);
 388               break;
 389
 390
 391             case 's':                   //single precision
 392               ParsecaseS();
 393               break;
 394             case 'd':                   //if -d then double presion
 395               ParsecaseD(argc,argv);            // else specifications for destination
 396               break;
 397             case 'l':                   //long double precision
 398               ParsecaseL();
 399               break;
 400             case 'q':                   //quadruple precision
 401               ParsecaseQ();
 402               break;
 403             case 'm':                   //multi precision
 404               ParsecaseM();
 405               break;
 406             case 'e':                   //specify no. of bits for exponent
 407               ParsecaseE(argc,argv);
 408               break;
 409             case 't':                   //specify <int> bits precision
 410               ParsecaseT(argc,argv);
 411               break;
 412
 413             case 'r':                   //specify rounding modes
 414               ParsecaseR(argc,argv);
 415            break;
 416
 417             case 'i':                   //test conversions within range specified by ieee
 418               ieee=0;
 419               break;
 420
 421             case 'n':                   //do not test specified exceptions, infinity,...
 422               ParsecaseN(argc,argv);
 423               break;
 424
 425             case 'j':                   //jump/skip test vectors raising overflow, ...
 426                 ParsecaseJ(argc,argv);
 427               break;
 428
 429             case 'h':                   //specify if hidden bit
 430               ParsecaseH();
 431               break;
 432
 433             case 'f':                    // name of log file
 434               CmndLinePos++;
 435               logfilename = argv[CmndLinePos];
 436               break;
 437
 438             default:
 439               falt = 1;
 440             } // end switch
 441         }
 442       else
 443         {
 444           if ( CmndLinePos == argc-1 )          //check if test vector file exists
 445             {
 446               input = fopen( argv[CmndLinePos], "r" );
 447               if ( input == NULL )
 448                 {
 449                   printf( "ERROR: can't open %s.\n",argv[CmndLinePos] );
 450                   falt =1;
 451                 } else
 452                   printf( "Taking input from %s.\n",argv[CmndLinePos] );
 453             } else
 454               falt = 1;
 455         }
 456
 457       if ( falt )
 458         {
 459           cout << USAGE;
 460           exit( 1 );
 461         }
 462
 463     CmndLinePos++;
 464     } // end while
 465
 466
 467   if ( ( dmant == 0 )  || ( dexpon == 0 ) )   //if destination precision still has initial
 468     {                                         //values then equate them to that of source.
 469       dmant = smant;                          //(could be changed only for conversions)
 470       dexpon = sexpon;
 471       dhidden = shidden;
 472     } // if
 473
 474
 475   PostParseAction(argc,argv);
 476
 477 }
 478
 479
 480
 481
 482 /*-----------------------------------------------------------------------------
 483 name        :ParsecaseO()
 484 description :change coonen to ucb format
 485 called from :ParseArgs()
 486 call fns    :\
 487 Global var. used :check,CmndLinePos,falt
 488 -----------------------------------------------------------------------------*/
 489 void ParsecaseO( int argc, char** argv)
 490 {
 491   if ( check == -1 )
 492     {
 493       if ( ++CmndLinePos != argc )
 494         {
 495           if ( argv[CmndLinePos][ 0 ]  != '-' )
 496             {
 497               output = fopen( argv[CmndLinePos], "w" );
 498               if ( output == NULL )
 499                 {
 500                   printf( "ERROR: can't create %s.\n",argv[CmndLinePos] );
 501                   falt =1;
 502                 } else
 503                   printf( "Output sent to %s.\n",argv[CmndLinePos] );
 504             } else
 505               falt = 1;
 506         } else
 507           falt = 1;
 508       check = 0;
 509     } else
 510       falt = 1;
 511 }
 512
 513
 514 /*-----------------------------------------------------------------------------
 515 name        :ParsecaseS()
 516 description :specify single precision
 517 called from :ParseArgs()
 518 call fns    :\
 519 Global var. used :sexpon,smant,shidden,def,falt
 520 -----------------------------------------------------------------------------*/
 521 void ParsecaseS()
 522 {
 523   if ( !def )
 524     {
 525       sexpon = 8;
 526       smant = 23;
 527       shidden = 1;
 528       def = 1;
 529     } else
 530       falt =1;
 531 }
 532
 533
 534 /*-----------------------------------------------------------------------------
 535 name        :ParsecaseD()
 536 description :default is for double precision else it gets precision
 537         specifications for the destination (only used in conversions)
 538 called from :ParseArgs()
 539 call fns    :\
 540 Global var. used :CmndLinePos,falt
 541 -----------------------------------------------------------------------------*/
 542 void ParsecaseD( int argc, char** argv)
 543 {
 544   switch ( argv[CmndLinePos][ 2 ] )
 545     {
 546     case 'h':
 547       dhidden= 1;
 548       dmant = dmant -1;
 549       break;
 550     case 'e':
 551       if ( ++CmndLinePos != argc )
 552         {
 553           dexpon = atoi( argv[CmndLinePos] );
 554           if ( dexpon == 0 )
 555             {
 556               cout << "Error: not a number or zero" << endl;
 557               falt =1;
 558             }
 559         }
 560       break;
 561     case 't':
 562       if ( ++CmndLinePos != argc )
 563         {
 564           dmant = dmant + atoi( argv[CmndLinePos] );
 565           if ( dmant == 0 )
 566             cout << "Error: not a number or zero" << endl;
 567         }
 568       break;
 569     case 's':
 570       dexpon = 8;
 571       dmant = 23;
 572       dhidden = 1;
 573       break;
 574     case 'd':
 575       dexpon = 11;
 576       dmant = 52;
 577       dhidden = 1;
 578       break;
 579     case 'l':
 580       dexpon = 15;
 581       dmant = 64;
 582       dhidden = 0;
 583       break;
 584     case 'q':
 585       dexpon = 15;
 586       dmant = 112;
 587       dhidden = 1;
 588       break;
 589     case 'm':
 590       dexpon = 15;
 591       dmant = 240;
 592       dhidden = 0;
 593       break;
 594     default:
 595       if ( !def )
 596         {
 597           sexpon = 11;
 598           smant = 52;
 599           shidden = 1;
 600           def = 1;
 601         } else
 602           falt =1;
 603       break;
 604     } // switch
 605 }
 606
 607
 608 /*-----------------------------------------------------------------------------
 609 name        :ParsecaseL()
 610 description :long doulbe precision
 611 called from :ParseArgs()
 612 call fns    :\
 613 Global var. used :
 614 -----------------------------------------------------------------------------*/
 615 void ParsecaseL()
 616 {
 617   if ( !def )
 618     {
 619       sexpon = 15;
 620       smant = 64;
 621       shidden = 0;
 622       def = 1;
 623     } else
 624       falt =1;
 625 }
 626
 627
 628 /*-----------------------------------------------------------------------------
 629 name        :ParsecaseQ()
 630 description :quadruple precision
 631 called from :ParseArgs()
 632 call fns    :\
 633 Global var. used :
 634 -----------------------------------------------------------------------------*/
 635 void ParsecaseQ()
 636 {
 637   if ( !def )
 638     {
 639       sexpon = 15;
 640       smant = 112;
 641       shidden = 1;
 642       def = 1;
 643     } else
 644       falt =1;
 645 }
 646
 647
 648 /*-----------------------------------------------------------------------------
 649 name        :ParsecaseM()
 650 description :Multi precision
 651 called from :ParseArgs()
 652 call fns    :\
 653 Global var. used :
 654 -----------------------------------------------------------------------------*/
 655 void ParsecaseM()
 656 {
 657   if ( !def )
 658     {
 659       sexpon = 15;
 660       smant = 240;
 661       shidden = 0;
 662       def = 1;
 663     } else
 664       falt =1;
 665 }
 666
 667
 668 /*-----------------------------------------------------------------------------
 669 name        :ParsecaseE()
 670 description :specify no. of bits for exponent
 671 called from :ParseArgs()
 672 call fns    :\
 673 Global var. used :CmndLinePos,falt
 674 -----------------------------------------------------------------------------*/
 675 void ParsecaseE( int argc, char** argv)
 676 {
 677   if ( !def )
 678     {
 679       if ( ++CmndLinePos != argc )
 680         {
 681           sexpon = atoi( argv[CmndLinePos] );
 682           if ( sexpon == 0 )
 683             {
 684               cout << "Error: not a number or zero" << endl;
 685               falt =1;
 686             }
 687           defe =1;
 688         } else
 689           falt = 1;
 690     } else
 691       falt =1;
 692 }
 693
 694
 695 /*-----------------------------------------------------------------------------
 696 name        :ParsecaseT()
 697 description :specify <int> bits precision
 698 called from :ParseArgs()
 699 call fns    :\
 700 Global var. used :CmndLinePos,falt
 701 -----------------------------------------------------------------------------*/
 702 void ParsecaseT( int argc, char** argv)
 703 {
 704   if( !def )
 705     {
 706       if ( ++CmndLinePos != argc )
 707         {
 708           smant = smant + atoi( argv[CmndLinePos] );
 709           if ( smant == 0 )
 710             {
 711               cout << "Error: not a number or zero" << endl;
 712               falt =1;
 713             }
 714           deft=1;
 715         } else
 716           falt =1;
 717     } else
 718       falt =1;
 719 }
 720
 721
 722 /*-----------------------------------------------------------------------------
 723 name        :ParsecaseR()
 724 description :specify rounding modes
 725 called from :ParseArgs()
 726 call fns    :\
 727 Global var. used :
 728 -----------------------------------------------------------------------------*/
 729 void ParsecaseR( int argc, char** argv)
 730 {
 731   int j;
 732   Round = 0;
 733   if ( ++CmndLinePos != argc )
 734     {
 735       if ( ( argv[CmndLinePos][ 0 ]  != '-' )  && ( CmndLinePos != argc -1 ) )
 736         {
 737           j=0;
 738           while ( argv[CmndLinePos][ j ]  != '\0' )
 739             {
 740               switch( argv[CmndLinePos][ j ] )
 741                 {
 742                 case 'n':
 743                   Round |= ROUND_NEAREST;
 744                   break;
 745                 case 'p':
 746                   Round |= ROUND_UP;
 747                   break;
 748                 case 'm':
 749                   Round |= ROUND_DOWN;
 750                   break;
 751                 case 'z':
 752                   Round |= ROUND_ZERO;
 753                   break;
 754                 }
 755               j++;
 756             } // while
 757         }
 758       else
 759         {
 760           Round = ROUND_ALL;     // Test all
 761           CmndLinePos--;
 762         }
 763     } else
 764       falt = 1;
 765 }
 766
 767
 768 /*-----------------------------------------------------------------------------
 769 name        :ParsecaseN()
 770 description :do not test specified exceptions, denormalized nos., NaNs,
 771                  signed infinty, signed zeros
 772 called from :ParseArgs()
 773 call fns    :\
 774 Global var. used :
 775 -----------------------------------------------------------------------------*/
 776 void ParsecaseN( int argc, char** argv)
 777 {
 778   int j;
 779   if ( ++CmndLinePos != argc )
 780     {
 781       if ( ( argv[CmndLinePos][ 0 ]  != '-' )  && ( CmndLinePos != argc -1 ) )
 782         {
 783           j=0;
 784           while ( argv[CmndLinePos][ j ]  != '\0' )
 785             {
 786               switch( argv[CmndLinePos][ j ] )
 787                 {
 788                   /*
 789                     case 'i':
 790                     noFlags |= NO_FLAGS_INVALID;
 791                     break;
 792                   */
 793                 case 'o':
 794                   noFlags |= NO_FLAGS_OVERFLOW ;
 795                   break;
 796                 case 'x':
 797                   noFlags |= NO_FLAGS_INEXACT;
 798                   break;
 799                 case 'z':
 800                   noFlags |= NO_FLAGS_DIV_BY_ZERO;
 801                   break;
 802                 case 'u':
 803                   noFlags |= NO_FLAGS_UNDERFLOW;
 804                   break;
 805                 case 't':
 806                   if ( ( argv[CmndLinePos][ j+1 ]  == 'i' )  &&
 807                        ( argv[CmndLinePos][ j+2 ]  == 'n' )  &&
 808                        ( argv[CmndLinePos][ j+3 ]  == 'y' ) )
 809                     {
 810                       tiny = 0;
 811                       j += 3;
 812                     } else
 813                       falt = 1;
 814                   break;
 815                 case 'i':
 816                   if ( ( argv[CmndLinePos][ j+1 ]  == 'n' )  &&
 817                        ( argv[CmndLinePos][ j+2 ]  == 'f' )  )  {
 818                     inf = 0;
 819                     j += 2;
 820                   } else
 821                     {
 822                       noFlags |= NO_FLAGS_INVALID;
 823                     }
 824                   break;
 825                 case 'n':
 826                   if ( ( argv[CmndLinePos][ j+1 ]  == 'a' )  &&
 827                        ( argv[CmndLinePos][ j+2 ]  == 'n' ) )
 828                     {
 829                       NaN = 0;
 830                       j += 2;
 831                     } else
 832                       falt = 1;
 833                   break;
 834                 case 's':
 835                   if ( ( argv[CmndLinePos][ j+1 ]  == 'n' )  &&
 836                        ( argv[CmndLinePos][ j+2 ]  == 'z' ) )
 837                     {
 838                       signedZero = 0;
 839                       j += 2;
 840                     } else
 841                       falt = 1;
 842                   break;
 843                 }
 844               j++;
 845             } // while
 846         }
 847       else
 848         {
 849           noFlags = NO_FLAGS;     //Test no flags (all values selected)
 850           CmndLinePos--;
 851         }
 852     } else
 853       falt = 1;
 854 }
 855
 856
 857 /*-----------------------------------------------------------------------------
 858 name        :ParsecaseJ()
 859 description :jump/skip test vectors raising overflow, underflow, invalid or
 860                 divide by zero exception
 861 called from :ParseArgs()
 862 call fns    :\
 863 Global var. used :
 864 -----------------------------------------------------------------------------*/
 865 void ParsecaseJ( int argc, char** argv)
 866 {
 867   int j;
 868   if ( ++CmndLinePos != argc )
 869     {
 870       if ( ( argv[CmndLinePos][ 0 ]  != '-' )  && ( CmndLinePos != argc -1 ) )
 871         {
 872           j=0;
 873           while ( argv[CmndLinePos][ j ]  != '\0' )
 874             {
 875               switch( argv[CmndLinePos][ j ] ) {
 876               case 'i':
 877                 jumpInvalid= 1;
 878                 break;
 879               case 'o':
 880                 jumpOverflow= 1;
 881                 break;
 882               case 'u':
 883                 jumpUnderflow= 1;
 884                 break;
 885               case 'z':
 886                 jumpDivZero= 1;
 887                 break;
 888               default:
 889                 falt = 1;
 890                 break;
 891               }
 892               j++;
 893             } // while
 894         }
 895     } else
 896       falt = 1;
 897 }
 898
 899
 900 /*-----------------------------------------------------------------------------
 901 name        :ParsecaseH()
 902 description :specify if hidden bit
 903 called from :ParseArgs()
 904 call fns    :\
 905 Global var. used :
 906 -----------------------------------------------------------------------------*/
 907 void ParsecaseH()
 908 {
 909   if( !def )
 910     {shidden= 1;
 911     smant = smant -1;}
 912   else
 913     falt =1;
 914 }
 915
 916
 917
 918
 919 /*-----------------------------------------------------------------------------
 920 name        :PostParseAction()
 921 description :if change to ucb then calls required fns,
 922                 else check for errors and opens log file.
 923 called from :ParseArgs()
 924 call fns    :ucb.OpenLogFile(),ucb.ReadLine(),ucb.DoLine(),ucb.Close
 925 Global var. used :check, def,defe,deft,flat,ucbInput
 926 -----------------------------------------------------------------------------*/
 927 void PostParseAction( int argc, char** argv )
 928 {
 929
 930   if ( !ucbInput )
 931     { // check float formats
 932       if ( ( sexpon > 32 )  || ( dexpon> 32 ) )
 933         {
 934           cout << "Exit: exponent size too big (> 32)" << endl;
 935           exit( 1 );
 936         } else if ( ( smant+shidden <24 )  || ( dmant+dhidden <24 )  || ( sexpon < 8 )  || ( dexpon < 8 ) )
 937           {
 938             cout << "Exit: floating-point format too small" << endl;
 939             exit( 1 );
 940           } else if ( ( smant + shidden < sexpon )  || ( dmant + dhidden < dexpon ) )
 941             {
 942               cout << "Exit: Precision t must at least equal exponent size" << endl;
 943               exit( 1 );
 944             } else if ( ( smant + shidden > ((unsigned long)1 << (sexpon-1)))  || ( dmant + dhidden > ((unsigned long)1 << (dexpon-1)) ) )
 945               {
 946                 cout << "Exit: Precision t > Bias + 1 = 2^(e-1)" << endl;
 947                 exit( 1 );
 948               }
 949     } // if
 950
 951   if ( ( !def )  & ( !deft | !defe ) )
 952     {
 953       cout << USAGE;
 954       exit( 1 );
 955     }
 956
 957
 958   if ( check )
 959     ucb.OpenLogFile( logfilename,argc,argv );
 960
 961
 962   if ( ucbInput )
 963     {
 964       if ( input )
 965       {
 966         ucb.OpenInput( argv[ argc-1 ] );
 967         while ( ucb.ReadLine( NULL,signedZero, noFlags ) )
 968           {
 969             ucb.DoLine( tiny,inf,NaN );
 970           }
 971         ucb.Close( tiny );
 972       }
 973       exit( 1 );
 974     }
 975
 976
 977   if ( ( input == NULL )  | ( ( output == NULL )  & ( check == -1 ) ) )
 978     {
 979       if ( input == NULL )
 980         cout << "Error: No input \n";
 981       else
 982         cout << "Error: No output\n";
 983
 984       cout << USAGE;
 985       exit( 1 );
 986     }
 987 }
 988
 989
 990
 991 /*-----------------------------------------------------------------------------
 992 name        :Calc_General_Stuff()
 993 description :calculates the sizes of arrays to store numbers according to the
 994                 precision specifications in the command line
 995 called from :main()
 996 call fns    :/
 997 exceptions  :/
 998 algorithm   :
 999 Global var. used :slength,srest,sexpon,smant,delength,drest,dexpon,dmant,
1000                         ssign_exp_length,ssign_exp_rest,skl_biased_exp,dkl_biased_exp,
1001                         dsign_exp_length,dsign_exp_rest.
1002 -----------------------------------------------------------------------------*/
1003
1004 void Calc_General_Stuff( )
1005 {
1006   //to calculate length of arrays used
1007   // change BV 17/01/01
1008
1009   /*  comments regarding the calculated variables
1010                 s for source , d for destination.
1011         slength ---- length of FP format in multiples of <no. of bits>
1012         srest   ---- no. of bits required for FP format in
1013                 the last of these allocated array elemnt
1014         ssign_exp_length ---number of array elements to
1015                 represent all the bits in exponent + sign bit
1016         ssign_exp_rest --- no. of bits required for sign + exponent in
1017                 the last of these allocated array elemnt - 1
1018         NO_OF_BITS--- number of bits in array element
1019   */
1020
1021   slength = ( unsigned long ) ceil( ( sexpon+smant+1 ) /float(NO_OF_BITS));
1022   srest = ( sexpon+smant+1 )  % NO_OF_BITS;
1023
1024   dlength = ( unsigned long ) ceil( ( dexpon+dmant+1 ) /float(NO_OF_BITS));
1025   drest = ( dexpon+dmant+1 )  % NO_OF_BITS;
1026
1027   ssign_exp_length = ( sexpon/NO_OF_BITS )  + 1;
1028   ssign_exp_rest = sexpon % NO_OF_BITS;
1029
1030   dsign_exp_length = ( dexpon/NO_OF_BITS )  + 1;
1031   dsign_exp_rest = dexpon % NO_OF_BITS;
1032
1033   skl_biased_exp = 1L << ( NO_OF_BITS-ssign_exp_rest-1 );
1034   dkl_biased_exp = 1L << ( NO_OF_BITS-dsign_exp_rest-1 );
1035
1036 }
1037
1038
1039
1040
1041
1042 /*-----------------------------------------------------------------------------
1043 name        :readAline()
1044 description :reads a test vector from the file and checks if the test
1045                 software/hardware is compatible with IEEE spcifications
1046 called from :main()
1047 call fns    :getbinary(),getNumber(),getdecimal()isspace(),mygetc(),myungetc(),
1048                 getModes(),getException(),skipSpace(),WriteNumber(),
1049                 initialize_values(),writeUCB().
1050 exceptions  :/
1051 algorithm   :/
1052 Global var. used :lots
1053 -----------------------------------------------------------------------------*/
1054
1055 int readAline( )
1056 {
1057   char c;
1058   err1=0,err2=0,err3=0, isIeeeVector=0;
1059   comp_sort=0;
1060   err_form=0;
1061   killerVector=0;
1062   int NaN1, NaN2, NaN3;
1063   int source=0;
1064
1065   char *tmp;
1066
1067   pre = post = 0;
1068   NaN1 = NaN2 = NaN3 = 0;
1069
1070   while ( ( ( c = mygetc( input ) )  == '!' )  || ( c == '\n' )  || ( c == ' ' ) )
1071     {
1072       if ( c == '!' )
1073         /* Skip rest of comment line. */
1074         while ( ( c = mygetc( input ) )  != '\n' )
1075           ;
1076       if ( c == '\n' )
1077         lines_in++;
1078     }
1079
1080   if ( c == EOF )
1081     return 0;
1082
1083   opcode = mygetc( input ); // 2nd character is the first character of the opcode
1084   op2code=op3code = op4code = ' ';
1085   testFormats=0;
1086   c = mygetc( input );
1087
1088   if ( !isspace( c ) )
1089     {
1090       op2code = c;
1091     }
1092   // seek further opcode bytes
1093   while ( !isspace( c ) ) c = mygetc( input ); // and throw them away because we only need first char
1094
1095
1096   // Get the precision specifications
1097   while ( isspace( c ) ) c= mygetc( input );
1098   switch ( c ) {
1099     case 's':
1100       op3code='s';
1101       break;
1102     case 'd':
1103       op3code='d';
1104       break;
1105     case 'l':
1106       op3code='l';
1107       break;
1108     case 'q':
1109       op3code='q';
1110       break;
1111     case 'm':
1112       op3code='m';
1113       break;
1114     case 'e':
1115       testFormats= TESTEVEN;
1116       break;
1117     case 'o':
1118       testFormats= TESTODD;
1119       break;
1120   default:
1121     op3code='z'; // op3code z-> no precisions specification
1122     myungetc( c, input); // put the current char back because it does not belong to precision specs
1123     break;
1124   }
1125   // check if test-vector precision specs correspond with command-line provided precision specs
1126   if ( op3code !='z' && op3code !='e' && op3code!='o' )
1127     {
1128       if ( ( (sexpon==8)  && ( (shidden+smant) ==24 )  && (op3code != 's') )  ||
1129            ( ( ( (sexpon!=8)  ||  (shidden+smant)!=24 ) )  && (op3code == 's') )  ||
1130            ( ( sexpon==11 )  && ( (shidden+smant)==53)  && (op3code != 'd') )   ||
1131            ( ( (sexpon!=11)  || ( (shidden+smant)!=53) )  && (op3code == 'd') ) ||
1132            ( ( sexpon==15 )  && ((shidden+smant)==64 )  && ( op3code != 'l' ) )  ||
1133            ( ( ( sexpon!=15 )  || ((shidden+smant)!=64 ) )  && ( op3code == 'l' ) )  ||
1134            ( ( sexpon==15 )  && ((shidden+smant)==113 )  && ( op3code != 'q' ) )  ||
1135            ( ( ( sexpon!=15 )  || ((shidden+smant)!=113 ) )  && ( op3code == 'q' ) )  ||
1136            ( ( sexpon==15 )  && ((shidden+smant)==240 )  && ( op3code != 'm' ) )  ||
1137            ( ( ( sexpon!=15 )  || ((shidden+smant)!=240 ) )   && ( op3code == 'm' ) )  )  {
1138         while ( ( c = mygetc( input ) )  != '\n' )  /* Skip rest of line. */
1139           ;
1140         return 1; // precision dependent test vectors
1141       }
1142     }
1143   c = mygetc( input );
1144   if ( op3code !='z' && op3code !='e' && op3code!='o' )
1145     {
1146       switch ( c )
1147         {
1148         case 'i':
1149           isIeeeVector= 1;
1150           while ( !isspace( c ) ) c= mygetc( input ); // read the eee and throw it away!
1151           break;
1152         }
1153     }
1154
1155   myungetc( c, input);
1156   skipSpace( );
1157   getModes( );
1158   skipSpace( );
1159
1160   if ( opcode == 'b' )
1161     {
1162         source=1;                                       //for source
1163         initialize_values(source);
1164       getbinary(operand1);
1165       for ( long i=0; i<slength;i++ )  // init operand2
1166         operand2[ i ] =0;
1167     } else if ( opcode == 'd' )  {
1168       getdecimal( );
1169       for ( long i=0; i<slength;i++ )  // init operand2
1170         operand2[ i ] =0;
1171     } else
1172       {
1173         source=1;                                       //for source
1174         initialize_values(source);
1175         getNumber(operand1);
1176
1177         if ( NaNresult )
1178           NaN1 = 1;
1179
1180         if ( wrong_input )
1181           err1=1;
1182
1183         skipSpace( );
1184         source=1;                                       //for source
1185         initialize_values(source);
1186         getNumber(operand2);
1187
1188         if ( NaNresult )
1189           NaN2 = 1;
1190
1191         if ( wrong_input )
1192           err2=1;
1193       }
1194
1195   skipSpace( );
1196   getExceptions( );
1197   skipSpace( );
1198
1199   if ( opcode == 'b' )
1200     {
1201       getdecimal( );
1202     }
1203   else if ( opcode == 'd' )
1204   {
1205         source=1;                                       //for source
1206           initialize_values(source);
1207         getbinary(result);
1208   }
1209        else
1210          {
1211         source=0;                                       //for destination
1212         initialize_values(source);
1213
1214            getNumber(result);
1215            if ( NaNresult )
1216              NaN3 = 1;
1217            if ( wrong_input )
1218              err3=1;
1219            if ( opcode != 'C' )
1220              comp_sort=6;
1221          }
1222   writeUCB();
1223   if ( err_form )
1224     printf( "\nWrong input on line %ld\n",lines_in+readings+1 );
1225   while ( c != '\n' )
1226     c = mygetc( input );
1227   return 1;
1228 }
1229
1230
1231 /*-----------------------------------------------------------------------------
1232   name        :writeUCB
1233   description :writes the test vector in UCB fromat
1234   called from :readAline()
1235   calls fns   :writeUCBopcode(),writeUCBprecision(),writeUCBround(),
1236                 writeUCBexcep(),writeUCBoperand,writeUCBresult()
1237                 writeUCBcomutaitive()
1238   Global var. used :/
1239   -----------------------------------------------------------------------------*/
1240 void writeUCB()
1241 {
1242   int m;
1243   if (  (  ( ( smant + shidden )  % 2 == 0 )  && ( testFormats == TESTEVEN ) )  |
1244         (  ( ( smant + shidden )  % 2 != 0 )  && ( testFormats == TESTODD ) )   |
1245         ( testFormats == ALL_FORMATS )  )
1246     {
1247       for ( m=0;m<4;m++ )
1248         {
1249           ostrstream stream;
1250
1251           if ( ( testModes >> m )  & 0x1 )
1252             {
1253                 writeUCBopcode(stream);
1254                 writeUCBprecision(stream);
1255                 writeUCBround(m,stream);
1256                 writeUCBexcep(stream);
1257                 writeUCBoperand(stream);
1258                 writeUCBresult(m,stream);
1259                 stream<<endl<<ends;
1260                 writeUCBcomutative(stream);
1261              }
1262         }
1263       if ( intstr != NULL )
1264         {
1265           delete [ ]  intstr;
1266           intstr = NULL;
1267         }
1268     }
1269
1270 }
1271
1272
1273 /*-----------------------------------------------------------------------------
1274   name        :writeUCBopcode()
1275   description :writes opcode of testvector in the UCB format
1276   called from :writeUCB
1277   calls fns   :
1278   Global var. used :/
1279   -----------------------------------------------------------------------------*/
1280 void writeUCBopcode(ostrstream &stream)
1281 {
1282
1283   switch ( opcode )  {
1284   case '+' :
1285     stream << "add";
1286     break;
1287   case '-' :
1288     stream << "sub";
1289     break;
1290   case '*' :
1291     stream << "mul";
1292     break;
1293   case '/' :
1294     stream << "div";
1295     break;
1296   case '%' :
1297     stream << "rem";
1298     break;
1299   case 'V' :
1300     stream <<"sqrt";
1301     break;
1302   case 'S' :
1303     stream <<"sqrt";
1304     break;
1305   case 'r':
1306     switch ( op2code )  {
1307     case ' ':
1308       stream <<"rt";
1309       break;
1310     case 'i':
1311       stream <<"ri";
1312       break;
1313     case 'u':
1314       stream <<"ru";
1315       break;
1316     case 'I':
1317       stream <<"rI";
1318       break;
1319     case 'U':
1320       stream <<"rU";
1321       break;
1322     } // switch
1323     break;
1324   case 'c':
1325     switch ( op2code )  {
1326     case ' ':
1327       stream <<"ct";
1328       break;
1329     case 'i':
1330       stream <<"ci";
1331       break;
1332     case 'u':
1333       stream <<"cu";
1334       break;
1335     case 'I':
1336       stream <<"cI";
1337       break;
1338     case 'U':
1339       stream <<"cU";
1340       break;
1341     } // switch
1342     break;
1343   case 'i':
1344     stream <<"i";
1345     break;
1346   case 'b':
1347     stream <<"b2d";
1348     break;
1349   case 'd':
1350     stream << "d2b";
1351     break;
1352   default:
1353     printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
1354     printf( "\nERROR: unrecognizable operation: %c\n", opcode );
1355     printf( "\n\nDone with %ld readings.\n", readings );
1356     exit( EXIT_SUCCESS );
1357   }
1358
1359 }
1360
1361 /*-----------------------------------------------------------------------------
1362   name        :writeUCBprecision()
1363   description :writes precision of testvector in the UCB format
1364   called from :writeUCB
1365   calls fns   :
1366   Global var. used :/
1367   -----------------------------------------------------------------------------*/
1368 void writeUCBprecision(ostrstream &stream)
1369 {
1370   if ( ( shidden )  && ( sexpon==8 )  && ( smant==23 ) )
1371     if ( !isIeeeVector )
1372       stream <<"s ";
1373     else
1374       stream <<"S ";
1375   else if ( ( shidden )  && ( sexpon==11 )  && ( smant==52 ) )
1376     if ( !isIeeeVector )
1377       stream <<"d ";
1378     else
1379       stream <<"D ";
1380   else if ( !( shidden )  && ( sexpon==15 )  && ( smant==64 ) )
1381     if ( !isIeeeVector )
1382       stream <<"l ";
1383     else
1384       stream <<"L ";
1385   else if ( ( shidden )  && ( sexpon==15 )  && ( smant==112 ) )
1386     if ( !isIeeeVector )
1387       stream <<"q ";
1388     else
1389       stream <<"Q ";
1390   else if ( !( shidden )  && ( sexpon==15 )  && ( smant==240 ) )
1391     if ( !isIeeeVector )
1392       stream <<"m ";
1393     else
1394       stream <<"M ";
1395   else
1396     stream << sexpon << " "<< shidden <<" "<< smant <<" " ;
1397
1398
1399   if ( ( ( opcode == 'r' )  || ( opcode == 'c' ) )  && ( op2code == ' ' ) )
1400     if ( ( dhidden )  && ( dexpon==8 )  && ( dmant==23 ) )
1401       stream <<"s ";
1402     else if ( ( dhidden )  && ( dexpon==11 )  && ( dmant==52 ) )
1403       stream <<"d ";
1404     else if ( !( dhidden )  && ( dexpon==15 )  && ( dmant==64 ) )
1405       stream <<"l ";
1406     else if ( ( dhidden )  && ( dexpon==15 )  && ( dmant==112 ) )
1407       stream <<"q ";
1408     else if ( !( dhidden )  && ( dexpon==15 )  && ( dmant==240 ) )
1409       stream <<"m ";
1410     else
1411       stream << dexpon << " "<< dhidden <<" "<< dmant <<" " ;
1412
1413 }
1414
1415
1416 /*-----------------------------------------------------------------------------
1417   name        :writeUCBround
1418   description :writes rounding mode of the test vector in the UCB format
1419   called from :writeUCB
1420   calls fns   :
1421   Global var. used :/
1422   -----------------------------------------------------------------------------*/
1423 void writeUCBround(int m,ostrstream &stream)
1424 {
1425   switch( m )  {
1426   case 0 :
1427     stream <<"n ";
1428     break;
1429   case 1 :
1430     stream <<"z ";
1431     break;
1432   case 2 :
1433     stream <<"p ";
1434     break;
1435   case 3 :
1436     stream <<"m ";
1437     break;
1438   }
1439 }
1440
1441
1442 /*-----------------------------------------------------------------------------
1443   name        :writeUCBexcep()
1444   description :writes exceptions of the testvector in the UCB format
1445   called from :
1446   calls fns   :
1447   Global var. used :/
1448   -----------------------------------------------------------------------------*/
1449 void writeUCBexcep(ostrstream &stream)
1450 {
1451   if ( !killerVector )
1452     stream <<"eq ";
1453   else
1454     stream <<"uo ";
1455
1456   if ( xcp==0 )
1457     stream <<"-";
1458   if ( xcp & long( TD_INVALID ) )
1459     stream <<"v";
1460   if ( xcp & long( TD_INEXACT ) )
1461     stream <<"x";
1462   if ( xcp & long( TD_OVERFLOW ) )
1463     stream <<"o";
1464   if ( xcp & long( TD_UNDERFLOW ) )
1465     stream <<"u";
1466   if ( xcp & long( TD_UNDERFLOWV ) )
1467     stream <<"a";
1468   if ( xcp & long( TD_UNDERFLOWW ) )
1469     stream <<"b";
1470   if ( xcp & long( TD_DIVBYZERO ) )
1471     stream <<"d";
1472 }
1473
1474
1475 /*-----------------------------------------------------------------------------
1476   name        :writeUCBoperand()
1477   description :writes operannds or source for the operation in UCB format
1478   called from :writeUCB
1479   calls fns   :
1480   Global var. used :/
1481   -----------------------------------------------------------------------------*/
1482 void writeUCBoperand(ostrstream &stream)
1483 {
1484 char *tmp;
1485   if ( ( intstr != NULL )  && ( opcode == 'c' ) )
1486     {
1487       stream << " ";
1488       stream << intstr;
1489       if ( NaN2 )
1490         NaNresult = 1;
1491       else
1492         NaNresult = 0;
1493       tmp = WriteNumber( operand2,sexpon,shidden,smant,slength,err1 );
1494       stream << tmp;
1495       delete [ ]  tmp;
1496
1497     } else
1498       {
1499
1500         if ( NaN1 )
1501           NaNresult = 1;
1502         else
1503           NaNresult = 0;
1504
1505         if ( opcode == 'd' )
1506           {
1507             stream << " ";
1508             stream << decimal;
1509           } else
1510             {
1511               tmp = WriteNumber( operand1,sexpon,shidden,smant,slength,err1 );
1512               stream << tmp;
1513               delete [ ]  tmp;
1514             }
1515
1516         if ( NaN2 )
1517           NaNresult = 1;
1518         else
1519           NaNresult = 0;
1520         tmp = WriteNumber( operand2,sexpon,shidden,smant,slength,err1 );
1521         stream << tmp;
1522         delete [ ]  tmp;
1523       }
1524 }
1525
1526
1527 /*-----------------------------------------------------------------------------
1528   name        :writeUCBresult()
1529   description :writes result or destination of the operation in UCB format
1530   called from :writeUCB
1531   calls fns   :writeUCBresultBin(),writeUCBresultDec(),writeUCBresultNum()
1532   Global var. used :/
1533   -----------------------------------------------------------------------------*/
1534 void writeUCBresult(int m,ostrstream &stream)
1535 {
1536   if ( ( intstr != NULL )  && ( opcode == 'r' ) )
1537     {
1538       stream << " ";
1539       if ( pre )
1540         stream << "?";
1541       stream << intstr;
1542       if ( post )
1543         stream << "?";
1544     } else if ( opcode == 'b' )
1545       {
1546         writeUCBresultBin(m,stream);
1547       }else if ( opcode == 'd' )
1548         {
1549           writeUCBresultDec(m,stream);
1550         }else
1551           {
1552             writeUCBresultNum(stream);
1553           }
1554 }
1555
1556
1557 /*-----------------------------------------------------------------------------
1558   name        :writeUCBresultBin()
1559   description :round and write the destination for binary to decimal conversion
1560   called from :writeUCBresult()
1561   calls fns   :
1562   Global var. used :/
1563   -----------------------------------------------------------------------------*/
1564 void writeUCBresultBin(int m,ostrstream &stream)
1565 {
1566   // check rounding mode
1567   stream << " ";
1568   switch( m )  {
1569   case 0 :
1570     switch ( nextdig )  {
1571     case 0:
1572     case 1:
1573     case 2:
1574     case 3:
1575     case 4:
1576       stream << decimal;
1577       break;
1578     default:
1579       stream << updecimal;
1580       break;
1581     }
1582     break;
1583   case 1 :
1584     stream << decimal;
1585     break;
1586   case 2 :
1587     if ( mysignbit )  // negative
1588       stream << decimal;
1589     else
1590       stream << updecimal;
1591     break;
1592   case 3 :
1593     if ( mysignbit )
1594       stream << updecimal;
1595     else
1596       stream << decimal;
1597     break;
1598   default:
1599     exit( 1 );
1600   }
1601 }
1602
1603 /*-----------------------------------------------------------------------------
1604   name        :writeUCBresultDec()
1605   description :writes deatination for decimal to binary conversion
1606   called from :writeUCBresult()
1607   calls fns   :
1608   Global var. used :/
1609   -----------------------------------------------------------------------------*/
1610 void writeUCBresultDec(int m,ostrstream &stream)
1611 {
1612 int source;
1613 char *tmp;
1614   // check rounding mode
1615   if ( pre )
1616     stream << "?";
1617
1618   if ( NaN3 )
1619     NaNresult = 1;
1620   else
1621     NaNresult = 0;
1622
1623   for ( long i = 0; i < dlength; i++ )
1624     operand1[ i ]  = result[ i ];
1625
1626   if ( nextdig )
1627     {
1628       source=0;                         //for destination
1629       initialize_values(source);
1630       Incr( operand1, 1);
1631     }
1632   switch( m )  {
1633   case 0 :
1634     switch ( nextdig )  {
1635     case 0:
1636     case 1:
1637       tmp = WriteNumber( result,dexpon,dhidden,dmant,dlength,err3 );
1638       stream << tmp;
1639       delete [ ]  tmp;
1640       break;
1641     default:
1642       tmp = WriteNumber( operand1,dexpon,dhidden,dmant,dlength,err3 );
1643       stream << tmp;
1644       delete [ ]  tmp;
1645       break;
1646     }
1647     break;
1648   case 1 :
1649     tmp = WriteNumber( result,dexpon,dhidden,dmant,dlength,err3 );
1650     stream << tmp;
1651     delete [ ]  tmp;
1652     break;
1653   case 2 :
1654     if ( mysignbit )
1655       {// negative
1656         tmp = WriteNumber( result,dexpon,dhidden,dmant,dlength,err3 );
1657         stream << tmp;
1658         delete [ ]  tmp;
1659       } else
1660         {
1661           tmp = WriteNumber( operand1,dexpon,dhidden,dmant,dlength,err3 );
1662           stream << tmp;
1663           delete [ ]  tmp;
1664         }
1665     break;
1666   case 3 :
1667     if ( mysignbit )
1668       {
1669         tmp = WriteNumber( operand1,dexpon,dhidden,dmant,dlength,err3 );
1670         stream << tmp;
1671         delete [ ]  tmp;
1672       } else
1673         {
1674           tmp = WriteNumber( result,dexpon,dhidden,dmant,dlength,err3 );
1675           stream << tmp;
1676           delete [ ]  tmp;
1677         }
1678     break;
1679   default:
1680     exit( 1 );
1681   }
1682   if ( post )
1683     stream << "?";
1684 }
1685
1686 /*-----------------------------------------------------------------------------
1687   name        :writeUCBresultNum()
1688   description :writes desination for operations other than b2d and d2b
1689   called from :writeUCBresult()
1690   calls fns   :
1691   Global var. used :/
1692   -----------------------------------------------------------------------------*/
1693 void writeUCBresultNum(ostrstream &stream)
1694 {
1695   char *tmp;
1696   if ( pre )
1697     stream << "?";
1698   if ( NaN3 )
1699     NaNresult = 1;
1700   else
1701     NaNresult = 0;
1702   tmp = WriteNumber( result,dexpon,dhidden,dmant,dlength,err3 );
1703
1704   stream << tmp;
1705   delete [ ]  tmp;
1706   if ( post )
1707     stream << "?";
1708 }
1709
1710 /*-----------------------------------------------------------------------------
1711   name        :writeUCBcomutative()
1712   description :if + or * then writes test vector to check comutativity and
1713                 calls functions to check the test vectors generated
1714   called from :writeUCB
1715   calls fns   :ucb.ReadLine(),ucb.DoLine()
1716   Global var. used :/
1717   -----------------------------------------------------------------------------*/
1718 void writeUCBcomutative(ostrstream &stream)
1719 {
1720   char* outputStream = stream.str( );
1721   string outputStream2;
1722   string op1,op2;
1723
1724   if (  opcode=='+' || opcode=='*' )
1725     {
1726       // if the operation is add or mult we have to
1727       // generate second line with operands switched places
1728       // in order to check correctness with concern to
1729       // commutativity law.
1730
1731       int i,j,k;
1732       outputStream2 = string( stream.str( ) );
1733
1734       switch (outputStream2[3])
1735         {
1736         case 's':
1737           i= FindNextSpace( outputStream2,10,1 );
1738           j= FindNextSpace( outputStream2,i+1,1 );
1739           op1= outputStream2.substr( i+1,j-i );
1740           k= FindNextSpace( outputStream2, j+1,1 );
1741           op2= outputStream2.substr( j+1,k-j );
1742           outputStream2.replace( i+1,op2.length( ),op2 );
1743           outputStream2.replace( i+1+op2.length( ),op1.length( ),op1 );
1744           break;
1745
1746         case 'd':
1747           i= FindNextSpace( outputStream2,10,1 );
1748           j= FindNextSpace( outputStream2,i+1,2 );
1749           op1= outputStream2.substr( i+1,j-i );
1750           k= FindNextSpace( outputStream2, j+1,2 );
1751           op2= outputStream2.substr( j+1,k-j );
1752           outputStream2.replace( i+1,op2.length( ),op2 );
1753           outputStream2.replace( i+1+op2.length( ),op1.length( ),op1 );
1754           break;
1755
1756         case 'q':
1757           i= FindNextSpace( outputStream2,10,1 );
1758           j= FindNextSpace( outputStream2,i+1,4 );
1759           op1= outputStream2.substr( i+1,j-i );
1760           k= FindNextSpace( outputStream2, j+1,4 );
1761           op2= outputStream2.substr( j+1,k-j );
1762           outputStream2.replace( i+1,op2.length( ),op2 );
1763           outputStream2.replace( i+1+op2.length( ),op1.length( ),op1 );
1764           break;
1765
1766         case 'l':
1767           i= FindNextSpace( outputStream2,10,1 );
1768           j= FindNextSpace( outputStream2,i+1,3 );
1769           op1= outputStream2.substr( i+1,j-i );
1770           k= FindNextSpace( outputStream2, j+1,3 );
1771           op2= outputStream2.substr( j+1,k-j );
1772           outputStream2.replace( i+1,op2.length( ),op2 );
1773           outputStream2.replace( i+1+op2.length( ),op1.length( ),op1 );
1774           break;
1775
1776         case 'm':
1777           i= FindNextSpace( outputStream2,10,1 );
1778           j= FindNextSpace( outputStream2,i+1,8 );
1779           op1= outputStream2.substr( i+1,j-i );
1780           k= FindNextSpace( outputStream2, j+1,8 );
1781           op2= outputStream2.substr( j+1,k-j );
1782           outputStream2.replace( i+1,op2.length( ),op2 );
1783           outputStream2.replace( i+1+op2.length( ),op1.length( ),op1 );
1784
1785         default:
1786           // mantsize and expsize were given explicitly
1787
1788           int opLength=smant+sexpon+1;
1789           int blocks;
1790
1791                                 // Calculate operand length as multiple of 32
1792           opLength= ( opLength/32 ) *32 + ( !( !( opLength%32 ) ) ) * 32;
1793
1794                                 // Calculate nr of 32 bit blocks for operand
1795           blocks=opLength/32;
1796
1797                                 // grep operand 1 out of the outputstream
1798           i= FindNextSpace( outputStream2,0,6 ); // find 6th space
1799           j= FindNextSpace( outputStream2,i+1,blocks );
1800           op1= outputStream2.substr( i+1,j-i );
1801
1802                                 // grep operand 2 out of the outputstream
1803           k= FindNextSpace( outputStream2, j+1, blocks );
1804           op2= outputStream2.substr( j+1,k-j );
1805
1806                                 // switch operand1 and operand2
1807           outputStream2.replace( i+1,op2.length( ),op2 );
1808           outputStream2.replace( i+1+op2.length( ),op1.length( ),op1 );
1809
1810           break;
1811         }
1812     }
1813   if ( check == 0 && !skipVector )
1814     {
1815       fprintf( output,outputStream );
1816       if (  opcode=='+' || opcode=='*' )
1817         {
1818           if (  !( op1==op2 )  )
1819             {
1820               // generate line with operands switched to test commutivity
1821               fprintf( output,outputStream2.c_str( ) );
1822             }
1823         }
1824     } else
1825       {
1826         if (!skipVector)
1827           {
1828             int tempLines=( int ) lines_in+readings+1;
1829
1830             ucb.ReadLine( outputStream,signedZero,noFlags,&tempLines );
1831             ucb.DoLine( tiny,inf,NaN );
1832
1833             if (  opcode=='+' || opcode=='*' )
1834               {
1835                 if (  !( op1==op2 )  )
1836                   {
1837                     // generate line with operands switched to test commutivity
1838                     ucb.ReadLine(  ( char* )  outputStream2.c_str( ), signedZero, noFlags, &tempLines );
1839                     ucb.DoLine( tiny,inf,NaN);
1840                   }
1841               }
1842           }
1843       }
1844   delete [ ]  outputStream;
1845 }
1846
1847
1848
1849
1850
1851 /*-----------------------------------------------------------------------------
1852 name        :putDot()
1853 description :Slam a dot out to show some life
1854 called from :readAline()
1855 call fns    :/
1856 exceptions  :/
1857 algorithm   :/
1858 Global var. used :/
1859 -----------------------------------------------------------------------------*/
1860
1861
1862 void putDot( )
1863 {
1864   cout << '.' << flush;
1865   ++dots;
1866   if ( dots > 78 )
1867     {
1868       cout << endl;
1869       dots = 0;
1870     }
1871 }
1872
1873
1874
1875
1876
1877 /*-----------------------------------------------------------------------------
1878 name        :FindNextSpace()
1879 description :to find location to stream with proper spaces
1880 called from :readAline()
1881 call fns    :/
1882 exceptions  :/
1883 algorithm   :
1884 Global var. used :/
1885 -----------------------------------------------------------------------------*/
1886
1887 int FindNextSpace(  string pstring , int beginpos, int occurance )
1888 {
1889
1890   int i=beginpos;
1891   while (  i<pstring.length( )  )
1892     {
1893       if (  pstring[ i ] ==' ' )
1894         {
1895           occurance--;
1896           if (  occurance==0 )
1897             return i;
1898         }
1899       i++;
1900     }
1901   return i;
1902 }
1903
1904
1905
1906 /*-----------------------------------------------------------------------------
1907 name        :skipSpace()
1908 description :avoid spaces while reading a vector
1909 called from :readAline()
1910 call fns    :mygetc(),myungetc()
1911 exceptions  :/
1912 algorithm   :/
1913 Global var. used :/
1914 -----------------------------------------------------------------------------*/
1915 void skipSpace( )
1916 {
1917   char c;
1918   while ( isspace( c = mygetc( input ) ) ) {
1919     if ( c=='\n' )
1920       lines_in++;
1921     ;
1922   }
1923   myungetc( c, input );
1924 }
1925
1926
1927
1928
1929
1930 /*-----------------------------------------------------------------------------
1931 name        :getModes()
1932 description :get rounding modes from the test vector
1933 called from :readAline()
1934 call fns    :mygetc(),myungetc()
1935 exceptions  :/
1936 algorithm   :
1937 Global var. used :round,testmodes,flushVector
1938 -----------------------------------------------------------------------------*/
1939
1940 void getModes( )
1941 {
1942   char c;
1943   testModes = flushVector = 0;
1944   while ( !isspace( c = mygetc( input ) ) )
1945     switch ( c )  {
1946       case 'A': // all rounding modes
1947         c = mygetc( input ); /* skip 'LL' */
1948       case 'U':
1949         if ( Round & ROUND_ZERO )
1950           testModes |= ( 1 << TD_TOWARDZERO );
1951         if ( Round & ROUND_NEAREST )
1952           testModes |= ( 1 << TD_TONEAREST );
1953         if ( Round & ROUND_DOWN )
1954           testModes |= ( 1 << TD_DOWNWARD );
1955         if ( Round & ROUND_UP )
1956           testModes |= ( 1 << TD_UPWARD );
1957         c = mygetc( input );
1958         break;
1959       case '0':
1960         if ( Round & ROUND_ZERO )
1961           testModes |= ( 1 << TD_TOWARDZERO );
1962         break;
1963       case '=':
1964         if ( Round & ROUND_NEAREST )
1965           testModes |= ( 1 << TD_TONEAREST );
1966         break;
1967       case '<':
1968         if ( Round & ROUND_DOWN )
1969           testModes |= ( 1 << TD_DOWNWARD );
1970         break;
1971       case '>':
1972         if ( Round & ROUND_UP )
1973           testModes |= ( 1 << TD_UPWARD );
1974         break;
1975       case 'o':
1976         testFormats = TESTODD;
1977         break;
1978       case 'e':
1979         testFormats = TESTEVEN;
1980         break;
1981       default:
1982         printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
1983         printf( "\nERROR: unrecognizable mode character: %c\n", c );
1984         printf( "\n\nDone with %ld readings.\n", readings );
1985         exit( EXIT_SUCCESS );
1986         break;
1987     }
1988   myungetc( c, input );
1989   if ( testFormats == 0 )
1990     testFormats = ALL_FORMATS;
1991 }
1992
1993
1994 /*-----------------------------------------------------------------------------
1995 name        :initialize_values()
1996 description :assigns values to global variables for precision parsed
1997                 from the command line.
1998 called from :readAline()
1999 call fns    :/
2000 exceptions  :/
2001 algorithm   :
2002 Global var. used :expon,hidden,mant,length,rest,sign_exp_length,
2003                 kl_biased_exp,exp_rest.
2004 -----------------------------------------------------------------------------*/
2005 void initialize_values( int source)
2006 {
2007   if (source ==1)
2008     {
2009       expon = sexpon;
2010       hidden = shidden;
2011       mant = smant;
2012       length = slength;
2013       rest = srest;
2014       sign_exp_length = ssign_exp_length;
2015       kl_biased_exp = skl_biased_exp;
2016       sign_exp_rest = ssign_exp_rest ;
2017     }
2018
2019   else                          //destnation
2020     {
2021       expon = dexpon;
2022       hidden = dhidden;
2023       mant = dmant;
2024       length = dlength;
2025       rest = drest;
2026       sign_exp_length = dsign_exp_length;
2027       kl_biased_exp = dkl_biased_exp;
2028       sign_exp_rest = dsign_exp_rest ;
2029     }
2030 }
2031
2032
2033
2034
2035
2036 /*-----------------------------------------------------------------------------
2037 name        :getExceptions()
2038 description :gets exceptions raised for the test vector
2039 called from :readAline()
2040 call fns    :
2041 exceptions  :/
2042 algorithm   :
2043 Global var. used :skipVector,jumpOverflow,jumpUnderflow,jumpInvalid,jumpDivZero
2044 -----------------------------------------------------------------------------*/
2045 void getExceptions( )
2046 {
2047   char c;
2048
2049   xcp = 0;
2050   while ( !isspace( c = mygetc( input ) ) ) {
2051     switch ( c ) {
2052       case 'x':
2053         xcp |= ( long )  TD_INEXACT;
2054         break;
2055         /*
2056          * In the Metro case, checking tininess before rounding to detect
2057          * underflow means that all possible cases of underflow apply.
2058          * (By contrast, checking AFTER rounding lets some boundary cases
2059          * slip by, and checking for EXTRAORDINARY error due to subnormalization
2060          * lets some tiny cases through that happen to round the same in
2061          * spite of subnormalization.
2062          */
2063
2064       case 'u':
2065         xcp |= ( long )  TD_UNDERFLOW;
2066         if (jumpUnderflow)
2067             skipVector=1;
2068         break;
2069       case 'v':
2070         xcp |= ( long )  TD_UNDERFLOWV;
2071         if (jumpUnderflow)
2072             skipVector=1;
2073         break;
2074       case 'w':
2075         xcp |= ( long )  TD_UNDERFLOWW;
2076         if (jumpUnderflow)
2077             skipVector=1;
2078         break;
2079       case 'o':
2080         xcp |= ( long )  TD_OVERFLOW;
2081         if (jumpOverflow)
2082             skipVector=1;
2083         break;
2084       case 'i':
2085         xcp |= ( long )  TD_INVALID;
2086         if (jumpInvalid)
2087             skipVector=1;
2088         break;
2089       case 'z':
2090         xcp |= ( long )  TD_DIVBYZERO;
2091         if (jumpDivZero)
2092             skipVector=1;
2093         break;
2094       case 'O':
2095         c = mygetc( input );
2096         /* skip the 'OK' */
2097         break;
2098       default:
2099         printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
2100         printf( "\nERROR: unrecognizable exception character: %c\n", c );
2101         printf( "\n\nDone with %ld readings.\n", readings );
2102         exit( EXIT_SUCCESS );
2103         break;
2104     }
2105   }
2106   myungetc( c, input );
2107 }
2108
2109
2110 /*-----------------------------------------------------------------------------
2111 name        :getbinary()
2112 description :reads the binary coded number form the test set //ONLY FOR CONVERSIONS
2113 called from :readAline()
2114 call fns    :mygetc(),myungetc(),skipSpaces()
2115 exceptions  :/
2116 algorithm   :
2117 Global var. used :length,rest,sign_exp_length,exp_rest,pre,binary,upbinary,nextdig
2118 -----------------------------------------------------------------------------*/
2119 void getbinary( unsigned long *number)
2120 {
2121   long i= 0;
2122   NaNresult = 0;
2123
2124   for ( i=0; i<length;i++ )  // init number
2125     number[ i ] =0;
2126
2127   mysignbit=getSign();
2128
2129   getBsignificant(number);
2130
2131   getBbias(number);
2132
2133 /* Stuff the sign nonarithmetically, to protect signaling NaNs. */
2134   if ( mysignbit )
2135     number[ 0 ]  |= 0x80000000L;
2136 }
2137
2138
2139 /*-----------------------------------------------------------------------------
2140   name        :getBsignificant()
2141   description :get significant of the binary number
2142   called from :getbinary()
2143   calls fns   :
2144   exceptions  :/
2145   algorithm   :
2146 Global var. used :/
2147   -----------------------------------------------------------------------------*/
2148 void getBsignificant(unsigned long *number)
2149 {
2150   int i,j,k=0,val=0;
2151   char c=mygetc(input);
2152   binary = new char[ maxstr ];
2153
2154   i = 0;
2155   while ( !isspace( c )  && ( c != '_' ) )
2156     {
2157       binary[ i++ ]  = c;
2158       c = mygetc( input );
2159     }
2160   nextdig = 0; // exact
2161   if ( c == '_' )
2162     {
2163       if ( mygetc( input )  == '1' )
2164         {
2165           nextdig = 2;
2166           if ( mygetc( input )  == '1' )
2167             nextdig++;
2168         } else
2169           nextdig = 1;
2170     }
2171   while ( !isspace( c ) )  c = mygetc( input ); // skip tail
2172   skipSpace( );
2173
2174   // copy binary per 4 bits for each hex input character
2175   i--; //  binary length (times 4 bits)
2176   j = length - 1;
2177
2178   while ( i - 7 >= 0 )
2179     {
2180       if ( isdigit( binary[ i-7 ] ) )
2181         val = binary[ i-7 ]  - '0';
2182       else
2183         val = 0xa + ( binary[ i-7 ]  - 'a' );
2184       number[ j ]  += val;
2185       for ( k = 6;k >= 0;k-- )
2186         {
2187           number[ j ]  *= 16;
2188           if ( isdigit( binary[ i-k ] ) )
2189             val = binary[ i-k ]  - '0';
2190           else
2191             val = 0xa + ( binary[ i-k ]  - 'a' );
2192           number[ j ]  += val;
2193         }
2194
2195       i -= 8;
2196       j--;
2197     }
2198   if ( i >= 0 )
2199     {
2200       for ( k = 0;k < i;k++ )
2201         {
2202           if ( isdigit( binary[ k ] ) )
2203             val = binary[ k ]  - '0';
2204           else
2205             val = 0xa + ( binary[ k ]  - 'a' );
2206           number[ j ]  += val;
2207           number[ j ]  *= 16;
2208         }
2209
2210       if ( isdigit( binary[ k ] ) )
2211         val = binary[ k ]  - '0';
2212       else
2213         val = 0xa + ( binary[ k ]  - 'a' );
2214       number[ j ]  += val;
2215     }
2216   if ( rest != 0 )
2217     {
2218       i = sign_exp_length - 1;
2219       for ( ; i<length-1;i++ )
2220         for ( j = 0;j < NO_OF_BITS - rest;j++ )
2221           {
2222             number[ i ]  *= 2;
2223             if ( number[ i+1 ]  & ( 1L << ( (NO_OF_BITS-1) - j ) ) )
2224               number[ i ] ++;
2225           }
2226       for ( j = 0;j < NO_OF_BITS - rest;j++ )
2227         number[ i ]  *= 2;
2228     }
2229
2230   if ( ( hidden )  && ( number[ 0 ]  != 0 ) )  // remove leading bit, if not zero!
2231     {
2232       if ( sign_exp_rest==(NO_OF_BITS-1) )
2233         number[ sign_exp_length ]  -= ( 1L<<(NO_OF_BITS-1) );
2234       else
2235         number[ sign_exp_length-1 ]  -= ( 1L<<( NO_OF_BITS-sign_exp_rest-1 ) );
2236     }
2237 }
2238
2239
2240 /*-----------------------------------------------------------------------------
2241   name        :getBbias()
2242   description :gets bias for the binary number
2243   called from :getbinary()
2244   calls fns   :getBexpon()
2245   exceptions  :/
2246   algorithm   :
2247 Global var. used :/
2248   -----------------------------------------------------------------------------*/
2249 void getBbias(unsigned long *number)
2250 {
2251   int i,j,k;
2252   long u=0;
2253
2254   u=getBexpon();
2255
2256   unsigned long * bias = new unsigned long[ sign_exp_length ];
2257   if ( bias == NULL )
2258     exitmem( );
2259   for ( i=1; i<sign_exp_length;i++ )
2260     bias[ i ] =0xFFFFFFFFL;
2261   bias[ 0 ] =0x3FFFFFFFL;
2262   if ( sign_exp_length > 1 )
2263     bias[ sign_exp_length-1 ]  <<= ( NO_OF_BITS-sign_exp_rest-1 );
2264   else
2265     bias[ 0 ] =( unsigned long ) ( ((unsigned long)1 << (expon-1)) -1 )  << ( NO_OF_BITS-sign_exp_rest-1 );
2266   // add exp to bias
2267   while ( u > 0 )
2268     { // positive exp
2269       bias[ sign_exp_length-1 ]  += ( 1L << ( NO_OF_BITS-sign_exp_rest-1 ) );
2270       k=sign_exp_length-1;
2271       while ( bias[ k ] ==0 )
2272         {
2273           k--;
2274           bias[ k ]  += 1;
2275         }
2276       k=sign_exp_length-1;
2277       u--;
2278     }
2279   while ( u < 0 )
2280     { // negative exp
2281       k=sign_exp_length-1;
2282       while ( bias[ k ] ==0 )
2283         {
2284           k--;
2285           bias[ k ]  -= 1;
2286         }
2287       bias[ sign_exp_length-1 ]  -= ( 1L << ( NO_OF_BITS-sign_exp_rest-1 ) );
2288       u++;
2289     }
2290
2291   // copy bias into number
2292   for ( i=0;i<sign_exp_length-1;i++ )
2293     number[ i ]  = bias[ i ];
2294   /*
2295     if (kl_biased_exp > 0) // change kl_biased_exp
2296     number[sign_exp_length-1] &= kl_biased_exp-1;
2297   */
2298   number[ sign_exp_length-1 ]  |= bias[ sign_exp_length-1 ];
2299   delete[ ]  bias;
2300
2301 }
2302
2303
2304 /*-----------------------------------------------------------------------------
2305   name        :getBexpon
2306   description :gets exponent for the binary number
2307   called from :getBbias()
2308   calls fns   :
2309   exceptions  :/
2310   algorithm   :
2311   Global var. used :/
2312   -----------------------------------------------------------------------------*/
2313 long int getBexpon()
2314 {
2315   char c;
2316   unsigned long signexp = 0;
2317   long u=0;
2318   c=mygetc( input ); // skip 'E' character
2319   c=mygetc( input );
2320   if ( c == '-' )
2321     {
2322       signexp = 1;
2323       c = mygetc( input );
2324     } else if ( c == '+' )
2325       {
2326         c = mygetc( input );
2327       }
2328
2329   while ( isdigit( c ) )
2330     {
2331       u = ( u * 10 )  + ( c - '0' );
2332       c = mygetc( input );
2333     } // while
2334   myungetc( c, input ); // for EOL
2335   if ( signexp )
2336     u = -u;
2337   return u;
2338 }
2339
2340
2341
2342
2343
2344 /*-----------------------------------------------------------------------------
2345 name        :getdecimal()
2346 description :get decimal coded number form testset  //  ONLY FOR CONVERSION
2347 called from :readAline()
2348 call fns    :getDsignificand(),getDexpon
2349 exceptions  :/
2350 algorithm   :
2351 Global var. used :decimal,updecimal,nextdig
2352 -----------------------------------------------------------------------------*/
2353
2354 void getdecimal( )
2355 {
2356   getDsignificand();            //i is the position till which
2357                                 //significant has occupied decimal[]
2358   getDexpon();
2359 }
2360
2361
2362
2363 /*-----------------------------------------------------------------------------
2364 name        :getDsignificand()
2365 description :get significand of decimal coded number
2366 called from :getdecimal()
2367 call fns    :mygetc(),myungetc()
2368 exceptions  :/
2369 algorithm   :
2370 Global var. used :decimal,updecimal,nextdig,incexp
2371 -----------------------------------------------------------------------------*/
2372
2373 void getDsignificand()
2374 {                                       // get significand
2375   posDec = 0,posUpDec = 0;
2376   incexp = 0;
2377   char c;
2378   c = mygetc( input );
2379
2380   while ( !isspace( c )  && ( c != '_' ) )
2381     {
2382       decimal[ posDec ]  = c;
2383       updecimal[ posDec++ ]  = c;
2384       c = mygetc( input );
2385     }
2386   if ( c == '_' )               //find next number for case of 9's in end
2387     {
2388       posUpDec = posDec-1;                      // or all 9's .
2389       while ( ( posUpDec > 1 )   && ( updecimal[ posUpDec ]  == '9' ) )
2390         updecimal[ posUpDec-- ]  = '0';
2391       if ( updecimal[ posUpDec ]  == '9' )
2392         {
2393           updecimal[ 1 ]  = '1';
2394           for ( posUpDec = 2; posUpDec <= posDec; posUpDec++ )
2395             updecimal[ posUpDec ]  = '0';
2396           incexp = 1;
2397         } else
2398           updecimal[ posUpDec ]  = ( int )  updecimal[ posUpDec ]  + 1;
2399       nextdig = mygetc( input )  - '0';
2400       c = mygetc( input );
2401       while ( !isspace( c ) )  c = mygetc( input ); // skip tail
2402     }  else
2403       {
2404         nextdig = 0;
2405       }
2406 }
2407
2408 /*-----------------------------------------------------------------------------
2409 name        :getDexpon()
2410 description :get exponent of decimal coded number
2411 called from :getdecimal()
2412 call fns    :mygetc(),myungetc(),exceptions  :/
2413 algorithm   :
2414 Global var. used :decimal,updecimal,nextdig,incexp
2415 -----------------------------------------------------------------------------*/
2416 void getDexpon()
2417 {
2418   char signexp,c;
2419   c = mygetc( input );
2420
2421   if ( c == 'E' )                       // get exponent
2422     {
2423       decimal[ posDec ]  = c;
2424       updecimal[ posDec++ ]  = c;
2425       c = mygetc( input );
2426       signexp = c;
2427       while ( !isspace( c ) )
2428         {
2429           decimal[ posDec ]  = c;
2430           updecimal[ posDec++ ]  = c;
2431           c = mygetc( input );
2432         }
2433       decimal[ posDec ]  = '\0';
2434       posUpDec = posDec;
2435
2436       if ( incexp )
2437         {
2438           if ( signexp != '-' )                 // inc exponent
2439             {
2440               posDec--;
2441               while ( updecimal[ posDec ]  == '9' )
2442                 {
2443                   updecimal[ posDec ]  = 0;
2444                   posDec--;
2445                 }
2446               if ( updecimal[ posDec ]  == 'E' )
2447                 {
2448                   updecimal[ posDec+1 ]  = '1';
2449                   updecimal[ posUpDec++ ]  = '0';
2450                 } else
2451                   {
2452                     updecimal[ posDec ]  = ( int )  updecimal[ posDec ]  + 1;
2453                   }
2454             } else                               // dec exponent
2455               {
2456                 posDec--;
2457                 while ( updecimal[ posDec ]  == '0' )
2458                   {
2459                     updecimal[ posDec ]  = 9;
2460                     posDec--;
2461                   }
2462                 updecimal[ posDec ]  = ( int )  updecimal[ posDec ]  - 1;
2463               }
2464         } // if
2465     }else
2466       {
2467         myungetc( c,input );
2468       }
2469   updecimal[ posUpDec ]  = '\0';
2470   myungetc( c, input );
2471 }
2472
2473
2474 /*-----------------------------------------------------------------------------
2475 name        :getNumber
2476 description :takes the coded number and precsion to get the number
2477 called from :readAline()
2478 call fns    :getSign(), getFPRootNumber(),getModifier()
2479 exceptions  :/
2480 algorithm   :
2481 Global var. used :HexNoIsGiven,cant_infinity,small_norm,NaNresult
2482 -----------------------------------------------------------------------------*/
2483
2484
2485 void getNumber( unsigned long *number)
2486 {
2487
2488   unsigned long i, k, pos;
2489   char c,capc,d,ch;
2490
2491   int signbit=0;
2492   HexNoIsGiven=0;               //initialize global variable
2493   cant_infinity=0, small_norm=0;        //initialize global variable
2494   NaNresult = 0;                        // change DV
2495
2496   for ( i=0; i<length;i++ )               // number initialized to 0
2497     number[ i ] =0;
2498
2499   signbit=getSign();
2500
2501   getFPRootNumber(number);
2502
2503   if (HexNoIsGiven) return ;
2504   /* if hexadecimal number is given as operand then save
2505      the number to the global char *intstr and exit this function */
2506
2507   getModifier(number);
2508
2509   /* Stuff the sign nonarithmetically, to protect signaling NaNs. */
2510   if (signbit)             // checksign() is to get the sign of sign of number
2511     number[ 0 ]  |= 0x80000000L;
2512 }
2513
2514
2515
2516 /*-----------------------------------------------------------------------------
2517 name        :WriteNumber()
2518 description :Inserts correct number of zeros at most significant posiotion.
2519 called from :readAline()
2520 call fns    :/
2521 Global var. used :length,nanresult
2522 -----------------------------------------------------------------------------*/
2523
2524
2525 char* WriteNumber( const unsigned long *number, unsigned long expon,
2526                    unsigned long hidden, unsigned long mant, unsigned long length,const int err )
2527 {
2528   unsigned long i;
2529   ostrstream stream;
2530
2531   if ( !err ) {
2532     for ( i=0; i<length;i++ )
2533       // stream << " "<< setw(8) << setfill('0') << hex << number[i];
2534       if ( number[ i ]  >= 0x10000000 )
2535         stream << " " << hex << number[ i ];
2536       else if ( number[ i ]  >= 0x01000000 )  {
2537         stream << " 0" << hex << number[ i ];
2538       } else if ( number[ i ]  >= 0x00100000 )  {
2539         stream << " 00" << hex << number[ i ];
2540       } else if ( number[ i ]  >= 0x00010000 )  {
2541         stream << " 000" << hex << number[ i ];
2542       } else if ( number[ i ]  >= 0x00001000 )  {
2543         stream << " 0000" << hex << number[ i ];
2544       } else if ( number[ i ]  >= 0x00000100 )  {
2545         stream << " 00000" << hex << number[ i ];
2546       } else if ( number[ i ]  >= 0x00000010 )  {
2547         stream << " 000000" << hex << number[ i ];
2548       } else {
2549         stream << " 0000000" << hex << number[ i ];
2550       }
2551   } else {
2552     stream<< " Fout invoernumber";
2553     err_form=1;
2554   }
2555   if ( NaNresult )
2556     stream << '?';
2557
2558   stream << ends;
2559   return stream.str( );
2560 }
2561
2562
2563
2564 /*-----------------------------------------------------------------------------
2565 name        :getSign()
2566 description :checks if the sign of a number is + or negative ... by default it returns + .
2567 called from :getNumber()
2568 call fns    :/                                //mygetc()
2569 exceptions  :/
2570 algorithm   :
2571 Global var. used :pre
2572 -----------------------------------------------------------------------------*/
2573 int getSign(){
2574   char c;
2575   int signBit=0;                //default takes positive
2576   c = mygetc( input );
2577   switch(c){
2578   case '-' :
2579     signBit=1;
2580     return signBit;
2581     break;
2582   case '+' :
2583     signBit=0;
2584     return signBit;
2585     break;
2586   case '?' :
2587     pre=1;
2588     break;
2589
2590   default :
2591     myungetc(c,input);
2592     return signBit;
2593     break;
2594   }
2595 }
2596
2597 /*-----------------------------------------------------------------------------
2598 name        :getFPRootNumber
2599 description :extracts the rootnumber from the coded number
2600 called from :getNumber()
2601 calls fns   :Infinity(),Quiet_NaN(),Signaling_NaN(),smallest_Norm(),getFPHexNo(),getFPInteger()
2602 exceptions  :/
2603 algorithm   :
2604 Global var. used :cant_infinity,killerVector,small_norm,NaNresult,HexNoIsGiven
2605 -----------------------------------------------------------------------------*/
2606 void getFPRootNumber(  unsigned long *number)
2607 {
2608   char c,d;
2609   c = mygetc( input );
2610   switch (c)
2611     {
2612     case 'H':
2613       Infinity( number);
2614       cant_infinity=1;
2615       break;
2616     case 'Q':
2617       Quiet_NaN( number);
2618       killerVector = 1;   /* change BV */
2619       NaNresult = 1; // change DV
2620       break;
2621     case 'S':
2622       Signaling_NaN( number);
2623       killerVector = 1;  /* change BV */
2624       NaNresult = 1;
2625       break;
2626     case 'E':              // this 'E' is given for compatibility with Coonen vectors
2627     case 'T':
2628       Smallest_Norm( number);
2629       small_norm=1;       // indien er een number volgt op T, dan kan via deze
2630       // variabele gezien worden of er een plus (voor T)
2631       // of een minus (voor H, zie aldaar) moet gebeuren
2632       break;
2633     case '0':
2634       d = mygetc( input );
2635       if ( d == 'x' )           //hexadecimal number
2636         {
2637           getFPHexNo();         //save Hexadecimal no. in intstr[]
2638           HexNoIsGiven=1;
2639         }
2640       else
2641         {
2642           myungetc(d,input);
2643           myungetc(c,input);
2644           getFPInteger( number);
2645         }
2646       break;
2647     default:
2648       if ( isdigit( c ) )
2649         {
2650           myungetc(c,input);
2651           getFPInteger( number);
2652         } else
2653           {
2654             printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
2655             printf( "\ngetNumber() ERROR: %c is an unknown number specifier in getNumber\n", c );
2656             printf( "\n\nDone with %ld readings.\n", readings );
2657             exit( EXIT_SUCCESS );
2658           }
2659     }
2660 }
2661
2662
2663 /*-----------------------------------------------------------------------------
2664   name        :getModifier
2665   description :Modifies the rootnumber
2666   called from :getNumber()
2667   calls fns   :ModifierPorM(),ModifierIorD(),ModifierU()
2668   exceptions  :/
2669   algorithm   :
2670 Global var. used :/
2671   -----------------------------------------------------------------------------*/
2672
2673
2674
2675 void getModifier(unsigned long *number)
2676 {
2677   unsigned long k;
2678   int McasePorM,McaseIorD;
2679   char c;
2680
2681   while ( !isspace( c = mygetc( input ) ) )
2682     {
2683       switch ( c )  {
2684       case 'p':
2685         McasePorM=0;
2686         ModifierPorM(number,McasePorM);
2687         break;
2688       case 'm' :
2689         McasePorM=1;
2690         ModifierPorM(number,McasePorM);
2691         break;
2692       case 'i' :
2693         McaseIorD=0;
2694         ModifierIorD(number,McaseIorD);
2695         break;
2696       case 'd' :
2697         McaseIorD=1;
2698         ModifierIorD(number,McaseIorD);
2699         break;
2700       case 'u' :
2701         k = getinteger( );
2702         ModifierU( number,k);
2703         break;
2704
2705
2706       case '0':
2707       case '1':
2708       case '2':
2709       case '3':             // indien er hier een cijfer wordt ingelezen, dan
2710
2711       case '4':             // kan dat enkel voorafgegaan zijn door een
2712
2713       case '5':             // H (infinity) of een T (smallest normalized)
2714
2715       case '6':             // if here a figure is read in , then that can
2716         // can only be preceded by H(infinty )
2717         //or T (smallest normalised)
2718       case '7':
2719             case '8':
2720       case '9':
2721       case 't':
2722       case 'B':
2723       case 'C':
2724         if ( ( c != 't' )  && ( c != 'u' )  && ( c!='B' )  && ( c != 'C' ) )
2725           k=( long ) atoi( &c );
2726         else if ( c == 't' )
2727           k = mant + hidden;
2728         else if ( c == 'B' )
2729           {
2730             // k = (unsigned long) (((unsigned long)1 << (expon-1)) - 1);
2731                   k = ULONG_MAX;
2732                   cout << "ULONG_MAX :" << ULONG_MAX << endl;
2733           } else if ( c == 'u' )
2734             k = dmant + dhidden;
2735         else if ( c == 'C' )
2736           // k = (unsigned long) (((unsigned long)1 << (dexpon-1)) - 1);
2737           k = ULONG_MAX;
2738         if ( cant_infinity )          // if H
2739         Minus(number,k);
2740                 //        PlusBk( number,k);
2741         if ( small_norm )        // if T
2742           Plus( number,k);
2743         break;
2744
2745
2746       default:
2747         printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
2748         printf( "\ngetNumber() ERROR: unrecognizable modifier character: %c\n", c );
2749         printf( "\n\nDone with %ld readings.\n", readings );
2750         exit( EXIT_SUCCESS );
2751       }
2752     }
2753   myungetc( c, input );
2754 }
2755
2756
2757
2758
2759
2760 /*-----------------------------------------------------------------------------
2761 name        :getFPHexNo()
2762 description :this is called if the coded number is hexadecimal . This stores hexadecimal number in intstr[] and quits getNumber().
2763 called from :getFPRootNumber()
2764 calls fns   :/
2765 exceptions  :/
2766 algorithm   :
2767 Global var. used : intstr[]
2768 -----------------------------------------------------------------------------*/
2769 void getFPHexNo(){
2770   long i;
2771   char c, capc ;
2772   intstr = new char[ 64 ];
2773   intstr[ 0 ]  = '0';
2774   intstr[ 1 ]  = 'x';
2775   i = 2;
2776   for ( capc = toupper( c = mygetc( input ) );
2777         ( '0' <= capc && capc <= '9' )  || ( 'A' <= capc && capc <= 'F' );
2778         capc = toupper( c = mygetc( input ) ) )
2779     {
2780       intstr[ i ]  = c;
2781       i++;
2782     }
2783   intstr[ i ]  = '\0';
2784  myungetc(c,input);
2785 }
2786
2787
2788
2789 /*-----------------------------------------------------------------------------
2790   name        :ModifierIorD
2791   description :Modifies the rootnumber for case 'i'-increment or case 'd'-decrement
2792   called from :getModifier()
2793   calls fns   :Incr_at_pos(), Incr(),Decr_at_pos(),Decr(),getInteger()
2794   exceptions  :/
2795   algorithm   :
2796   Global var. used :/
2797   -----------------------------------------------------------------------------*/
2798
2799 void ModifierIorD(unsigned long *number,int McaseIorD) {
2800   char c;
2801
2802   unsigned long pos,k;
2803   c = mygetc( input );
2804   switch ( c )  {
2805
2806   case '(':
2807     pos = getinteger( );
2808     c = mygetc( input );
2809     switch ( c )  {
2810
2811     case '+':
2812       pos += getinteger( );
2813       c = mygetc( input );
2814       if ( c == ')' )
2815         {
2816           k = getinteger( );
2817           if (McaseIorD)
2818             Decr_at_pos( number,pos,k);
2819           else
2820             Incr_at_pos( number,pos,k);
2821         } else
2822           {
2823             printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
2824             printf( "\ngetNumber() ERROR: unexpected character: %c\n", c );
2825             printf( "\n\nDone with %ld readings.\n", readings );
2826             exit( EXIT_SUCCESS );
2827                       } // end if
2828                   break;
2829
2830     case '-':
2831       pos -= getinteger( );
2832       c = mygetc( input );
2833       if ( c == ')' )
2834         {
2835           k = getinteger( );
2836           if (McaseIorD)
2837             Decr_at_pos( number,pos,k);
2838           else
2839             Incr_at_pos( number,pos,k);
2840         } else
2841           {
2842             printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
2843             printf( "\ngetNumber() ERROR: unexpected character: %c\n", c );
2844             printf( "\n\nDone with %ld readings.\n", readings );
2845             exit( EXIT_SUCCESS );
2846           } // end if
2847       break;
2848
2849     case ')':
2850       k = getinteger( );
2851       if (McaseIorD)
2852         Decr_at_pos( number,pos,k);
2853       else
2854         Incr_at_pos( number,pos,k);
2855       break;
2856
2857     default:
2858       printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
2859       printf( "\ngetNumber() ERROR: unexpected character: %c\n", c );
2860       printf( "\n\nDone with %ld readings.\n", readings );
2861       exit( EXIT_SUCCESS );
2862     } // end if
2863     break;
2864
2865   default:
2866     myungetc( c, input );
2867     k = getinteger( );
2868     if (McaseIorD)
2869       Decr( number,k);
2870     else
2871       Incr( number,k);
2872     break;
2873   } // switch case 'i'
2874
2875 }
2876
2877
2878
2879 /*-----------------------------------------------------------------------------
2880   name        :ModifierPorM
2881   description :Modifies the rootnumber for case 'p'-plus or case 'm'-minus
2882   called from :getModifier()
2883   calls fns   :Plus(),PLusBk(),Minus(),MinusBk()
2884   exceptions  :/
2885   algorithm   :
2886   Global var. used :/
2887   -----------------------------------------------------------------------------*/
2888 void ModifierPorM(unsigned long *number,int McasePorM)
2889 {
2890 char c,ch;
2891 unsigned long k=0;
2892
2893 c = mygetc( input );
2894   if (c== 'B' || c== 'C' )              //Bias (C is only for conversions)
2895     {
2896       ch=c;
2897       c= mygetc(input);
2898       if ((c - '0')>0 && (c - '0')<10)
2899         {
2900           k=c-'0' ;
2901           if (McasePorM)
2902             MinusBk( number,k);
2903           else
2904             PlusBk( number,k);
2905           return;
2906         } else  if ((c - '0') != 0 )    //makes B0 as B
2907           myungetc(c,input) ;
2908         c=ch;
2909     }
2910   myungetc( c,input );
2911   k = getinteger( );
2912   if (McasePorM)
2913     Minus(number,k);
2914   else
2915     Plus(number,k);
2916 }
2917
2918
2919
2920
2921 /*-----------------------------------------------------------------------------
2922 name        :ModifierU()
2923 description :operates u<digit>
2924 called from :getModifier()
2925 call fns    :null_exp(),one_exp()
2926 exceptions  :/
2927 algorithm   :
2928 Global var. used :sign_exp_length,exp_rest,length,rest,hidden,kl_biased_exp,mant
2929 -----------------------------------------------------------------------------*/
2930
2931 void ModifierU( unsigned long *number, unsigned long ulps)
2932 {
2933   unsigned long *e, i, exp_all_null, bits_ulps, k, shift_part;
2934   int not_done;
2935
2936   e = new unsigned long[ sign_exp_length ];
2937
2938   if (  e == NULL ) {
2939     printf( "\nNot enough memory in procedure ModifierU!!!\n" );
2940     printf( "Program aborted.\n" );
2941     fclose( input );
2942     fclose( output );
2943     delete[ ]  operand1;
2944     delete[ ]  operand2;
2945     delete[ ]  result;
2946     printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
2947     printf( "\n\nDone with %ld readings.\n", readings );
2948     exit( EXIT_SUCCESS );
2949   }
2950
2951   for ( i=0;i<sign_exp_length-1;i++ )
2952     e[ i ] =number[ i ];
2953   e[ sign_exp_length-1 ]  = ( number[ sign_exp_length-1 ]  >> ( (NO_OF_BITS-1)-sign_exp_rest ) );
2954   e[ sign_exp_length-1 ]  <<= ( (NO_OF_BITS-1)-sign_exp_rest );
2955
2956
2957
2958   if (  null_exp( e)  ||
2959         one_exp( e)  )
2960     {
2961       for ( i=0;i<length;i++ )
2962         number[ i ] =0;
2963       if ( ulps != 0 )
2964         {
2965           bits_ulps = ( unsigned long )  floor( ::log10( ulps ) / ::log10( 2 ) ) +1;
2966           if ( bits_ulps > mant )
2967             {
2968               ulps=0;
2969               wrong_input=1;
2970             }
2971         }
2972       number[ length-1 ] =ulps;
2973
2974       if ( rest!=0 )
2975         {
2976           if ( ( ( bits_ulps > rest )  && hidden )  || ( ( bits_ulps+1 > rest )  && ( !hidden ) ) )
2977             number[ length-2 ]  = ( ulps >> rest );
2978           number[ length-1 ]  <<= ( NO_OF_BITS-rest );
2979         }
2980     }
2981   else
2982     {
2983
2984       for ( i=0;i<length;i++ )
2985         number[ i ] =0;
2986       if ( ulps != 0 )
2987         {
2988           bits_ulps = ( unsigned long )  floor( ::log10( ulps ) / ::log10( 2 ) ) +1;
2989         if ( bits_ulps > mant )
2990           {
2991             ulps=0;
2992             wrong_input=1;
2993           }
2994         }
2995       number[ length-1 ] =ulps;
2996
2997       if ( rest!=0 )
2998         {
2999           if ( ( ( bits_ulps > rest )  && hidden )  || ( ( bits_ulps+1 > rest )  && ( !hidden ) ) )
3000             number[ length-2 ]  = ( ulps >> rest );
3001           number[ length-1 ]  <<= ( NO_OF_BITS-rest );
3002         }
3003
3004       k=sign_exp_length-1;
3005       shift_part=length-1;
3006
3007       while ( ( hidden && !( number[ sign_exp_length-1 ]  & kl_biased_exp ) )  ||
3008               ( !hidden &&  (  ( ( !( number[ sign_exp_length ]  & ( 1L<<(NO_OF_BITS-1) ) ) )  &&
3009                                  ( sign_exp_rest == (NO_OF_BITS-1) ) )  ||
3010                                ( !( number[ sign_exp_length-1 ]  & ( kl_biased_exp>>1 ) )  &&
3011                                  ( sign_exp_rest != (NO_OF_BITS-1) ) )  ) )
3012               )  {
3013         e[ sign_exp_length-1 ]  -= ( 1L << ( (NO_OF_BITS-1)-sign_exp_rest ) );
3014         k=sign_exp_length-1;
3015         while ( k==sign_exp_length-1 ? e[ k ] ==( 0xFFFFFFFFL << ( NO_OF_BITS-sign_exp_rest-1 ) )  : e[ k ] ==0xFFFFFFFFL )
3016           {
3017             k--;
3018             e[ k ]  -= 1;
3019           }
3020         // shift left
3021         if ( shift_part!=0 )
3022           number[ shift_part-1 ]  = ( number[ shift_part-1 ]  << 1 )  | ( ( number[ shift_part ]  >> (NO_OF_BITS-1) )  & 0x01L );
3023         number[ shift_part ]  <<= 1;
3024         if ( number[ shift_part ] ==0 )
3025           shift_part--;
3026         i++;
3027       } // while
3028
3029       for ( i=0; i<sign_exp_length-1; i++ )
3030         number[ i ]  = e[ i ];
3031       number[ i ]  &= kl_biased_exp-1;
3032       number[ i ]  |= e[ sign_exp_length-1 ];
3033     }
3034
3035   delete[ ]  e;
3036 }
3037
3038
3039
3040
3041
3042 /*-----------------------------------------------------------------------------
3043 name        :getinteger()
3044 description :gets integer value from the string (also checks for hex no.)
3045 called from :Modifier(),ModifierIorD(),ModifierPorM()
3046 call fns    :mygetc(),myungetc()
3047 exceptions  :/
3048 algorithm   :
3049 Global var. used :
3050 -----------------------------------------------------------------------------*/
3051 /* Peruse the input stream for an unsigned integer.
3052  * Allow the 0xddddd syntax of C integers, and also decimals.
3053  */
3054
3055 unsigned long getinteger()
3056 {
3057   char c;
3058   unsigned long u = 0;
3059
3060   u=getHex();                   //returns integer value if Hex
3061   if (u != 0) return u;         // hex no. is given
3062
3063           //Now the input character c is either a decimal
3064           //digit or a terminal.  A null number returns as zero.
3065
3066   c=mygetc(input);
3067   switch(c){
3068   case 't':
3069     u = smant + shidden;
3070     break;
3071
3072
3073   case 'B' :
3074     u = ( unsigned long )  ( ((unsigned long)1 << (sexpon-1))  - 1 );
3075     break;
3076
3077   case 'u' :
3078     u = dmant + dhidden;
3079     break;
3080
3081
3082   case 'C' :
3083     u = ( unsigned long )  ( ((unsigned long)1 << (dexpon-1))  - 1 );
3084     // u = ULONG_MAX;
3085     break;
3086
3087   case 'h' :
3088     u = ( mant + hidden + 1 )  / 2 - 1;
3089     break;
3090
3091   default:
3092     while ( isdigit( c ) ) {
3093       u = ( u * 10 )  + ( c - '0' );
3094       c = mygetc( input );
3095     }
3096     myungetc( c, input );
3097   }
3098     return u;
3099 }
3100
3101
3102 /*-----------------------------------------------------------------------------
3103 name        :getHex()
3104 description :gets integer value of Hex no.
3105 called from :getinteger()
3106 call fns    :mygetc(),myungetc()
3107 Global var. used :
3108 -----------------------------------------------------------------------------*/
3109 unsigned long getHex()
3110 {
3111 char c,capc;
3112 unsigned long u=0;
3113 c = mygetc( input );
3114   if ( c == '0' ) {
3115     c = mygetc( input );
3116     if ( c == 'x' ) {
3117       for ( capc = toupper( c = mygetc( input ) );
3118             ( '0' <= capc && capc <= '9' )  || ( 'A' <= capc && capc <= 'F' ); capc = toupper( c = mygetc( input ) ) ) {
3119         u <<= 4;
3120         if ( isdigit( capc ) )
3121           u |= capc - '0';
3122         else
3123           u |= 0xa + ( capc - 'A' );
3124       }
3125       myungetc( c, input );
3126       return u;
3127     }
3128   }
3129 myungetc(c,input);
3130 return u;
3131
3132 }
3133
3134
3135 /*-----------------------------------------------------------------------------
3136 name        :getFPInteger()
3137 description :get floating ponit integer in IEEE format
3138 called from :getFPRootNumber()
3139 call fns    :mygetc(),myungetc()
3140 exceptions  :/
3141 algorithm   :
3142 Global var. used :wrong_input
3143 -----------------------------------------------------------------------------*/
3144
3145 void getFPInteger( unsigned long *number)
3146 {
3147   char c;
3148   unsigned long bits_u, ex, *bias, u=0;
3149   long k, i;
3150
3151   wrong_input=0;
3152   c=mygetc( input );
3153   while ( isdigit( c ) )
3154     {
3155       u = ( u * 10 )  + ( c - '0' );
3156       c = mygetc( input );
3157     } // while
3158
3159   myungetc( c, input );
3160
3161   if ( u != 0 )
3162     {  // determine no. of bits for u
3163       bits_u = ( unsigned long ) floor( ::log10( u ) / ::log10( 2 ) )  + 1;
3164       if ( bits_u> mant+hidden )
3165         {
3166           u=0;
3167           wrong_input=1;            // u greater than maximal value of significand
3168         }
3169
3170     }
3171   u &= 0x00FFFFFFL;                             // limited to 24-bit integer
3172   if ( u != 0 )
3173     {
3174       //determine no. of bits for u (now less than 24)
3175       bits_u = ( unsigned long ) floor( ::log10( u ) / ::log10( 2 ) )  + 1;
3176
3177       if ( !hidden )
3178         u |= ( 1L << bits_u );                   // put hidden bit in front
3179
3180       number[ length-1 ] =u;
3181
3182       if ( rest != 0 )
3183         {
3184           if ( ( ( bits_u > rest )  && hidden )  || ( ( bits_u+1 > rest )  && ( !hidden ) ) )
3185             number[ length-2 ] =( u >> rest );          // if the 2 parts overlap
3186                                                         // then the (bits_u-rest) bits that
3187                                                         //follow will be added to the previous part
3188           number[ length-1 ]  <<= ( NO_OF_BITS-rest );
3189         }
3190
3191       if ( hidden )
3192         ex=mant;
3193       else
3194         ex=mant-1;
3195
3196       k=length-1;
3197
3198       while ( ( number[ sign_exp_length-1 ]  & kl_biased_exp ) ==0 )
3199         {
3200           if ( k!=0 )
3201             number[ k-1 ]  = ( number[ k-1 ]  << 1 )  | ( ( number[ k ]  >> (NO_OF_BITS-1) )  & 0x01L );
3202
3203           number[ k ]  <<= 1;
3204           ex--;
3205
3206           if ( number[ k ] ==0 )
3207             k--;
3208         }
3209
3210       bias = new unsigned long[ sign_exp_length ];
3211
3212       if ( bias == NULL )
3213         {
3214           printf( "\nNot enough memory in procedure getFPInteger!!!\n" );
3215           printf( "Program aborted.\n" );
3216           fclose( input );
3217           fclose( output );
3218           delete[ ]  operand1;
3219           delete[ ]  operand2;
3220           delete[ ]  result;
3221           printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
3222           printf( "\n\nDone with %ld readings.\n", readings );
3223           exit( EXIT_SUCCESS );
3224         }
3225
3226       for ( i=0; i<sign_exp_length;i++ )
3227         bias[ i ] =0;
3228
3229       for ( i=1; i<sign_exp_length;i++ )
3230         bias[ i ] =0xFFFFFFFFL;
3231       bias[ 0 ] =0x3FFFFFFFL;
3232
3233       if ( sign_exp_length > 1 )
3234         bias[ sign_exp_length-1 ]  <<= ( NO_OF_BITS-sign_exp_rest-1 );
3235       else
3236         bias[ 0 ] =( unsigned long ) ( ((unsigned long)1 << (expon-1)) -1 )  << ( NO_OF_BITS-sign_exp_rest-1 );
3237
3238                                       //add rest of exponent to bias
3239       k=sign_exp_length-1;
3240
3241       for ( i=0;i<ex;i++ )
3242         {
3243           bias[ sign_exp_length-1 ]  += ( 1L << ( NO_OF_BITS-sign_exp_rest-1 ) );
3244           while ( bias[ k ] ==0 )
3245             {
3246               k--;
3247               bias[ k ]  += 1;
3248             }
3249           k=sign_exp_length-1;
3250         }
3251
3252       // determine Exponent of number[] by assigning number[] =bias[]
3253       for ( i=0;i<sign_exp_length-1;i++ )
3254         number[ i ] =bias[ i ];
3255
3256       if ( kl_biased_exp > 0 )  // change kl_biased_exp
3257         number[ sign_exp_length-1 ]  &= kl_biased_exp-1;
3258
3259       number[ sign_exp_length-1 ]  |= bias[ sign_exp_length-1 ];
3260
3261       delete[ ]  bias;
3262     }
3263 }
3264
3265
3266
3267 /*-----------------------------------------------------------------------------
3268 name        :Infinity()
3269 description :Sets a number to infinty according to precision
3270 called from :getFPRootNumber()
3271 call fns    :
3272 exceptions  :/
3273 algorithm   :
3274 Global var. used :sign_exp_length,exp_rest,hidden
3275 -----------------------------------------------------------------------------*/
3276
3277
3278 void Infinity( unsigned long *number)
3279 {
3280   long j;
3281
3282   if ( sign_exp_length==1 )
3283     number[ 0 ] =0x7FFFFFFFL<<( NO_OF_BITS-sign_exp_rest-1 ) ;
3284   else
3285     number[ 0 ] =0x7FFFFFFFL;
3286
3287   if ( ( sign_exp_length==1 )  && ( sign_exp_rest!=(NO_OF_BITS-1) ) )
3288     number[ 0 ]  -= ( 1L << (NO_OF_BITS-1) );
3289
3290   for ( j=sign_exp_length-2;j>0;j-- )
3291     number[ j ]  += 0xFFFFFFFFL;
3292
3293   if ( sign_exp_length != 1 )
3294     number[ sign_exp_length-1 ]  = ( 0xFFFFFFFFL<<( NO_OF_BITS-sign_exp_rest-1 ) );
3295
3296   if ( !hidden ) {
3297     if ( sign_exp_rest==(NO_OF_BITS-1) )
3298       number[ sign_exp_length ]  += ( 1L<<(NO_OF_BITS-1) );
3299     else
3300       number[ sign_exp_length-1 ]  += ( 1L<<( NO_OF_BITS-sign_exp_rest-2 ) );
3301   }
3302 }
3303
3304
3305
3306 /*-----------------------------------------------------------------------------
3307 name        :Quiet_NaN()
3308 description :sets a number to Q
3309 called from :getFPRootNumber()
3310 call fns    :/
3311 exceptions  :/
3312 algorithm   :
3313 Global var. used :sign_exp_length,exp_rest,hidden
3314 -----------------------------------------------------------------------------*/
3315
3316
3317 void Quiet_NaN( unsigned long *number)
3318 {
3319   long j;
3320   if ( sign_exp_length==1 )
3321     number[ 0 ] =0x7FFFFFFFL<<( NO_OF_BITS-sign_exp_rest-1 );
3322   else
3323     number[ 0 ] =0x7FFFFFFFL;
3324
3325   if ( ( sign_exp_length==1 )  && ( sign_exp_rest!=(NO_OF_BITS-1) ) )
3326     number[ 0 ]  -= ( 1L << (NO_OF_BITS-1) );
3327
3328   for ( j=sign_exp_length-2;j>0;j-- )
3329     number[ j ]  += 0xFFFFFFFFL;
3330   if ( sign_exp_length != 1 )
3331     number[ sign_exp_length-1 ]  = ( 0xFFFFFFFFL<<( NO_OF_BITS-sign_exp_rest-1 ) );
3332
3333   if ( !hidden ) {
3334     if ( sign_exp_rest==(NO_OF_BITS-1) )
3335       number[ sign_exp_length ]  += ( 1L<<30 );
3336     else if ( sign_exp_rest==30 )
3337       number[ sign_exp_length ]  += ( 1L<<(NO_OF_BITS-1) );
3338     else
3339       number[ sign_exp_length-1 ]  += ( 1L<<( NO_OF_BITS-sign_exp_rest-3 ) );
3340   } else {
3341     if ( sign_exp_rest==(NO_OF_BITS-1) )
3342       number[ sign_exp_length ]  += ( 1L<<(NO_OF_BITS-1) );
3343     else
3344       number[ sign_exp_length-1 ]  += ( 1L<<( NO_OF_BITS-sign_exp_rest-2 ) );
3345   }
3346   if ( !hidden ) {
3347     if ( sign_exp_rest==(NO_OF_BITS-1) )
3348       number[ sign_exp_length ]  += ( 1L<<(NO_OF_BITS-1) );
3349     else
3350       number[ sign_exp_length-1 ]  += ( 1L<<( NO_OF_BITS-sign_exp_rest-2 ) );
3351   }
3352 }
3353
3354 /*-----------------------------------------------------------------------------
3355 name        :Signaling_NaN()
3356 description :sets number to S
3357 called from :getFPRootNumber()
3358 call fns    :/
3359 exceptions  :/
3360 algorithm   :
3361 Global var. used :sign_exp_length,exp_rest,hidden
3362 -----------------------------------------------------------------------------*/
3363
3364 void Signaling_NaN( unsigned long *number)
3365 {
3366   long j;
3367   if (  sign_exp_length==1 )
3368     number[ 0 ] =0x7FFFFFFFL<<( NO_OF_BITS-sign_exp_rest-1 );
3369   else
3370     number[ 0 ] =0x7FFFFFFFL;
3371
3372   if ( ( sign_exp_length==1 )  && ( sign_exp_rest!=(NO_OF_BITS-1) ) )
3373     number[ 0 ]  -= ( 1L << (NO_OF_BITS-1) );
3374   for ( j=sign_exp_length-2;j>0;j-- )
3375     number[ j ]  += 0xFFFFFFFFL;
3376   if ( sign_exp_length != 1 )
3377     number[ sign_exp_length-1 ]  = ( 0xFFFFFFFFL<<( NO_OF_BITS-sign_exp_rest-1 ) );
3378
3379   if ( rest )
3380     number[ length-1 ]  += ( 1L<<( NO_OF_BITS-rest ) );
3381   else
3382     number[ length-1 ]  += 1;
3383
3384   if ( !hidden ) {
3385     if ( sign_exp_rest==(NO_OF_BITS-1) )
3386       number[ sign_exp_length ]  += ( 1L<<(NO_OF_BITS-1) );
3387     else
3388       number[ sign_exp_length-1 ]  += ( 1L<<( NO_OF_BITS-sign_exp_rest-2 ) );
3389   }
3390 }
3391
3392
3393 /*-----------------------------------------------------------------------------
3394 name        :Smallest_Norm()
3395 description :sets number to samllest normal (or tiny T)
3396 called from :getFPRootNumber()
3397 call fns    :/
3398 exceptions  :/
3399 algorithm   :
3400 Global var. used :sign_exp_length,exp_rest,hidden
3401 -----------------------------------------------------------------------------*/
3402
3403 void Smallest_Norm( unsigned long *number)
3404 {
3405   // if (hidden) // change DV
3406   number[ sign_exp_length-1 ]  += ( 1L<<( NO_OF_BITS-sign_exp_rest-1 ) );
3407   if ( ( !hidden )  && ( sign_exp_rest==(NO_OF_BITS-1) ) )
3408     number[ sign_exp_length ]  += ( 1L<<(NO_OF_BITS-1) );
3409   else if ( !hidden )
3410     number[ sign_exp_length-1 ]  += ( 1L<<( NO_OF_BITS-sign_exp_rest-2 ) );
3411 }
3412
3413
3414
3415 /*-----------------------------------------------------------------------------
3416 name        :Plus()
3417 description :operates p<digit> operation in the coded number
3418 called from :ModifierPorM()
3419 call fns    :null_exp()
3420 exceptions  :/
3421 algorithm   :
3422 Global var. used :sign_exp_length,exp_rest
3423 -----------------------------------------------------------------------------*/
3424
3425 void Plus( unsigned long *number, unsigned long k)
3426 {
3427   unsigned long pre_zero, j, i;
3428   pre_zero = null_exp( number);
3429
3430   //bij de exponent 1'tje optellen en dit k keer
3431   for ( i=0;i<k;i++ ) {
3432     j = sign_exp_length-1;
3433     number[ j ]  += ( 1L<<( (NO_OF_BITS-1)-sign_exp_rest ) );
3434     if ( ( ( number[ j ] >>( NO_OF_BITS-sign_exp_rest-1 ) )  == 0 )  && ( j!=0 ) )      //enkel exponent_stuk
3435     {
3436       j--;
3437       number[ j ]  += 1;
3438     }
3439     if ( j>0 )
3440       if ( number[ j ] ==0 )  {
3441         j--;
3442         number[ j ]  += 1;
3443       } // if
3444   }
3445 }
3446
3447 /*-----------------------------------------------------------------------------
3448 name        :Minus()
3449 description :operates m<digits> operation in coded number
3450 called from :ModifierPorM
3451 call fns    :null_exp()
3452 exceptions  :/
3453 algorithm   :
3454 Global var. used :sign_exp_length,exp_rest
3455 -----------------------------------------------------------------------------*/
3456
3457 void Minus( unsigned long *number, unsigned long k)
3458 {
3459   unsigned long pre_zero, j, i;
3460
3461   pre_zero = null_exp( number);
3462   //bij de exponent 1'tje aftrekken en dit k keer
3463   // Only works for e<=31
3464   for ( i=0;i<k;i++ ) {
3465     j = sign_exp_length-1;
3466     number[ j ]  -= ( 1L<<( NO_OF_BITS-sign_exp_rest-1 ) );
3467     if ( ( number[ j ] >>( NO_OF_BITS-sign_exp_rest-1 ) )  == ( 0xFFFFFFFFL>>( NO_OF_BITS-sign_exp_rest-1 ) )  && ( j!=0 ) ) {
3468       j--;
3469       number[ j ]  -= 1;
3470     }
3471     while ( ( j>0 )  && ( number[ j ]  == 0xFFFFFFFFL ) ) {
3472       j--;
3473       number[ j ]  -= 1;
3474     }
3475   }
3476 }
3477
3478 /*-----------------------------------------------------------------------------
3479 name        :MinusBk()
3480 description :divide Bias by 2^k ... operation mB<digits>
3481 called from :ModifierPorM
3482 call fns    :/
3483 exceptions  :/
3484 algorithm   :
3485 Global var. used :
3486 -----------------------------------------------------------------------------*/
3487
3488
3489 void MinusBk( unsigned long *number, unsigned long k)
3490 {
3491   // incorrect, only for e <= 31    // for exp <= 2^31
3492   unsigned long Bk = 0;
3493   int i;
3494
3495   // Calc Bk
3496   Bk =  0x40000000L;
3497   /*   Divide (Bias + 1) by k */
3498   for ( i=1;i<=k;i++ )  // shift exponent k places to the right
3499     Bk = Bk >> 1;
3500
3501   // Sub Bk
3502   number[ 0 ]  -= Bk;
3503 }
3504
3505
3506 /*-----------------------------------------------------------------------------
3507 name        :PlusBk()
3508 description :multiply B by 2^k ... operation pB<digits>
3509 called from :ModifierPorM()
3510 call fns    :/
3511 exceptions  :/
3512 algorithm   :
3513 Global var. used :
3514 -----------------------------------------------------------------------------*/
3515
3516
3517 void PlusBk( unsigned long *number, unsigned long k)
3518 {
3519   // incorrect only works for e<=31   // for exp <= 2^31
3520   unsigned long Bk = 0;
3521   int i;
3522
3523   // Calc Bk
3524   Bk =  0x40000000L;
3525   /*   Divide (Bias + 1) by k */
3526   for ( i=1;i<=k;i++ )  // shift exponent k places to the right
3527     Bk = Bk >> 1;
3528
3529   // Add Bk
3530   number[ 0 ]  += Bk;
3531
3532 }
3533
3534
3535 /*-----------------------------------------------------------------------------
3536 name        :Incr_at_pos()
3537 description :operates i(<pos>)<digit>
3538 called from :ModifierIorD()
3539 call fns    :null_exp()
3540 exceptions  :/
3541 algorithm   :
3542 Global var. used :sign_exp_length,exp_rest,hidden,expon
3543 -----------------------------------------------------------------------------*/
3544
3545
3546 void Incr_at_pos( unsigned long *number, unsigned long pos, unsigned long k)
3547 {
3548   unsigned long leading, j, i,position,pos_part,pos_rest;
3549
3550   // determine leading bit if not hidden
3551   if ( ( sign_exp_rest!=(NO_OF_BITS-1) )  && ( !hidden ) )
3552     leading = ( number[ sign_exp_length-1 ] >>( NO_OF_BITS-sign_exp_rest-2 ) )  & 1L;
3553   else  if ( !hidden )
3554     leading = ( number[ sign_exp_length ] >>(NO_OF_BITS-1) )  & 1L;
3555
3556   // calculate position
3557   position = expon + pos + 1 - hidden;
3558   pos_part = position / NO_OF_BITS + 1;
3559   pos_rest = position % NO_OF_BITS + 1;
3560
3561   //Increment significand by 1, k times
3562   for ( i=0;i<k;i++ )  {
3563     if ( pos_rest==0 )  {
3564       number[ pos_part-2 ]  += 1;
3565       j = pos_part-3;
3566       while ( ( j>0 )  && ( number[ j ] ==0 ) ) {
3567         j--;
3568         number[ j ]  += 1;
3569       }
3570     } // if
3571     else {
3572       number[ pos_part-1 ]  += ( 1L<<( NO_OF_BITS-pos_rest ) );
3573       j = pos_part-1;
3574       if ( ( number[ j ] >>( NO_OF_BITS-pos_rest ) )  == 0 )  {
3575         j--;
3576         number[ j ]  += 1;
3577         while ( ( j>0 )  && ( number[ j ] ==0 ) ) {
3578           j--;
3579           number[ j ]  += 1;
3580         }
3581       } // if
3582     } // else
3583   } //for
3584
3585   // replace hidden bit
3586
3587   if ( ( !hidden )  && ( j<sign_exp_length ) ) {
3588     if ( ( sign_exp_rest==(NO_OF_BITS-1) )  && ( ( number[ sign_exp_length ] ) ==0 ) ) {
3589       if ( leading==0 )
3590         number[ sign_exp_length-1 ]  += 1;          //exp+1
3591       else if ( !null_exp( number) )
3592         number[ sign_exp_length ]  |= ( 1L<<(NO_OF_BITS-1) );     //leading 1
3593     }
3594     else if ( ( number[ sign_exp_length-1 ]  << ( sign_exp_rest+2 ) ) ==0 ) {
3595       if ( leading==0 )
3596         number[ sign_exp_length-1 ]  += 1L<<( NO_OF_BITS-sign_exp_rest-1 );      //exp+1
3597       else if ( !null_exp( number) )
3598         number[ sign_exp_length-1 ]  |= ( 1L<<( NO_OF_BITS-sign_exp_rest-2 ) );    //leading 1
3599     }
3600
3601   } else if ( ( !hidden )  && ( j==sign_exp_length ) ) {
3602     if ( ( sign_exp_rest==(NO_OF_BITS-1) )  && ( ( number[ sign_exp_length ] ) ==0 ) ) {
3603       if ( leading==0 )
3604         number[ sign_exp_length-1 ]  += 1;
3605       else if ( !null_exp( number) )
3606         number[ sign_exp_length-1 ]  |= ( 1L << (NO_OF_BITS-1) );
3607     }
3608   }
3609 }
3610
3611
3612 /*-----------------------------------------------------------------------------
3613 name        :Incr()
3614 description :operates i<digit> on the coded number
3615 called from :ModifierIorD()
3616 call fns    :null_exp()
3617 exceptions  :/
3618 algorithm   :
3619 Global var. used :length,hidden,sign_exp_length,exp_rest,rest
3620 -----------------------------------------------------------------------------*/
3621
3622 void Incr( unsigned long *number, unsigned long k)
3623 {
3624   unsigned long leidende, j, i;
3625
3626   //bepalen van leidende bit als deze niet hidden is
3627   if ( ( sign_exp_rest!=(NO_OF_BITS-1) )  && ( !hidden ) )
3628     leidende = ( number[ sign_exp_length-1 ] >>( NO_OF_BITS-sign_exp_rest-2 ) )  & 1L;
3629   else  if ( !hidden )
3630     leidende = ( number[ sign_exp_length ] >>(NO_OF_BITS-1) )  & 1L;
3631
3632   //1 optellen bij de significand en dit k keer
3633   for ( i=0;i<k;i++ ) {
3634     if ( rest==0 )
3635       number[ length-1 ]  += 1;
3636     else
3637       number[ length-1 ]  += ( 1L<<( NO_OF_BITS-rest ) );
3638     j = length-1;
3639     while ( ( j>0 )  && ( number[ j ] ==0 ) )  // carry propagation
3640     {
3641       j--;
3642       number[ j ]  += 1;
3643     }
3644   }
3645
3646   if ( ( !hidden )  && ( j<sign_exp_length ) ) {
3647     if ( ( sign_exp_rest==(NO_OF_BITS-1) )  && ( ( number[ sign_exp_length ] ) ==0 ) ) {
3648       if ( leidende==0 )
3649         number[ sign_exp_length-1 ]  += 1;          //exp+1
3650       else if ( !null_exp( number) )
3651         number[ sign_exp_length ]  |= ( 1L<<(NO_OF_BITS-1) );     //leidende 1
3652     }
3653     else if ( ( number[ sign_exp_length-1 ]  << ( sign_exp_rest+2 ) ) ==0 ) {
3654       if ( leidende==0 )
3655         number[ sign_exp_length-1 ]  += 1L<<( NO_OF_BITS-sign_exp_rest-1 );      //exp+1
3656       else if ( !null_exp( number ) )
3657         number[ sign_exp_length-1 ]  |= ( 1L<<( NO_OF_BITS-sign_exp_rest-2 ) );    //leidende 1
3658     }
3659
3660   } else if ( ( !hidden )  && ( j==sign_exp_length ) ) {
3661     if ( ( sign_exp_rest==(NO_OF_BITS-1) )  && ( ( number[ sign_exp_length ] ) ==0 ) ) {
3662       if ( leidende==0 )
3663         number[ sign_exp_length-1 ]  += 1;
3664       else if ( !null_exp( number) )
3665         number[ sign_exp_length-1 ]  |= ( 1L << (NO_OF_BITS-1) );
3666     }
3667   }
3668 }
3669
3670
3671 /*-----------------------------------------------------------------------------
3672 name        :Decr_at_pos()
3673 description :operates d(<pos>)<digit>  on the coded number
3674 called from :ModifierIorD()
3675 call fns    :null_exp()
3676 exceptions  :/
3677 algorithm   :
3678 Global var. used :hidden,sign_exp_length,exp_rest,expon
3679 -----------------------------------------------------------------------------*/
3680
3681
3682 void Decr_at_pos( unsigned long *number, unsigned long pos, unsigned long k)
3683 {
3684   unsigned long leading, j, i,position,pos_part,pos_rest;
3685
3686   // determine leading bit if not hidden
3687   if ( ( sign_exp_rest!=(NO_OF_BITS-1) )  && ( !hidden ) )
3688     leading = ( number[ sign_exp_length-1 ] >>( NO_OF_BITS-sign_exp_rest-2 ) )  & 1L;
3689   else  if ( !hidden )
3690     leading = ( number[ sign_exp_length ] >>(NO_OF_BITS-1) )  & 1L;
3691
3692   // calculate position
3693   position = expon + pos + 1 - hidden;
3694   pos_part = position / NO_OF_BITS + 1;
3695   pos_rest = position % NO_OF_BITS + 1;
3696
3697   for ( i=0;i<k;i++ )  {
3698     if ( pos_rest==0 )  {
3699       number[ pos_part-2 ]  -= 1;
3700       j = pos_part-3;
3701       while ( ( j>0 )  && ( number[ j ]  == 0xFFFFFFFFL ) )  {
3702         j--;
3703         number[ j ]  -= 1;
3704       } // while
3705     }
3706     else {
3707       number[ pos_part-1 ]  -= ( 1L<<( NO_OF_BITS-pos_rest ) );
3708       j = pos_part-1;
3709       if ( ( number[ j ] >>( NO_OF_BITS-pos_rest ) )  == ( 0xFFFFFFFFL>>( NO_OF_BITS-pos_rest ) ) )  {
3710         j--;
3711         number[ j ]  -= 1;
3712         while ( ( j>0 )  && ( number[ j ]  == 0xFFFFFFFFL ) )  {
3713           j--;
3714           number[ j ]  -= 1;
3715         } // while
3716       } // if
3717     } // else
3718   } // for
3719
3720   // bepalen van leidende bit als deze niet hidden is
3721   // change DV
3722
3723   if ( ( !hidden )  && ( !null_exp( number) ) )  {
3724     j = sign_exp_length-1;
3725
3726     if ( ( number[ j ]  >> ( 30-sign_exp_rest ) )  % 2 == 0 )  {
3727       j = sign_exp_length-1;
3728       number[ j ]  -= ( 1L<<( (NO_OF_BITS-1) - sign_exp_rest ) );
3729       if ( number[ j ] >>( (NO_OF_BITS-1)-sign_exp_rest )  == 0xFFFFFFFFL>>( (NO_OF_BITS-1)-sign_exp_rest ) ) {
3730         j--;
3731         number[ j ]  -= 1;
3732       }
3733       while ( ( j>0 )  && ( number[ j ]  == 0xFFFFFFFFL ) ) {
3734         j--;
3735         number[ j ]  -= 1;
3736       }
3737       // replace leading bit if necessary
3738       if ( sign_exp_rest==(NO_OF_BITS-1) )
3739         number[ sign_exp_length ]  |= 1L<<(NO_OF_BITS-1);
3740       else
3741         number[ sign_exp_length-1 ]  |= 1L<<( 30-sign_exp_rest );
3742     }
3743   } // if
3744 }
3745
3746 /*-----------------------------------------------------------------------------
3747 name        :Decr()
3748 description :operates d<digit> on the coded number
3749 called from :ModifierIorD()
3750 call fns    :null_exp()
3751 exceptions  :/
3752 algorithm   :
3753 Global var. used :length,hidden,sign_exp_length,exp_rest,rest
3754 -----------------------------------------------------------------------------*/
3755
3756
3757 void Decr( unsigned long *number, unsigned long k)
3758 {
3759   unsigned j, i;
3760   int carry = 0;
3761
3762   //decrement significand k times
3763   for ( i=0;i<k;i++ ) {
3764     j = length-1;
3765     if ( rest==0 )
3766       number[ j ]  -= 1;
3767     else {
3768       number[ j ]  -= ( 1L<<( NO_OF_BITS-rest ) );
3769       if ( ( number[ j ] >>( NO_OF_BITS-rest ) )  == ( 0xFFFFFFFFL>>( NO_OF_BITS-rest ) )  && ( j!=0 ) ) {
3770         j--;
3771         number[ j ]  -= 1;
3772       }
3773     }
3774     while ( ( j>0 )  && ( number[ j ]  == 0xFFFFFFFFL ) ) {
3775       j--;
3776       number[ j ]  -= 1;
3777     }
3778   }
3779
3780   // change DV
3781   if ( ( !hidden )  &&
3782        ( !null_exp( number) ) )  {
3783     // leading bit = zero ?
3784     j = sign_exp_length-1;
3785
3786     if (  ( ( sign_exp_rest < (NO_OF_BITS-1) )  && !( number[ j ]  & ( 1L << 30 - sign_exp_rest ) ) )  ||
3787           ( ( sign_exp_rest == (NO_OF_BITS-1) )  && !( number[ j ]  & ( 1L << (NO_OF_BITS-1) ) ) )  )  {
3788       // leading bit == 0
3789
3790       if ( sign_exp_rest == (NO_OF_BITS-1) )  {
3791         j++;
3792         number[ j ]  -=  1L;
3793         if ( ( j > 0 )  && ( number[ j ]  == 0xFFFFFFFFL ) )  {
3794           j--;
3795           number[ j ]  -= 1;
3796         }
3797       } else {
3798         /// change DV 28/01/2000
3799
3800         if ( hidden && ( !( ( number[ j ]  >> ( (NO_OF_BITS-1) - sign_exp_rest ) )  & 1L ) ) )  {
3801           // trailing exp bit == 0
3802           number[ j ]  += 1L << ( (NO_OF_BITS-1) - sign_exp_rest );
3803           j--;
3804           number [ j ]  -= 1;
3805
3806           while ( ( j>0 )  && ( number[ j ]  == 0xFFFFFFFFL ) ) {
3807             j--;
3808             number[ j ]  -= 1;
3809           }
3810         } else if ( sign_exp_rest == 0 )  { // change DV 8/3/2000
3811           if ( !( number[ j ]  & ( 1L << ( (NO_OF_BITS-1) - sign_exp_rest ) ) ) )  { // leading exp. bit = 0
3812             number[ j-1 ]  -= 1; // borrow carry
3813             number[ j ]  += ( 1L << ( (NO_OF_BITS-1) - sign_exp_rest ) );
3814           } else
3815             number[ j ]  -= ( 1L << ( (NO_OF_BITS-1) - sign_exp_rest ) );
3816
3817         }
3818         else {
3819           number[ j ]  -= ( 1L << ( (NO_OF_BITS-1) - sign_exp_rest ) );
3820
3821         }
3822
3823       }
3824     }
3825   }
3826   // replace leading bit if necessary
3827   if ( ( !hidden )  && ( !null_exp( number) ) )  {
3828     if ( sign_exp_rest==(NO_OF_BITS-1) )
3829       number[ sign_exp_length ]  |= 1L<<(NO_OF_BITS-1);
3830     else
3831       number[ sign_exp_length-1 ]  |= 1L<<( 30-sign_exp_rest );
3832   }
3833 }
3834
3835
3836 /*-----------------------------------------------------------------------------
3837 name        :exitmem()
3838 description :exit if less memory to store number
3839                                 //but why only getbinary() (have to figure out)
3840 called from :getbinary()
3841 call fns    :
3842 Global var. used :
3843 -----------------------------------------------------------------------------*/
3844
3845
3846 void exitmem( )
3847 {
3848   printf( "\nNot enough memory\n" );
3849   printf( "Program aborted.\n" );
3850   fclose( input );
3851   fclose( output );
3852   delete[ ]  operand1;
3853   delete[ ]  operand2;
3854   delete[ ]  result;
3855   printf( "\n\nError on line %ld of inputfile.",lines_in+readings+1 );
3856   printf( "\n\nDone with %ld readings.\n", readings );
3857   exit( 1 );
3858 }
3859
3860
3861
3862 /*-----------------------------------------------------------------------------
3863 name        :null_exp()
3864 description :checks if exponent is zero
3865 called from :Plus(),Minus(),Incr_at_pos(),Incr(),Decr_at_pos(),
3866                 Decr(),ModifierU().
3867 call fns    :/
3868 exceptions  :/
3869 algorithm   :
3870 Global var. used :sign_exp_length,exp_rest,length
3871 -----------------------------------------------------------------------------*/
3872
3873 unsigned long null_exp( const unsigned long *number)
3874 {
3875   unsigned long m, null=1;
3876
3877   for ( m=1;m<sign_exp_length-1;m++ ) {
3878     if ( number[ m ] !=0 )
3879       null = 0;
3880   }
3881   if ( ( sign_exp_rest==0 )  && ( ( ( number[ 0 ]  & 0x7FFFFFFFL )  << 1 )  != 0 ) )
3882     null=0;
3883   if ( ( ( number[ 0 ]  & 0x7FFFFFFFL )  >> ( NO_OF_BITS-sign_exp_rest-1 ) )  != 0 )
3884     null=0;
3885   if ( ( length != 1 )  && ( sign_exp_length !=1 )  && ( ( number[ sign_exp_length-1 ]  >> ( NO_OF_BITS-sign_exp_rest-1 ) )  != 0 ) )
3886     null=0;
3887
3888   return null;
3889 }
3890
3891 /*-----------------------------------------------------------------------------
3892 name        :one_exp()
3893 description :checks if exponent is one
3894 called from :ModifierU()
3895 call fns    :/
3896 exceptions  :/
3897 algorithm   :
3898 Global var. used :sign_exp_length,exp_rest
3899 -----------------------------------------------------------------------------*/
3900
3901
3902 unsigned long one_exp( const unsigned long *number)
3903 {
3904   unsigned long m, one=0;
3905
3906   if ( ( ( number[ sign_exp_length-1 ]  >> ( NO_OF_BITS - sign_exp_rest - 1 ) )  | 1L )  == 1L )
3907     one = 1;
3908
3909   for ( m=0;m<sign_exp_length-1;m++ ) {
3910     if ( number[ m ]  != 0 )
3911       one = 0;
3912   }
3913   if ( ( sign_exp_length > 1 )  && ( ( number[ 0 ]  & 0x7FFFFFFFL )  != 0 ) )
3914     one = 0;
3915
3916   return one;
3917 }
3918
3919
3920
3921
3922 /*-----------------------------------------------------------------------------
3923 name        :mygetc()
3924 description :apply getc() of <stdio.h> to the present location of file in use
3925 called from :whole program
3926 call fns    :/
3927 exceptions  :/
3928 algorithm   :/
3929 Global var. used :/
3930 -----------------------------------------------------------------------------*/
3931
3932 int mygetc( FILE *f )
3933 {
3934   int c;
3935   c = getc( f );
3936   return c;
3937 }
3938
3939
3940 /*-----------------------------------------------------------------------------
3941 name        :myungetc()
3942 description :apply ungetc() of <stdio.h> to the present location of file in use
3943 called from :whole program
3944 call fns    :/
3945 exceptions  :/
3946 algorithm   :/
3947 Global var. used :/
3948 -----------------------------------------------------------------------------*/
3949
3950 void myungetc( int c, FILE *f )
3951 {
3952   ungetc( c, f );
3953 }
3954
3955
3956
3957 /* ============================ END OF FILE ============================*/