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