Subversion Repositories HelenOS-historic

Rev

Rev 734 | Rev 828 | 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.         };
  106.    
  107.     if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
  108.         ++exp1;
  109.         frac1 >>= 1;
  110.     };
  111.    
  112.     /* rounding - if first bit after fraction is set then round up */
  113.     frac1 += (0x1 << 5);
  114.    
  115.     if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
  116.         /* rounding overflow */
  117.         ++exp1;
  118.         frac1 >>= 1;
  119.     };
  120.    
  121.     if ((a.parts.exp == FLOAT32_MAX_EXPONENT ) || (a.parts.exp < exp1)) {
  122.             /* overflow - set infinity as result */
  123.             a.parts.exp = FLOAT32_MAX_EXPONENT;
  124.             a.parts.fraction = 0;
  125.             return a;
  126.             }
  127.    
  128.     a.parts.exp = exp1;
  129.    
  130.     /*Clear hidden bit and shift */
  131.     a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
  132.     return a;
  133. }
  134.  
  135.  
  136. /** Add two Float64 numbers with same signs
  137.  */
  138. float64 addFloat64(float64 a, float64 b)
  139. {
  140.     int expdiff;
  141.     __u32 exp1, exp2;
  142.     __u64 frac1, frac2;
  143.    
  144.     expdiff = a.parts.exp - b.parts.exp;
  145.     if (expdiff < 0) {
  146.         if (isFloat64NaN(b)) {
  147.             /* TODO: fix SigNaN */
  148.             if (isFloat64SigNaN(b)) {
  149.             };
  150.  
  151.             return b;
  152.         };
  153.        
  154.         /* b is infinity and a not */  
  155.         if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
  156.             return b;
  157.         }
  158.        
  159.         frac1 = b.parts.fraction;
  160.         exp1 = b.parts.exp;
  161.         frac2 = a.parts.fraction;
  162.         exp2 = a.parts.exp;
  163.         expdiff *= -1;
  164.     } else {
  165.         if (isFloat64NaN(a)) {
  166.             /* TODO: fix SigNaN */
  167.             if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
  168.             };
  169.             return a;
  170.         };
  171.        
  172.         /* a is infinity and b not */
  173.         if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
  174.             return a;
  175.         }
  176.        
  177.         frac1 = a.parts.fraction;
  178.         exp1 = a.parts.exp;
  179.         frac2 = b.parts.fraction;
  180.         exp2 = b.parts.exp;
  181.     };
  182.    
  183.     if (exp1 == 0) {
  184.         /* both are denormalized */
  185.         frac1 += frac2;
  186.         if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
  187.             /* result is not denormalized */
  188.             a.parts.exp = 1;
  189.         };
  190.         a.parts.fraction = frac1;
  191.         return a;
  192.     };
  193.    
  194.     /* add hidden bit - frac1 is sure not denormalized */
  195.     frac1 |= FLOAT64_HIDDEN_BIT_MASK;
  196.  
  197.     /* second operand ... */
  198.     if (exp2 == 0) {
  199.         /* ... is denormalized */
  200.         --expdiff; 
  201.     } else {
  202.         /* is not denormalized */
  203.         frac2 |= FLOAT64_HIDDEN_BIT_MASK;
  204.     };
  205.    
  206.     /* create some space for rounding */
  207.     frac1 <<= 6;
  208.     frac2 <<= 6;
  209.    
  210.     if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
  211.         frac2 >>= expdiff;
  212.         frac1 += frac2;
  213.         };
  214.    
  215.     if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
  216.         ++exp1;
  217.         frac1 >>= 1;
  218.     };
  219.    
  220.     /* rounding - if first bit after fraction is set then round up */
  221.     frac1 += (0x1 << 5);
  222.    
  223.     if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
  224.         /* rounding overflow */
  225.         ++exp1;
  226.         frac1 >>= 1;
  227.     };
  228.    
  229.     if ((a.parts.exp == FLOAT64_MAX_EXPONENT ) || (a.parts.exp < exp1)) {
  230.             /* overflow - set infinity as result */
  231.             a.parts.exp = FLOAT64_MAX_EXPONENT;
  232.             a.parts.fraction = 0;
  233.             return a;
  234.             }
  235.    
  236.     a.parts.exp = exp1;
  237.     /*Clear hidden bit and shift */
  238.     a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
  239.     return a;
  240. }
  241.  
  242.  
  243.