アットウィキロゴ
sakat @ ウィキ
アットウィキロゴ
sakat @ ウィキ
ページ一覧

IEEE_754単精度の足し算を行うルーチン

  1. #define BITSHIFT 5  /* additional bit shift to avoid rounding errors */
  2.  
  3. /* add:  returns z which is the sum of x and y */
  4. void add(float *zp, const float *xp, const float *yp)
  5. {
  6.   unsigned long x, sx, ex, fx;
  7.   unsigned long y, sy, ey, fy;
  8.   unsigned long z, sz, ez, fz;
  9.   int d, p;
  10.  
  11.   /* get the bit pattern of each floting-point */
  12.   x = *(unsigned long *)xp;
  13.   y = *(unsigned long *)yp;
  14.  
  15.   /* sign */
  16.   sx = x >> 31;
  17.   sy = y >> 31;
  18.   /* exponential part */
  19.   ex = (x & 0x7fffffffUL) >> 23;
  20.   ey = (y & 0x7fffffffUL) >> 23;
  21.   /* fractional part */
  22.   fx = (x & 0x007fffffUL) | (1UL << 23);
  23.   fy = (y & 0x007fffffUL) | (1UL << 23);
  24.  
  25.   /* additional bit shift to avoid rounding errors */
  26.   fx <<= BITSHIFT;
  27.   fy <<= BITSHIFT;
  28.  
  29.   /* shift the bits of fractional part to conform the exponential part. 
  30.       This may cause a round error. */
  31.   d = ex - ey;
  32.   if (d > 0) {
  33.     if (d > 24)
  34.       fy = 0;
  35.     else
  36.       fy >>= d;
  37.     ez = ex;
  38.   } else if (d) { /* d < 0 */
  39.     d = -d;
  40.     if (d > 24)
  41.       fx = 0;
  42.     else
  43.       fx >>= d;
  44.     ez = ey;
  45.   } else {
  46.     ez = ex;
  47.   }
  48.  
  49.   /* convert negative numbers in two-complement number representation */
  50.   if (sx)
  51.     fx = ~fx + 1;
  52.   if (sy)
  53.     fy = ~fy + 1;
  54.  
  55.   /* obtain the sum of fx and fy */
  56.   fz = fx + fy;
  57.  
  58.   /* zero check */
  59.   if (fz) {
  60.  
  61.     /* obtain the sign and convert negative numbers into positive */
  62.     sz = fz >> 31;
  63.     if (sz)
  64.       fz = ~(fz - 1);
  65.  
  66.     /* adjust exponential part so that the 1st 1-bit occurs at 23rd position from the right */
  67.     for (d = 31; d >= 0 && ((fz >> d) & 1) == 0; d--)
  68.       ;
  69.     if ((p = d - 23 - BITSHIFT) > 0) {
  70.       fz >>= p;
  71.       ez += p;
  72.     } else if (p) {  /* p < 0 */
  73.       p = -p;
  74.       fz <<= p;
  75.       ez -= p;
  76.     }
  77.  
  78.     /* take off the additional bit shift and round the last bit */
  79.     fz = (fz >> BITSHIFT) + ((fz >> (BITSHIFT-1)) & 1);
  80.  
  81.     /* overflow check */
  82.     if (ez >= 0xffUL)
  83.       fz = 0UL;
  84.  
  85.     /* build the complete an IEEE 754 single precision value */
  86.     z = (sz << 31) | ((ez & 0x000000ffUL) << 23) | (fz & 0x007fffffUL);
  87.  
  88.   } else {
  89.  
  90.     z = 0UL;
  91.  
  92.   }
  93.  
  94.   /* set the value */
  95.   *zp = *(float*)&z;
  96. }
  97.  
「IEEE_754単精度の足し算を行うルーチン」をウィキ内検索
最終更新：2008年03月11日 23:27
ツールボックス

下から選んでください:

新しいページを作成する
ヘルプ / FAQ もご覧ください。