Subversion Repositories HelenOS-historic

Rev

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