Subversion Repositories HelenOS

Rev

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