Subversion Repositories HelenOS

Rev

Rev 2071 | 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<mul.h>
  37. #include<comparison.h>
  38. #include<common.h>
  39.  
  40. /** Multiply two 32 bit float numbers
  41.  *
  42.  */
  43. float32 mulFloat32(float32 a, float32 b)
  44. {
  45.     float32 result;
  46.     uint64_t frac1, frac2;
  47.     int32_t exp;
  48.  
  49.     result.parts.sign = a.parts.sign ^ b.parts.sign;
  50.    
  51.     if (isFloat32NaN(a) || isFloat32NaN(b) ) {
  52.         /* TODO: fix SigNaNs */
  53.         if (isFloat32SigNaN(a)) {
  54.             result.parts.fraction = a.parts.fraction;
  55.             result.parts.exp = a.parts.exp;
  56.             return result;
  57.         };
  58.         if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
  59.             result.parts.fraction = b.parts.fraction;
  60.             result.parts.exp = b.parts.exp;
  61.             return result;
  62.         };
  63.         /* set NaN as result */
  64.         result.binary = FLOAT32_NAN;
  65.         return result;
  66.     };
  67.        
  68.     if (isFloat32Infinity(a)) {
  69.         if (isFloat32Zero(b)) {
  70.             /* FIXME: zero * infinity */
  71.             result.binary = FLOAT32_NAN;
  72.             return result;
  73.         }
  74.         result.parts.fraction = a.parts.fraction;
  75.         result.parts.exp = a.parts.exp;
  76.         return result;
  77.     }
  78.  
  79.     if (isFloat32Infinity(b)) {
  80.         if (isFloat32Zero(a)) {
  81.             /* FIXME: zero * infinity */
  82.             result.binary = FLOAT32_NAN;
  83.             return result;
  84.         }
  85.         result.parts.fraction = b.parts.fraction;
  86.         result.parts.exp = b.parts.exp;
  87.         return result;
  88.     }
  89.  
  90.     /* exp is signed so we can easy detect underflow */
  91.     exp = a.parts.exp + b.parts.exp;
  92.     exp -= FLOAT32_BIAS;
  93.    
  94.     if (exp >= FLOAT32_MAX_EXPONENT) {
  95.         /* FIXME: overflow */
  96.         /* set infinity as result */
  97.         result.binary = FLOAT32_INF;
  98.         result.parts.sign = a.parts.sign ^ b.parts.sign;
  99.         return result;
  100.     };
  101.    
  102.     if (exp < 0) {
  103.         /* FIXME: underflow */
  104.         /* return signed zero */
  105.         result.parts.fraction = 0x0;
  106.         result.parts.exp = 0x0;
  107.         return result;
  108.     };
  109.    
  110.     frac1 = a.parts.fraction;
  111.     if (a.parts.exp > 0) {
  112.         frac1 |= FLOAT32_HIDDEN_BIT_MASK;
  113.     } else {
  114.         ++exp;
  115.     };
  116.    
  117.     frac2 = b.parts.fraction;
  118.  
  119.     if (b.parts.exp > 0) {
  120.         frac2 |= FLOAT32_HIDDEN_BIT_MASK;
  121.     } else {
  122.         ++exp;
  123.     };
  124.  
  125.     frac1 <<= 1; /* one bit space for rounding */
  126.  
  127.     frac1 = frac1 * frac2;
  128. /* round and return */
  129.    
  130.     while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
  131.         /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
  132.         ++exp;
  133.         frac1 >>= 1;
  134.     };
  135.  
  136.     /* rounding */
  137.     /* ++frac1; FIXME: not works - without it is ok */
  138.     frac1 >>= 1; /* shift off rounding space */
  139.    
  140.     if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
  141.         ++exp;
  142.         frac1 >>= 1;
  143.     };
  144.  
  145.     if (exp >= FLOAT32_MAX_EXPONENT ) {
  146.         /* TODO: fix overflow */
  147.         /* return infinity*/
  148.         result.parts.exp = FLOAT32_MAX_EXPONENT;
  149.         result.parts.fraction = 0x0;
  150.         return result;
  151.     }
  152.    
  153.     exp -= FLOAT32_FRACTION_SIZE;
  154.  
  155.     if (exp <= FLOAT32_FRACTION_SIZE) {
  156.         /* denormalized number */
  157.         frac1 >>= 1; /* denormalize */
  158.         while ((frac1 > 0) && (exp < 0)) {
  159.             frac1 >>= 1;
  160.             ++exp;
  161.         };
  162.         if (frac1 == 0) {
  163.             /* FIXME : underflow */
  164.         result.parts.exp = 0;
  165.         result.parts.fraction = 0;
  166.         return result;
  167.         };
  168.     };
  169.     result.parts.exp = exp;
  170.     result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
  171.    
  172.     return result; 
  173.    
  174. }
  175.  
  176. /** Multiply two 64 bit float numbers
  177.  *
  178.  */
  179. float64 mulFloat64(float64 a, float64 b)
  180. {
  181.     float64 result;
  182.     uint64_t frac1, frac2;
  183.     int32_t exp;
  184.  
  185.     result.parts.sign = a.parts.sign ^ b.parts.sign;
  186.    
  187.     if (isFloat64NaN(a) || isFloat64NaN(b) ) {
  188.         /* TODO: fix SigNaNs */
  189.         if (isFloat64SigNaN(a)) {
  190.             result.parts.fraction = a.parts.fraction;
  191.             result.parts.exp = a.parts.exp;
  192.             return result;
  193.         };
  194.         if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
  195.             result.parts.fraction = b.parts.fraction;
  196.             result.parts.exp = b.parts.exp;
  197.             return result;
  198.         };
  199.         /* set NaN as result */
  200.         result.binary = FLOAT64_NAN;
  201.         return result;
  202.     };
  203.        
  204.     if (isFloat64Infinity(a)) {
  205.         if (isFloat64Zero(b)) {
  206.             /* FIXME: zero * infinity */
  207.             result.binary = FLOAT64_NAN;
  208.             return result;
  209.         }
  210.         result.parts.fraction = a.parts.fraction;
  211.         result.parts.exp = a.parts.exp;
  212.         return result;
  213.     }
  214.  
  215.     if (isFloat64Infinity(b)) {
  216.         if (isFloat64Zero(a)) {
  217.             /* FIXME: zero * infinity */
  218.             result.binary = FLOAT64_NAN;
  219.             return result;
  220.         }
  221.         result.parts.fraction = b.parts.fraction;
  222.         result.parts.exp = b.parts.exp;
  223.         return result;
  224.     }
  225.  
  226.     /* exp is signed so we can easy detect underflow */
  227.     exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
  228.    
  229.     frac1 = a.parts.fraction;
  230.  
  231.     if (a.parts.exp > 0) {
  232.         frac1 |= FLOAT64_HIDDEN_BIT_MASK;
  233.     } else {
  234.         ++exp;
  235.     };
  236.    
  237.     frac2 = b.parts.fraction;
  238.  
  239.     if (b.parts.exp > 0) {
  240.         frac2 |= FLOAT64_HIDDEN_BIT_MASK;
  241.     } else {
  242.         ++exp;
  243.     };
  244.  
  245.     frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
  246.     frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
  247.  
  248.     mul64integers(frac1, frac2, &frac1, &frac2);
  249.  
  250.     frac2 |= (frac1 != 0);
  251.     if (frac2 & (0x1ll << 62)) {
  252.         frac2 <<= 1;
  253.         exp--;
  254.     }
  255.  
  256.     result = finishFloat64(exp, frac2, result.parts.sign);
  257.     return result;
  258. }
  259.  
  260. /** Multiply two 64 bit numbers and return result in two parts
  261.  * @param a first operand
  262.  * @param b second operand
  263.  * @param lo lower part from result
  264.  * @param hi higher part of result
  265.  */
  266. void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi)
  267. {
  268.     uint64_t low, high, middle1, middle2;
  269.     uint32_t alow, blow;
  270.  
  271.     alow = a & 0xFFFFFFFF;
  272.     blow = b & 0xFFFFFFFF;
  273.    
  274.     a >>= 32;
  275.     b >>= 32;
  276.    
  277.     low = ((uint64_t)alow) * blow;
  278.     middle1 = a * blow;
  279.     middle2 = alow * b;
  280.     high = a * b;
  281.  
  282.     middle1 += middle2;
  283.     high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
  284.     middle1 <<= 32;
  285.     low += middle1;
  286.     high += (low < middle1);
  287.     *lo = low;
  288.     *hi = high;
  289.    
  290.     return;
  291. }
  292.  
  293. /** @}
  294.  */
  295.