home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / minix / libsrc~1.z / libsrc~1 / norm.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-12-28  |  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. double norm( double m, int exp, int sign, int rbits )
  31. {
  32.     /* we'll use constraints later for values to end up in D reggies */
  33.     register unsigned long d0, d1;    /*  64 bit accumulator     */
  34.     register char rounding; /* rounding reg-signed char for convenience */
  35.     double result;            /* the normalized result    */
  36.     struct bitdouble *r = (struct bitdouble *)&result; /* alias    */
  37.  
  38.     /* init */
  39.     d0 = ((struct ldouble *)&m)->hi;
  40.     d1 = ((struct ldouble *)&m)->lo;
  41.     rounding = rbits;
  42.     
  43.     if( (d0 == 0) && (d1 == 0) && (rounding == 0) )
  44.     goto stuff_result; /* nothing to do, but stuff bits into result */
  45.  
  46.     /* keep shifting R & accumulating rounding bits until we have 53 bits */
  47.     /* lsb of rounding register is sticky for 1 bits shifted out          */
  48.  
  49.     asm volatile("ReNormalize:");    /* ok, i admit it :-) */
  50.     while( (d0 & 0xffe00000L) != 0)
  51.     {
  52. #if 0 /* only 5 operands allowed in asm() - why?? */
  53.     asm volatile (" lsrl    #1,%2 ; 
  54.                 roxrl    #1,%1 ;
  55.              roxrb    #1,%0 ;
  56.              jcc    NoCarry;
  57.              orb    #1,%0 ;
  58.              NoCarry:      "
  59.         : "=d" (rounding), "=d" (d1), "=d" (d0)       /* outputs */
  60.         : "0"  (rounding), "1"  (d1), "2"  (d0) ); /* inputs  */
  61. #else
  62.     asm volatile (" lsrl    #1,%1 ; 
  63.                 roxrl    #1,%0 "
  64.         : "=d" (d1), "=d" (d0)       /* outputs */
  65.         : "0"  (d1), "1"  (d0) ); /* inputs  */
  66.         asm volatile ("    roxrb    #1,%0 ;
  67.              jcc    NoCarry;
  68.              orb    #1,%0 ;
  69.              NoCarry:      "
  70.         : "=d" (rounding)       /* outputs */
  71.         : "0"  (rounding) ); /* inputs  */
  72. #endif
  73.     /* increment exponent to reflect the shift */
  74.     exp++;
  75.     }
  76.  
  77.     /* shift L until there is a 1 in the hidden bit pos, shifting
  78.        the rounding reggie into the lsb */
  79.     while( (d0 & 0x00100000L) == 0)
  80.     {
  81. #if 0 /* same here */
  82.     asm volatile ("    lslb    #1,%2 ;
  83.              roxll    #1,%1 ;
  84.              roxll    #1,%0 "
  85.         : "=d" (d0), "=d" (d1), "=d" (rounding)       /* outputs */
  86.         : "0"  (d0), "1"  (d1), "2"  (rounding) ); /* inputs  */
  87. #else
  88.     asm volatile ("    lslb    #1,%1 ;
  89.              roxll    #1,%0 "
  90.         : "=d" (d1), "=d" (rounding)       /* outputs */
  91.         : "0"  (d1), "1"  (rounding) );       /* inputs  */
  92.  
  93.     asm volatile ("    roxll    #1,%0 "
  94.         : "=d" (d0)       /* outputs */
  95.         : "0"  (d0));        /* inputs  */
  96. #endif
  97.     /* decrement exponent to reflect the shift */
  98.     --exp;
  99.     }
  100.  
  101.     /* check rounding register */
  102.     if (rounding < 0)
  103.     {    /* not round down */
  104.     if( (((unsigned char)rounding) & 0x80) == 0)
  105.     {    /* tie: round to even */
  106.         d1 &= 0xfffffffeL;    /* set lsb to 0 */
  107.     }
  108.     else
  109.     {    /* round up */
  110.         rounding = 0;
  111.     asm volatile ("
  112.         addql  #1,%1 ;
  113.         jcc    NoCarry1 ;
  114.         addql  #1,%0 ;
  115.             bra   ReNormalize; |just in case of overflow into hidden bit\n
  116.                 NoCarry1:           "
  117.         : "=d" (d0), "=d" (d1)       /* outputs */
  118.         : "0"  (d0), "1"  (d1) );  /* inputs  */
  119.     }
  120.     }
  121.     
  122.     /* exponent underflow ?? */
  123.     if(exp <= 0)
  124.     {
  125.     printf("Underflow %lx %lx %d\n", d0, d1, exp);
  126.     d0 = 0;
  127.     d1 = 0;
  128.     sign = 0;
  129.     exp = 0;
  130.     goto stuff_result;
  131.     }
  132.     /* exp overflow ? */
  133.     if(exp >= 2047)
  134.     {
  135.     /* cause overflow trap, but how ?? */
  136.     printf("Overflow %lx %lx %d\n", d0, d1, exp);
  137.     }
  138.     
  139.     stuff_result: /* stuff bit in result and ret */
  140.     r->sign = sign;
  141.     r->exp  = exp;
  142.     r->mant1 = d0;
  143.     r->mant2 = d1;
  144.     
  145.     return result;
  146. }
  147.  
  148. #ifdef TEST
  149. #ifdef __MSHORT__
  150. #error "please run this test in 32 bit int mode"
  151. #endif
  152.  
  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.