Subversion Repositories HelenOS

Rev

Rev 828 | Rev 1031 | Go to most recent revision | 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. #include<sftypes.h>
  30. #include<add.h>
  31. #include<comparison.h>
  32.  
  33. /** Add two Float32 numbers with same signs
  34.  */
  35. float32 addFloat32(float32 a, float32 b)
  36. {
  37.     int expdiff;
  38.     __u32 exp1, exp2,frac1, frac2;
  39.    
  40.     expdiff = a.parts.exp - b.parts.exp;
  41.     if (expdiff < 0) {
  42.         if (isFloat32NaN(b)) {
  43.             /* TODO: fix SigNaN */
  44.             if (isFloat32SigNaN(b)) {
  45.             };
  46.  
  47.             return b;
  48.         };
  49.        
  50.         if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
  51.             return b;
  52.         }
  53.        
  54.         frac1 = b.parts.fraction;
  55.         exp1 = b.parts.exp;
  56.         frac2 = a.parts.fraction;
  57.         exp2 = a.parts.exp;
  58.         expdiff *= -1;
  59.     } else {
  60.         if ((isFloat32NaN(a)) || (isFloat32NaN(b))) {
  61.             /* TODO: fix SigNaN */
  62.             if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
  63.             };
  64.             return (isFloat32NaN(a)?a:b);
  65.         };
  66.        
  67.         if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
  68.             return a;
  69.         }
  70.        
  71.         frac1 = a.parts.fraction;
  72.         exp1 = a.parts.exp;
  73.         frac2 = b.parts.fraction;
  74.         exp2 = b.parts.exp;
  75.     };
  76.    
  77.     if (exp1 == 0) {
  78.         /* both are denormalized */
  79.         frac1 += frac2;
  80.         if (frac1 & FLOAT32_HIDDEN_BIT_MASK ) {
  81.             /* result is not denormalized */
  82.             a.parts.exp = 1;
  83.         };
  84.         a.parts.fraction = frac1;
  85.         return a;
  86.     };
  87.    
  88.     frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
  89.  
  90.     if (exp2 == 0) {
  91.         /* second operand is denormalized */
  92.         --expdiff;
  93.     } else {
  94.         /* add hidden bit to second operand */
  95.         frac2 |= FLOAT32_HIDDEN_BIT_MASK;
  96.     };
  97.    
  98.     /* create some space for rounding */
  99.     frac1 <<= 6;
  100.     frac2 <<= 6;
  101.    
  102.     if (expdiff < (FLOAT32_FRACTION_SIZE + 2) ) {
  103.         frac2 >>= expdiff;
  104.         frac1 += frac2;
  105.     } else {
  106.         a.parts.exp = exp1;
  107.         a.parts.fraction = (frac1 >> 6) & (~(FLOAT32_HIDDEN_BIT_MASK));
  108.         return a;
  109.     }
  110.    
  111.     if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
  112.         ++exp1;
  113.         frac1 >>= 1;
  114.     };
  115.    
  116.     /* rounding - if first bit after fraction is set then round up */
  117.     frac1 += (0x1 << 5);
  118.    
  119.     if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
  120.         /* rounding overflow */
  121.         ++exp1;
  122.         frac1 >>= 1;
  123.     };
  124.    
  125.    
  126.     if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
  127.             /* overflow - set infinity as result */
  128.             a.parts.exp = FLOAT32_MAX_EXPONENT;
  129.             a.parts.fraction = 0;
  130.             return a;
  131.             }
  132.    
  133.     a.parts.exp = exp1;
  134.    
  135.     /*Clear hidden bit and shift */
  136.     a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
  137.     return a;
  138. }
  139.  
  140.  
  141. /** Add two Float64 numbers with same signs
  142.  */
  143. float64 addFloat64(float64 a, float64 b)
  144. {
  145.     int expdiff;
  146.     __u32 exp1, exp2;
  147.     __u64 frac1, frac2;
  148.    
  149.     expdiff = ((int )a.parts.exp) - b.parts.exp;
  150.     if (expdiff < 0) {
  151.         if (isFloat64NaN(b)) {
  152.             /* TODO: fix SigNaN */
  153.             if (isFloat64SigNaN(b)) {
  154.             };
  155.  
  156.             return b;
  157.         };
  158.        
  159.         /* b is infinity and a not */  
  160.         if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
  161.             return b;
  162.         }
  163.        
  164.         frac1 = b.parts.fraction;
  165.         exp1 = b.parts.exp;
  166.         frac2 = a.parts.fraction;
  167.         exp2 = a.parts.exp;
  168.         expdiff *= -1;
  169.     } else {
  170.         if (isFloat64NaN(a)) {
  171.             /* TODO: fix SigNaN */
  172.             if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
  173.             };
  174.             return a;
  175.         };
  176.        
  177.         /* a is infinity and b not */
  178.         if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
  179.             return a;
  180.         }
  181.        
  182.         frac1 = a.parts.fraction;
  183.         exp1 = a.parts.exp;
  184.         frac2 = b.parts.fraction;
  185.         exp2 = b.parts.exp;
  186.     };
  187.    
  188.     if (exp1 == 0) {
  189.         /* both are denormalized */
  190.         frac1 += frac2;
  191.         if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
  192.             /* result is not denormalized */
  193.             a.parts.exp = 1;
  194.         };
  195.         a.parts.fraction = frac1;
  196.         return a;
  197.     };
  198.    
  199.     /* add hidden bit - frac1 is sure not denormalized */
  200.     frac1 |= FLOAT64_HIDDEN_BIT_MASK;
  201.  
  202.     /* second operand ... */
  203.     if (exp2 == 0) {
  204.         /* ... is denormalized */
  205.         --expdiff; 
  206.     } else {
  207.         /* is not denormalized */
  208.         frac2 |= FLOAT64_HIDDEN_BIT_MASK;
  209.     };
  210.    
  211.     /* create some space for rounding */
  212.     frac1 <<= 6;
  213.     frac2 <<= 6;
  214.    
  215.     if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
  216.         frac2 >>= expdiff;
  217.         frac1 += frac2;
  218.     } else {
  219.         a.parts.exp = exp1;
  220.         a.parts.fraction = (frac1 >> 6) & (~(FLOAT64_HIDDEN_BIT_MASK));
  221.         return a;
  222.     }
  223.    
  224.     if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
  225.         ++exp1;
  226.         frac1 >>= 1;
  227.     };
  228.    
  229.     /* rounding - if first bit after fraction is set then round up */
  230.     frac1 += (0x1 << 5);
  231.    
  232.     if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
  233.         /* rounding overflow */
  234.         ++exp1;
  235.         frac1 >>= 1;
  236.     };
  237.    
  238.     if ((exp1 == FLOAT64_MAX_EXPONENT ) || (exp2 > exp1)) {
  239.             /* overflow - set infinity as result */
  240.             a.parts.exp = FLOAT64_MAX_EXPONENT;
  241.             a.parts.fraction = 0;
  242.             return a;
  243.             }
  244.    
  245.     a.parts.exp = exp1;
  246.     /*Clear hidden bit and shift */
  247.     a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
  248.    
  249.     return a;
  250. }
  251.  
  252.  
  253.