Subversion Repositories HelenOS-historic

Rev

Rev 1031 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed

  1. /*
  2.  * Copyright (C) 2005 Josef Cejka
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  *
  9.  * - Redistributions of source code must retain the above copyright
  10.  *   notice, this list of conditions and the following disclaimer.
  11.  * - Redistributions in binary form must reproduce the above copyright
  12.  *   notice, this list of conditions and the following disclaimer in the
  13.  *   documentation and/or other materials provided with the distribution.
  14.  * - The name of the author may not be used to endorse or promote products
  15.  *   derived from this software without specific prior written permission.
  16.  *
  17.  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  18.  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  19.  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  20.  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  21.  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  22.  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  23.  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  24.  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  25.  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  26.  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  27.  */
  28.  
  29.  /** @addtogroup softfloat 
  30.  * @{
  31.  */
  32. /** @file
  33.  */
  34.  
  35. #include<sftypes.h>
  36. #include<sub.h>
  37. #include<comparison.h>
  38.  
  39. /** Subtract two float32 numbers with same signs
  40.  */
  41. float32 subFloat32(float32 a, float32 b)
  42. {
  43.     int expdiff;
  44.     uint32_t exp1, exp2, frac1, frac2;
  45.     float32 result;
  46.  
  47.     result.f = 0;
  48.    
  49.     expdiff = a.parts.exp - b.parts.exp;
  50.     if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.fraction < b.parts.fraction))) {
  51.         if (isFloat32NaN(b)) {
  52.             /* TODO: fix SigNaN */
  53.             if (isFloat32SigNaN(b)) {
  54.             };
  55.             return b;
  56.         };
  57.        
  58.         if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
  59.             b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
  60.             return b;
  61.         }
  62.        
  63.         result.parts.sign = !a.parts.sign;
  64.        
  65.         frac1 = b.parts.fraction;
  66.         exp1 = b.parts.exp;
  67.         frac2 = a.parts.fraction;
  68.         exp2 = a.parts.exp;
  69.         expdiff *= -1;
  70.     } else {
  71.         if (isFloat32NaN(a)) {
  72.             /* TODO: fix SigNaN */
  73.             if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
  74.             };
  75.             return a;
  76.         };
  77.        
  78.         if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
  79.             if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
  80.                 /* inf - inf => nan */
  81.                 /* TODO: fix exception */
  82.                 result.binary = FLOAT32_NAN;
  83.                 return result;
  84.             };
  85.             return a;
  86.         }
  87.        
  88.         result.parts.sign = a.parts.sign;
  89.        
  90.         frac1 = a.parts.fraction;
  91.         exp1 = a.parts.exp;
  92.         frac2 = b.parts.fraction;
  93.         exp2 = b.parts.exp;
  94.     };
  95.    
  96.     if (exp1 == 0) {
  97.         /* both are denormalized */
  98.         result.parts.fraction = frac1-frac2;
  99.         if (result.parts.fraction > frac1) {
  100.             /* TODO: underflow exception */
  101.             return result;
  102.         };
  103.         result.parts.exp = 0;
  104.         return result;
  105.     };
  106.  
  107.     /* add hidden bit */
  108.     frac1 |= FLOAT32_HIDDEN_BIT_MASK;
  109.    
  110.     if (exp2 == 0) {
  111.         /* denormalized */
  112.         --expdiff; 
  113.     } else {
  114.         /* normalized */
  115.         frac2 |= FLOAT32_HIDDEN_BIT_MASK;
  116.     };
  117.    
  118.     /* create some space for rounding */
  119.     frac1 <<= 6;
  120.     frac2 <<= 6;
  121.    
  122.     if (expdiff > FLOAT32_FRACTION_SIZE + 1) {
  123.          goto done;
  124.          };
  125.    
  126.     frac1 = frac1 - (frac2 >> expdiff);
  127. done:
  128.     /* TODO: find first nonzero digit and shift result and detect possibly underflow */
  129.     while ((exp1 > 0) && (!(frac1 & (FLOAT32_HIDDEN_BIT_MASK << 6 )))) {
  130.         --exp1;
  131.         frac1 <<= 1;
  132.             /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
  133.     };
  134.    
  135.     /* rounding - if first bit after fraction is set then round up */
  136.     frac1 += 0x20;
  137.  
  138.     if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
  139.         ++exp1;
  140.         frac1 >>= 1;
  141.     };
  142.    
  143.     /*Clear hidden bit and shift */
  144.     result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
  145.     result.parts.exp = exp1;
  146.    
  147.     return result;
  148. }
  149.  
  150. /** Subtract two float64 numbers with same signs
  151.  */
  152. float64 subFloat64(float64 a, float64 b)
  153. {
  154.     int expdiff;
  155.     uint32_t exp1, exp2;
  156.     uint64_t frac1, frac2;
  157.     float64 result;
  158.  
  159.     result.d = 0;
  160.    
  161.     expdiff = a.parts.exp - b.parts.exp;
  162.     if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.fraction < b.parts.fraction))) {
  163.         if (isFloat64NaN(b)) {
  164.             /* TODO: fix SigNaN */
  165.             if (isFloat64SigNaN(b)) {
  166.             };
  167.             return b;
  168.         };
  169.        
  170.         if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
  171.             b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
  172.             return b;
  173.         }
  174.        
  175.         result.parts.sign = !a.parts.sign;
  176.        
  177.         frac1 = b.parts.fraction;
  178.         exp1 = b.parts.exp;
  179.         frac2 = a.parts.fraction;
  180.         exp2 = a.parts.exp;
  181.         expdiff *= -1;
  182.     } else {
  183.         if (isFloat64NaN(a)) {
  184.             /* TODO: fix SigNaN */
  185.             if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
  186.             };
  187.             return a;
  188.         };
  189.        
  190.         if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
  191.             if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
  192.                 /* inf - inf => nan */
  193.                 /* TODO: fix exception */
  194.                 result.binary = FLOAT64_NAN;
  195.                 return result;
  196.             };
  197.             return a;
  198.         }
  199.        
  200.         result.parts.sign = a.parts.sign;
  201.        
  202.         frac1 = a.parts.fraction;
  203.         exp1 = a.parts.exp;
  204.         frac2 = b.parts.fraction;
  205.         exp2 = b.parts.exp;
  206.     };
  207.    
  208.     if (exp1 == 0) {
  209.         /* both are denormalized */
  210.         result.parts.fraction = frac1 - frac2;
  211.         if (result.parts.fraction > frac1) {
  212.             /* TODO: underflow exception */
  213.             return result;
  214.         };
  215.         result.parts.exp = 0;
  216.         return result;
  217.     };
  218.  
  219.     /* add hidden bit */
  220.     frac1 |= FLOAT64_HIDDEN_BIT_MASK;
  221.    
  222.     if (exp2 == 0) {
  223.         /* denormalized */
  224.         --expdiff; 
  225.     } else {
  226.         /* normalized */
  227.         frac2 |= FLOAT64_HIDDEN_BIT_MASK;
  228.     };
  229.    
  230.     /* create some space for rounding */
  231.     frac1 <<= 6;
  232.     frac2 <<= 6;
  233.    
  234.     if (expdiff > FLOAT64_FRACTION_SIZE + 1) {
  235.          goto done;
  236.          };
  237.    
  238.     frac1 = frac1 - (frac2 >> expdiff);
  239. done:
  240.     /* TODO: find first nonzero digit and shift result and detect possibly underflow */
  241.     while ((exp1 > 0) && (!(frac1 & (FLOAT64_HIDDEN_BIT_MASK << 6 )))) {
  242.         --exp1;
  243.         frac1 <<= 1;
  244.             /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
  245.     };
  246.    
  247.     /* rounding - if first bit after fraction is set then round up */
  248.     frac1 += 0x20;
  249.  
  250.     if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
  251.         ++exp1;
  252.         frac1 >>= 1;
  253.     };
  254.    
  255.     /*Clear hidden bit and shift */
  256.     result.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK));
  257.     result.parts.exp = exp1;
  258.    
  259.     return result;
  260. }
  261.  
  262.  
  263.  /** @}
  264.  */
  265.  
  266.