Subversion Repositories HelenOS-historic

Rev

Rev 661 | Go to most recent revision | Blame | 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<arithmetic.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,mant1,mant2;
  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==0xFF) {
  51.             return b;
  52.         }
  53.        
  54.         mant1=b.parts.mantisa;
  55.         exp1=b.parts.exp;
  56.         mant2=a.parts.mantisa;
  57.         exp2=a.parts.exp;
  58.         expdiff*=-1;
  59.     } else {
  60.         if (isFloat32NaN(a)) {
  61.             //TODO: fix SigNaN
  62.             if ((isFloat32SigNaN(a))||(isFloat32SigNaN(b))) {
  63.             };
  64.             return a;
  65.         };
  66.        
  67.         if (a.parts.exp==0xFF) {
  68.             return a;
  69.         }
  70.        
  71.         mant1=a.parts.mantisa;
  72.         exp1=a.parts.exp;
  73.         mant2=b.parts.mantisa;
  74.         exp2=b.parts.exp;
  75.     };
  76.    
  77.     if (exp1==0) {
  78.         //both are denormalized
  79.         mant1+=mant2;
  80.         if (mant1&0xF00000) {
  81.             a.parts.exp=1;
  82.         };
  83.         a.parts.mantisa=mant1;
  84.         return a;
  85.     };
  86.    
  87.     // create some space for rounding
  88.     mant1<<=6;
  89.     mant2<<=6;
  90.    
  91.     mant1|=0x20000000; //add hidden bit
  92.    
  93.    
  94.     if (exp2==0) {
  95.         --expdiff; 
  96.     } else {
  97.         mant2|=0x20000000; //hidden bit
  98.     };
  99.    
  100.     if (expdiff>24) {
  101.          goto done;
  102.          };
  103.    
  104.     mant2>>=expdiff;
  105.     mant1+=mant2;
  106. done:
  107.     if (mant1&0x40000000) {
  108.         ++exp1;
  109.         mant1>>=1;
  110.     };
  111.    
  112.     //rounding - if first bit after mantisa is set then round up
  113.     mant1+=0x20;
  114.    
  115.     if (mant1&0x40000000) {
  116.         ++exp1;
  117.         mant1>>=1;
  118.     };
  119.    
  120.     a.parts.exp=exp1;
  121.     a.parts.mantisa = ((mant1&(~0x20000000))>>6); /*Clear hidden bit and shift */
  122.     return a;
  123. };
  124.  
  125. /** Subtract two float32 numbers with same signs
  126.  */
  127. float32 subFloat32(float32 a, float32 b)
  128. {
  129.     int expdiff;
  130.     __u32 exp1,exp2,mant1,mant2;
  131.     float32 result;
  132.  
  133.     result.f = 0;
  134.    
  135.     expdiff=a.parts.exp - b.parts.exp;
  136.     if ((expdiff<0)||((expdiff==0)&&(a.parts.mantisa<b.parts.mantisa))) {
  137.         if (isFloat32NaN(b)) {
  138.             //TODO: fix SigNaN
  139.             if (isFloat32SigNaN(b)) {
  140.             };
  141.             return b;
  142.         };
  143.        
  144.         if (b.parts.exp==0xFF) {
  145.             b.parts.sign=!b.parts.sign; /* num -(+-inf) = -+inf */
  146.             return b;
  147.         }
  148.        
  149.         result.parts.sign = !a.parts.sign;
  150.        
  151.         mant1=b.parts.mantisa;
  152.         exp1=b.parts.exp;
  153.         mant2=a.parts.mantisa;
  154.         exp2=a.parts.exp;
  155.         expdiff*=-1;
  156.     } else {
  157.         if (isFloat32NaN(a)) {
  158.             //TODO: fix SigNaN
  159.             if ((isFloat32SigNaN(a))||(isFloat32SigNaN(b))) {
  160.             };
  161.             return a;
  162.         };
  163.        
  164.         if (a.parts.exp==0xFF) {
  165.             if (b.parts.exp==0xFF) {
  166.                 /* inf - inf => nan */
  167.                 //TODO: fix exception
  168.                 result.binary = FLOAT32_NAN;
  169.                 return result;
  170.             };
  171.             return a;
  172.         }
  173.        
  174.         result.parts.sign = a.parts.sign;
  175.        
  176.         mant1=a.parts.mantisa;
  177.         exp1=a.parts.exp;
  178.         mant2=b.parts.mantisa;
  179.         exp2=b.parts.exp;
  180.        
  181.  
  182.        
  183.     };
  184.    
  185.     if (exp1==0) {
  186.         //both are denormalized
  187.         result.parts.mantisa=mant1-mant2;
  188.         if (result.parts.mantisa>mant1) {
  189.             //TODO: underflow exception
  190.             return result;
  191.         };
  192.         result.parts.exp=0;
  193.         return result;
  194.     };
  195.    
  196.     // create some space for rounding
  197.     mant1<<=6;
  198.     mant2<<=6;
  199.    
  200.     mant1|=0x20000000; //add hidden bit
  201.    
  202.    
  203.     if (exp2==0) {
  204.         --expdiff; 
  205.     } else {
  206.         mant2|=0x20000000; //hidden bit
  207.     };
  208.    
  209.     if (expdiff>24) {
  210.          goto done;
  211.          };
  212.    
  213.     mant1 = mant1-(mant2>>expdiff);
  214. done:
  215.    
  216.     //TODO: find first nonzero digit and shift result and detect possibly underflow
  217.     while ((exp1>0)&&(!(mant1&0x20000000))) {
  218.         exp1--;
  219.         mant1 <<= 1;
  220.         if(mant1 == 0) {
  221.             /* Realy is it an underflow? ... */
  222.             /* TODO: fix underflow */
  223.         };
  224.     };
  225.    
  226.     //rounding - if first bit after mantisa is set then round up   
  227.     mant1 += 0x20;
  228.  
  229.     if (mant1&0x40000000) {
  230.         ++exp1;
  231.         mant1>>=1;
  232.     };
  233.    
  234.     result.parts.mantisa = ((mant1&(~0x20000000))>>6); /*Clear hidden bit and shift */
  235.     result.parts.exp = exp1;
  236.    
  237.     return result;
  238. };
  239.  
  240. /** Multiply two 32 bit float numbers
  241.  *
  242.  */
  243. float32 mulFloat32(float32 a, float32 b)
  244. {
  245.     float32 result;
  246.     __u64 mant1, mant2;
  247.     __s32 exp;
  248.  
  249.     result.parts.sign = a.parts.sign ^ b.parts.sign;
  250.    
  251.     if ((isFloat32NaN(a))||(isFloat32NaN(b))) {
  252.         /* TODO: fix SigNaNs */
  253.         if (isFloat32SigNaN(a)) {
  254.             result.parts.mantisa = a.parts.mantisa;
  255.             result.parts.exp = a.parts.exp;
  256.             return result;
  257.         };
  258.         if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
  259.             result.parts.mantisa = b.parts.mantisa;
  260.             result.parts.exp = b.parts.exp;
  261.             return result;
  262.         };
  263.         /* set NaN as result */
  264.         result.parts.mantisa = 0x1;
  265.         result.parts.exp = 0xFF;
  266.         return result;
  267.     };
  268.        
  269.     if (isFloat32Infinity(a)) {
  270.         if (isFloat32Zero(b)) {
  271.             /* FIXME: zero * infinity */
  272.             result.parts.mantisa = 0x1;
  273.             result.parts.exp = 0xFF;
  274.             return result;
  275.         }
  276.         result.parts.mantisa = a.parts.mantisa;
  277.         result.parts.exp = a.parts.exp;
  278.         return result;
  279.     }
  280.  
  281.     if (isFloat32Infinity(b)) {
  282.         if (isFloat32Zero(a)) {
  283.             /* FIXME: zero * infinity */
  284.             result.parts.mantisa = 0x1;
  285.             result.parts.exp = 0xFF;
  286.             return result;
  287.         }
  288.         result.parts.mantisa = b.parts.mantisa;
  289.         result.parts.exp = b.parts.exp;
  290.         return result;
  291.     }
  292.  
  293.     /* exp is signed so we can easy detect underflow */
  294.     exp = a.parts.exp + b.parts.exp;
  295.     exp -= FLOAT32_BIAS;
  296.    
  297.     if (exp >= 0xFF ) {
  298.         /* FIXME: overflow */
  299.         /* set infinity as result */
  300.         result.parts.mantisa = 0x0;
  301.         result.parts.exp = 0xFF;
  302.         return result;
  303.     };
  304.    
  305.     if (exp < 0) {
  306.         /* FIXME: underflow */
  307.         /* return signed zero */
  308.         result.parts.mantisa = 0x0;
  309.         result.parts.exp = 0x0;
  310.         return result;
  311.     };
  312.    
  313.     mant1 = a.parts.mantisa;
  314.     if (a.parts.exp>0) {
  315.         mant1 |= 0x800000;
  316.     } else {
  317.         ++exp;
  318.     };
  319.    
  320.     mant2 = b.parts.mantisa;
  321.     if (b.parts.exp>0) {
  322.         mant2 |= 0x800000;
  323.     } else {
  324.         ++exp;
  325.     };
  326.  
  327.     mant1 <<= 1; /* one bit space for rounding */
  328.  
  329.     mant1 = mant1 * mant2;
  330. /* round and return */
  331.    
  332.     while ((exp < 0xFF )&&(mant1 > 0x1FFFFFF )) { /* 0xFFFFFF is 23 bits of mantisa + one more for hidden bit (all shifted 1 bit left)*/
  333.         ++exp;
  334.         mant1 >>= 1;
  335.     };
  336.  
  337.     /* rounding */
  338.     //++mant1; /* FIXME: not works - without it is ok */
  339.     mant1 >>= 1; /* shift off rounding space */
  340.    
  341.     if ((exp < 0xFF )&&(mant1 > 0xFFFFFF )) {
  342.         ++exp;
  343.         mant1 >>= 1;
  344.     };
  345.  
  346.     if (exp >= 0xFF ) {
  347.         /* TODO: fix overflow */
  348.         /* return infinity*/
  349.         result.parts.exp = 0xFF;
  350.         result.parts.mantisa = 0x0;
  351.         return result;
  352.     }
  353.    
  354.     exp -= FLOAT32_MANTISA_SIZE;
  355.  
  356.     if (exp <= FLOAT32_MANTISA_SIZE) {
  357.         /* denormalized number */
  358.         mant1 >>= 1; /* denormalize */
  359.         while ((mant1 > 0) && (exp < 0)) {
  360.             mant1 >>= 1;
  361.             ++exp;
  362.         };
  363.         if (mant1 == 0) {
  364.             /* FIXME : underflow */
  365.         result.parts.exp = 0;
  366.         result.parts.mantisa = 0;
  367.         return result;
  368.         };
  369.     };
  370.     result.parts.exp = exp;
  371.     result.parts.mantisa = mant1 & 0x7FFFFF;
  372.    
  373.     return result; 
  374.    
  375. };
  376.  
  377.  
  378.