home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / libsrc87 / norm.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-30  |  5.9 KB  |  216 lines
  1. /*
  2.  *  Ieee double normalization function
  3.  *
  4.  *  double norm( double d, int exp, int sign, int rbits )
  5.  *
  6.  *  inputs:
  7.  *    m    a 63 bit mantissa (passed as 8 bytes in a double's storage)
  8.  *    exp    BIASED exponent
  9.  *    sign    sign of the result (1 == -ve  0 == +ve)
  10.  *    rbits    8 bits of rounding info
  11.  *      on input radix point is at extreme right (exp should reflect this)
  12.  * output:
  13.  *    the normalized double in ieee double format
  14.  *
  15.  * algorithm:
  16.  *    see posting by p. housel housel@ecn.purdue.edu (thanks peter)
  17.  *        and  knuth v2 section 4.2.1    algorithm N
  18.  *
  19.  *    ++jrb    bammi@dsrgsun.ces.cwru.edu
  20.  * bugs:
  21.  *    strictly mc68X coding
  22.  *    (can be coded using gnu c `long long' type very easily)
  23.  */
  24. #include "flonum.h"
  25.  
  26. struct ldouble {
  27.     unsigned long hi, lo;
  28. };        /* yet another alias for a double */
  29.  
  30. #if 0
  31. double norm( double m, int exp, int sign, int rbits )
  32. #else
  33. double __normdf( double m, int exp, int sign, int rbits )
  34. #endif
  35. {
  36.     /* we'll use constraints later for values to end up in D reggies */
  37.     register unsigned long d0, d1;    /*  64 bit accumulator     */
  38.     register char rounding; /* rounding reg-signed char for convenience */
  39.     double result;            /* the normalized result    */
  40.     struct bitdouble *r = (struct bitdouble *)&result; /* alias    */
  41.  
  42.     /* init */
  43.     d0 = ((struct ldouble *)&m)->hi;
  44.     d1 = ((struct ldouble *)&m)->lo;
  45.     rounding = rbits;
  46.     
  47.     if( (d0 == 0) && (d1 == 0) && (rounding == 0) )
  48.     goto stuff_result; /* nothing to do, but stuff bits into result */
  49.  
  50.     /* keep shifting R & accumulating rounding bits until we have 53 bits */
  51.     /* lsb of rounding register is sticky for 1 bits shifted out          */
  52.  
  53.     asm volatile("ReNormalize:");    /* ok, i admit it :-) */
  54.     while( (d0 & 0xffe00000L) != 0)
  55.     {
  56. #if 0 /* only 5 operands allowed in asm() - why?? */
  57.     asm volatile (" lsrl    #1,%2 ; 
  58.                 roxrl    #1,%1 ;
  59.              roxrb    #1,%0 ;
  60.              jcc    NoCarry;
  61.              orb    #1,%0 ;
  62.              NoCarry:      "
  63.         : "=d" (rounding), "=d" (d1), "=d" (d0)       /* outputs */
  64.         : "0"  (rounding), "1"  (d1), "2"  (d0) ); /* inputs  */
  65. #else
  66.     asm volatile (" lsrl    #1,%1 ; 
  67.                 roxrl    #1,%0 "
  68.         : "=d" (d1), "=d" (d0)       /* outputs */
  69.         : "0"  (d1), "1"  (d0) ); /* inputs  */
  70.         asm volatile ("    roxrb    #1,%0 ;
  71.              jcc    NoCarry;
  72.              orb    #1,%0 ;
  73.              NoCarry:      "
  74.         : "=d" (rounding)       /* outputs */
  75.         : "0"  (rounding) ); /* inputs  */
  76. #endif
  77.     /* increment exponent to reflect the shift */
  78.     exp++;
  79.     }
  80.  
  81.     /* shift L until there is a 1 in the hidden bit pos, shifting
  82.        the rounding reggie into the lsb */
  83.     while( (d0 & 0x00100000L) == 0)
  84.     {
  85. #if 0 /* same here */
  86.     asm volatile ("    lslb    #1,%2 ;
  87.              roxll    #1,%1 ;
  88.              roxll    #1,%0 "
  89.         : "=d" (d0), "=d" (d1), "=d" (rounding)       /* outputs */
  90.         : "0"  (d0), "1"  (d1), "2"  (rounding) ); /* inputs  */
  91. #else
  92.     asm volatile ("    lslb    #1,%1 ;
  93.              roxll    #1,%0 "
  94.         : "=d" (d1), "=d" (rounding)       /* outputs */
  95.         : "0"  (d1), "1"  (rounding) );       /* inputs  */
  96.  
  97.     asm volatile ("    roxll    #1,%0 "
  98.         : "=d" (d0)       /* outputs */
  99.         : "0"  (d0));        /* inputs  */
  100. #endif
  101.     /* decrement exponent to reflect the shift */
  102.     --exp;
  103.     }
  104.  
  105.     /* check rounding register */
  106.     if (rounding < 0)
  107.     {    /* not round down */
  108.     if( (((unsigned char)rounding) & 0x80) == 0)
  109.     {    /* tie: round to even */
  110.         d1 &= 0xfffffffeL;    /* set lsb to 0 */
  111.     }
  112.     else
  113.     {    /* round up */
  114.         rounding = 0;
  115.     asm volatile ("
  116.         addql  #1,%1 ;
  117.         jcc    NoCarry1 ;
  118.         addql  #1,%0 ;
  119.             bra   ReNormalize; |just in case of overflow into hidden bit\n
  120.                 NoCarry1:           "
  121.         : "=d" (d0), "=d" (d1)       /* outputs */
  122.         : "0"  (d0), "1"  (d1) );  /* inputs  */
  123.     }
  124.     }
  125.     
  126.     /* exponent underflow ?? */
  127.     if(exp <= 0)
  128.     {
  129.     printf("Underflow %lx %lx %d\n", d0, d1, exp);
  130.     d0 = 0;
  131.     d1 = 0;
  132.     sign = 0;
  133.     exp = 0;
  134.     goto stuff_result;
  135.     }
  136.     /* exp overflow ? */
  137.     if(exp >= 2047)
  138.     {
  139.     /* cause overflow trap, but how ?? */
  140.     printf("Overflow %lx %lx %d\n", d0, d1, exp);
  141.     }
  142.     
  143.     stuff_result: /* stuff bit in result and ret */
  144.     r->sign = sign;
  145.     r->exp  = exp;
  146.     r->mant1 = d0;
  147.     r->mant2 = d1;
  148.     
  149.     return result;
  150. }
  151.  
  152. #ifdef TEST
  153. #include <stdio.h>
  154.  
  155. main()
  156. {
  157.     register unsigned long d0, d1;
  158.     double r1, r2, r3;
  159.     struct bitdouble *pr1 = (struct bitdouble *)&r1;
  160.     struct bitdouble *pr3 = (struct bitdouble *)&r3;
  161.     unsigned long *l = (unsigned long *)&r2;
  162.     
  163.     r2 = r1 = 3.14159265358979323846265e23;
  164.  
  165.     *l &= 0x000FFFFFL;    /* wallup sign and exponent  in r2 */
  166.     *l |= 0x00100000L;  /* add hidden bit               */
  167.     
  168.     /* try the straight case first */
  169.     r3 = norm(r2, (unsigned int)pr1->exp, (unsigned int)pr1->sign, 0);
  170.     
  171.     printf("%30.25e %05lx %08lx %4d %1d\n", r1, pr1->mant1, pr1->mant2,
  172.        pr1->exp, pr1->sign);
  173.     printf("%30.25e %05lx %08lx %4d %1d\n\n", r3, pr3->mant1, pr3->mant2,
  174.        pr3->exp, pr3->sign);
  175.  
  176.     /* shift L and try */
  177.     d0 = l[0];
  178.     d1 = l[1];
  179.     asm volatile ("    lsll    #1,%1 ;
  180.             roxll    #1,%0 "
  181.         : "=d" (d0), "=d" (d1)       /* outputs */
  182.         : "0"  (d0), "1"  (d1) );  /* inputs  */
  183.     l[0] = d0;
  184.     l[1] = d1;
  185.     
  186.     r3 = norm(r2, (unsigned int)pr1->exp - 1, (unsigned int)pr1->sign, 0);
  187.     
  188.     printf("%30.25e %05lx %08lx %4d %1d\n", r1, pr1->mant1, pr1->mant2,
  189.        pr1->exp, pr1->sign);
  190.     printf("%30.25e %05lx %08lx %4d %1d\n\n", r3, pr3->mant1, pr3->mant2,
  191.        pr3->exp, pr3->sign);
  192.  
  193.     /* how about a shift R */
  194.     r2 =r1;
  195.  
  196.     *l &= 0x000FFFFFL;    /* wallup sign and exponent  in r2 */
  197.     *l |= 0x00100000L;  /* add hidden bit               */
  198.  
  199.     d0 = l[0];
  200.     d1 = l[1];
  201.     asm volatile ("    lsrl    #1,%0 ;
  202.             roxrl    #1,%1 "
  203.         : "=d" (d0), "=d" (d1)       /* outputs */
  204.         : "0"  (d0), "1"  (d1) );  /* inputs  */
  205.     l[0] = d0;
  206.     l[1] = d1;
  207.     
  208.     r3 = norm(r2, (unsigned int)pr1->exp + 1, (unsigned int)pr1->sign, 0);
  209.     
  210.     printf("%30.25e %05lx %08lx %4d %1d\n", r1, pr1->mant1, pr1->mant2,
  211.        pr1->exp, pr1->sign);
  212.     printf("%30.25e %05lx %08lx %4d %1d\n\n", r3, pr3->mant1, pr3->mant2,
  213.        pr3->exp, pr3->sign);
  214. }
  215. #endif
  216.