Subversion Repositories HelenOS-historic

Rev

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.  /** @addtogroup softfloat 
  30.  * @{
  31.  */
  32. /** @file
  33.  */
  34.  
  35. #include<sftypes.h>
  36. #include<add.h>
  37. #include<div.h>
  38. #include<comparison.h>
  39. #include<mul.h>
  40. #include<common.h>
  41.  
  42.  
  43. float32 divFloat32(float32 a, float32 b)
  44. {
  45.     float32 result;
  46.     int32_t aexp, bexp, cexp;
  47.     uint64_t afrac, bfrac, cfrac;
  48.    
  49.     result.parts.sign = a.parts.sign ^ b.parts.sign;
  50.    
  51.     if (isFloat32NaN(a)) {
  52.         if (isFloat32SigNaN(a)) {
  53.             /*FIXME: SigNaN*/
  54.         }
  55.         /*NaN*/
  56.         return a;
  57.     }
  58.    
  59.     if (isFloat32NaN(b)) {
  60.         if (isFloat32SigNaN(b)) {
  61.             /*FIXME: SigNaN*/
  62.         }
  63.         /*NaN*/
  64.         return b;
  65.     }
  66.    
  67.     if (isFloat32Infinity(a)) {
  68.         if (isFloat32Infinity(b)) {
  69.             /*FIXME: inf / inf */
  70.             result.binary = FLOAT32_NAN;
  71.             return result;
  72.         }
  73.         /* inf / num */
  74.         result.parts.exp = a.parts.exp;
  75.         result.parts.fraction = a.parts.fraction;
  76.         return result;
  77.     }
  78.  
  79.     if (isFloat32Infinity(b)) {
  80.         if (isFloat32Zero(a)) {
  81.             /* FIXME 0 / inf */
  82.             result.parts.exp = 0;
  83.             result.parts.fraction = 0;
  84.             return result;
  85.         }
  86.         /* FIXME: num / inf*/
  87.         result.parts.exp = 0;
  88.         result.parts.fraction = 0;
  89.         return result;
  90.     }
  91.    
  92.     if (isFloat32Zero(b)) {
  93.         if (isFloat32Zero(a)) {
  94.             /*FIXME: 0 / 0*/
  95.             result.binary = FLOAT32_NAN;
  96.             return result;
  97.         }
  98.         /* FIXME: division by zero */
  99.         result.parts.exp = 0;
  100.         result.parts.fraction = 0;
  101.         return result;
  102.     }
  103.  
  104.    
  105.     afrac = a.parts.fraction;
  106.     aexp = a.parts.exp;
  107.     bfrac = b.parts.fraction;
  108.     bexp = b.parts.exp;
  109.    
  110.     /* denormalized numbers */
  111.     if (aexp == 0) {
  112.         if (afrac == 0) {
  113.         result.parts.exp = 0;
  114.         result.parts.fraction = 0;
  115.         return result;
  116.         }
  117.         /* normalize it*/
  118.        
  119.         afrac <<= 1;
  120.             /* afrac is nonzero => it must stop */ 
  121.         while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
  122.             afrac <<= 1;
  123.             aexp--;
  124.         }
  125.     }
  126.  
  127.     if (bexp == 0) {
  128.         bfrac <<= 1;
  129.             /* bfrac is nonzero => it must stop */ 
  130.         while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
  131.             bfrac <<= 1;
  132.             bexp--;
  133.         }
  134.     }
  135.  
  136.     afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
  137.     bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
  138.  
  139.     if ( bfrac <= (afrac << 1) ) {
  140.         afrac >>= 1;
  141.         aexp++;
  142.     }
  143.    
  144.     cexp = aexp - bexp + FLOAT32_BIAS - 2;
  145.    
  146.     cfrac = (afrac << 32) / bfrac;
  147.     if ((  cfrac & 0x3F ) == 0) {
  148.         cfrac |= ( bfrac * cfrac != afrac << 32 );
  149.     }
  150.    
  151.     /* pack and round */
  152.    
  153.     /* find first nonzero digit and shift result and detect possibly underflow */
  154.     while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
  155.         cexp--;
  156.         cfrac <<= 1;
  157.             /* TODO: fix underflow */
  158.     };
  159.    
  160.     cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
  161.    
  162.     if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
  163.         ++cexp;
  164.         cfrac >>= 1;
  165.         }  
  166.  
  167.     /* check overflow */
  168.     if (cexp >= FLOAT32_MAX_EXPONENT ) {
  169.         /* FIXME: overflow, return infinity */
  170.         result.parts.exp = FLOAT32_MAX_EXPONENT;
  171.         result.parts.fraction = 0;
  172.         return result;
  173.     }
  174.  
  175.     if (cexp < 0) {
  176.         /* FIXME: underflow */
  177.         result.parts.exp = 0;
  178.         if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
  179.             result.parts.fraction = 0;
  180.             return result;
  181.         }
  182.         cfrac >>= 1;
  183.         while (cexp < 0) {
  184.             cexp ++;
  185.             cfrac >>= 1;
  186.         }
  187.        
  188.     } else {
  189.         result.parts.exp = (uint32_t)cexp;
  190.     }
  191.    
  192.     result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
  193.    
  194.     return result; 
  195. }
  196.  
  197. float64 divFloat64(float64 a, float64 b)
  198. {
  199.     float64 result;
  200.     int64_t aexp, bexp, cexp;
  201.     uint64_t afrac, bfrac, cfrac;
  202.     uint64_t remlo, remhi;
  203.    
  204.     result.parts.sign = a.parts.sign ^ b.parts.sign;
  205.    
  206.     if (isFloat64NaN(a)) {
  207.        
  208.         if (isFloat64SigNaN(b)) {
  209.             /*FIXME: SigNaN*/
  210.             return b;
  211.         }
  212.        
  213.         if (isFloat64SigNaN(a)) {
  214.             /*FIXME: SigNaN*/
  215.         }
  216.         /*NaN*/
  217.         return a;
  218.     }
  219.    
  220.     if (isFloat64NaN(b)) {
  221.         if (isFloat64SigNaN(b)) {
  222.             /*FIXME: SigNaN*/
  223.         }
  224.         /*NaN*/
  225.         return b;
  226.     }
  227.    
  228.     if (isFloat64Infinity(a)) {
  229.         if (isFloat64Infinity(b) || isFloat64Zero(b)) {
  230.             /*FIXME: inf / inf */
  231.             result.binary = FLOAT64_NAN;
  232.             return result;
  233.         }
  234.         /* inf / num */
  235.         result.parts.exp = a.parts.exp;
  236.         result.parts.fraction = a.parts.fraction;
  237.         return result;
  238.     }
  239.  
  240.     if (isFloat64Infinity(b)) {
  241.         if (isFloat64Zero(a)) {
  242.             /* FIXME 0 / inf */
  243.             result.parts.exp = 0;
  244.             result.parts.fraction = 0;
  245.             return result;
  246.         }
  247.         /* FIXME: num / inf*/
  248.         result.parts.exp = 0;
  249.         result.parts.fraction = 0;
  250.         return result;
  251.     }
  252.    
  253.     if (isFloat64Zero(b)) {
  254.         if (isFloat64Zero(a)) {
  255.             /*FIXME: 0 / 0*/
  256.             result.binary = FLOAT64_NAN;
  257.             return result;
  258.         }
  259.         /* FIXME: division by zero */
  260.         result.parts.exp = 0;
  261.         result.parts.fraction = 0;
  262.         return result;
  263.     }
  264.  
  265.    
  266.     afrac = a.parts.fraction;
  267.     aexp = a.parts.exp;
  268.     bfrac = b.parts.fraction;
  269.     bexp = b.parts.exp;
  270.    
  271.     /* denormalized numbers */
  272.     if (aexp == 0) {
  273.         if (afrac == 0) {
  274.             result.parts.exp = 0;
  275.             result.parts.fraction = 0;
  276.             return result;
  277.         }
  278.         /* normalize it*/
  279.        
  280.         aexp++;
  281.             /* afrac is nonzero => it must stop */ 
  282.         while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
  283.             afrac <<= 1;
  284.             aexp--;
  285.         }
  286.     }
  287.  
  288.     if (bexp == 0) {
  289.         bexp++;
  290.             /* bfrac is nonzero => it must stop */ 
  291.         while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
  292.             bfrac <<= 1;
  293.             bexp--;
  294.         }
  295.     }
  296.  
  297.     afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
  298.     bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
  299.  
  300.     if ( bfrac <= (afrac << 1) ) {
  301.         afrac >>= 1;
  302.         aexp++;
  303.     }
  304.    
  305.     cexp = aexp - bexp + FLOAT64_BIAS - 2;
  306.    
  307.     cfrac = divFloat64estim(afrac, bfrac);
  308.    
  309.     if ((  cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
  310.         mul64integers( bfrac, cfrac, &remlo, &remhi);
  311.         /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ 
  312.         remhi = afrac - remhi - ( remlo > 0);
  313.         remlo = - remlo;
  314.        
  315.         while ((int64_t) remhi < 0) {
  316.             cfrac--;
  317.             remlo += bfrac;
  318.             remhi += ( remlo < bfrac );
  319.         }
  320.         cfrac |= ( remlo != 0 );
  321.     }
  322.    
  323.     /* round and shift */
  324.     result = finishFloat64(cexp, cfrac, result.parts.sign);
  325.     return result;
  326.  
  327. }
  328.  
  329. uint64_t divFloat64estim(uint64_t a, uint64_t b)
  330. {
  331.     uint64_t bhi;
  332.     uint64_t remhi, remlo;
  333.     uint64_t result;
  334.    
  335.     if ( b <= a ) {
  336.         return 0xFFFFFFFFFFFFFFFFull;
  337.     }
  338.    
  339.     bhi = b >> 32;
  340.     result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
  341.     mul64integers(b, result, &remlo, &remhi);
  342.    
  343.     remhi = a - remhi - (remlo > 0);
  344.     remlo = - remlo;
  345.  
  346.     b <<= 32;
  347.     while ( (int64_t) remhi < 0 ) {
  348.             result -= 0x1ll << 32; 
  349.             remlo += b;
  350.             remhi += bhi + ( remlo < b );
  351.         }
  352.     remhi = (remhi << 32) | (remlo >> 32);
  353.     if (( bhi << 32) <= remhi) {
  354.         result |= 0xFFFFFFFF;
  355.     } else {
  356.         result |= remhi / bhi;
  357.     }
  358.    
  359.    
  360.     return result;
  361. }
  362.  
  363.  
  364.  /** @}
  365.  */
  366.  
  367.