Subversion Repositories HelenOS-historic

Rev

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