Subversion Repositories HelenOS

Rev

Rev 828 | Rev 835 | 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.     __s32 aexp, bexp, cexp;
  41.     __u64 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 = (__u32)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.     __s32 aexp, bexp, cexp;
  195.     __u64 afrac, bfrac, cfrac;
  196.     __u64 remlo, remhi;
  197.    
  198.     result.parts.sign = a.parts.sign ^ b.parts.sign;
  199.    
  200.     if (isFloat64NaN(a)) {
  201.         if (isFloat64SigNaN(a)) {
  202.             /*FIXME: SigNaN*/
  203.         }
  204.         /*NaN*/
  205.         return a;
  206.     }
  207.    
  208.     if (isFloat64NaN(b)) {
  209.         if (isFloat64SigNaN(b)) {
  210.             /*FIXME: SigNaN*/
  211.         }
  212.         /*NaN*/
  213.         return b;
  214.     }
  215.    
  216.     if (isFloat64Infinity(a)) {
  217.         if (isFloat64Infinity(b)) {
  218.             /*FIXME: inf / inf */
  219.             result.binary = FLOAT64_NAN;
  220.             return result;
  221.         }
  222.         /* inf / num */
  223.         result.parts.exp = a.parts.exp;
  224.         result.parts.fraction = a.parts.fraction;
  225.         return result;
  226.     }
  227.  
  228.     if (isFloat64Infinity(b)) {
  229.         if (isFloat64Zero(a)) {
  230.             /* FIXME 0 / inf */
  231.             result.parts.exp = 0;
  232.             result.parts.fraction = 0;
  233.             return result;
  234.         }
  235.         /* FIXME: num / inf*/
  236.         result.parts.exp = 0;
  237.         result.parts.fraction = 0;
  238.         return result;
  239.     }
  240.    
  241.     if (isFloat64Zero(b)) {
  242.         if (isFloat64Zero(a)) {
  243.             /*FIXME: 0 / 0*/
  244.             result.binary = FLOAT64_NAN;
  245.             return result;
  246.         }
  247.         /* FIXME: division by zero */
  248.         result.parts.exp = 0;
  249.         result.parts.fraction = 0;
  250.         return result;
  251.     }
  252.  
  253.    
  254.     afrac = a.parts.fraction;
  255.     aexp = a.parts.exp;
  256.     bfrac = b.parts.fraction;
  257.     bexp = b.parts.exp;
  258.    
  259.     /* denormalized numbers */
  260.     if (aexp == 0) {
  261.         if (afrac == 0) {
  262.         result.parts.exp = 0;
  263.         result.parts.fraction = 0;
  264.         return result;
  265.         }
  266.         /* normalize it*/
  267.        
  268.         afrac <<= 1;
  269.             /* afrac is nonzero => it must stop */ 
  270.         while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
  271.             afrac <<= 1;
  272.             aexp--;
  273.         }
  274.     }
  275.  
  276.     if (bexp == 0) {
  277.         bfrac <<= 1;
  278.             /* bfrac is nonzero => it must stop */ 
  279.         while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
  280.             bfrac <<= 1;
  281.             bexp--;
  282.         }
  283.     }
  284.  
  285.     afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
  286.     bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
  287.  
  288.     if ( bfrac <= (afrac << 1) ) {
  289.         afrac >>= 1;
  290.         aexp++;
  291.     }
  292.    
  293.     cexp = aexp - bexp + FLOAT64_BIAS - 2;
  294.    
  295.     cfrac = divFloat64estim(afrac, bfrac);
  296.    
  297.     if ((  cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
  298.         mul64integers( bfrac, cfrac, &remlo, &remhi);
  299.         /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ 
  300.         remhi = afrac - remhi - ( remlo > 0);
  301.         remlo = - remlo;
  302.        
  303.         while ((__s64) remhi < 0) {
  304.             cfrac--;
  305.             remlo += bfrac;
  306.             remhi += ( remlo < bfrac );
  307.         }
  308.         cfrac |= ( remlo != 0 );
  309.     }
  310.    
  311.     /* round and shift */
  312.     result = finishFloat64(cexp, cfrac, result.parts.sign);
  313.     return result;
  314.  
  315. }
  316.  
  317. __u64 divFloat64estim(__u64 a, __u64 b)
  318. {
  319.     __u64 bhi;
  320.     __u64 remhi, remlo;
  321.     __u64 result;
  322.    
  323.     if ( b <= a ) {
  324.         return 0xFFFFFFFFFFFFFFFFull;
  325.     }
  326.    
  327.     bhi = b >> 32;
  328.     result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
  329.     mul64integers(b, result, &remlo, &remhi);
  330.    
  331.     remhi = a - remhi - (remlo > 0);
  332.     remlo = - remlo;
  333.  
  334.     b <<= 32;
  335.     while ( (__s64) remhi < 0 ) {
  336.             result -= 0x1ll << 32; 
  337.             remlo += b;
  338.             remhi += bhi + ( remlo < b );
  339.         }
  340.     remhi = (remhi << 32) | (remlo >> 32);
  341.     if (( bhi << 32) <= remhi) {
  342.         result |= 0xFFFFFFFF;
  343.     } else {
  344.         result |= remhi / bhi;
  345.     }
  346.    
  347.    
  348.     return result;
  349. }
  350.  
  351.