Subversion Repositories HelenOS

Rev

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