Back to index

glibc  2.9
s_remquof.c
Go to the documentation of this file.
00001 /* Compute remainder and a congruent to the quotient.
00002    Copyright (C) 1997 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
00005 
00006    The GNU C Library is free software; you can redistribute it and/or
00007    modify it under the terms of the GNU Lesser General Public
00008    License as published by the Free Software Foundation; either
00009    version 2.1 of the License, or (at your option) any later version.
00010 
00011    The GNU C Library is distributed in the hope that it will be useful,
00012    but WITHOUT ANY WARRANTY; without even the implied warranty of
00013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014    Lesser General Public License for more details.
00015 
00016    You should have received a copy of the GNU Lesser General Public
00017    License along with the GNU C Library; if not, write to the Free
00018    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00019    02111-1307 USA.  */
00020 
00021 #include <math.h>
00022 
00023 #include "math_private.h"
00024 
00025 
00026 static const float zero = 0.0;
00027 
00028 
00029 float
00030 __remquof (float x, float y, int *quo)
00031 {
00032   int32_t hx,hy;
00033   u_int32_t sx;
00034   int cquo, qs;
00035 
00036   GET_FLOAT_WORD (hx, x);
00037   GET_FLOAT_WORD (hy, y);
00038   sx = hx & 0x80000000;
00039   qs = sx ^ (hy & 0x80000000);
00040   hy &= 0x7fffffff;
00041   hx &= 0x7fffffff;
00042 
00043   /* Purge off exception values.  */
00044   if (hy == 0)
00045     return (x * y) / (x * y);                    /* y = 0 */
00046   if ((hx >= 0x7f800000)                  /* x not finite */
00047       || (hy > 0x7f800000))               /* y is NaN */
00048     return (x * y) / (x * y);
00049 
00050   if (hy <= 0x7dffffff)
00051     x = __ieee754_fmodf (x, 8 * y);              /* now x < 8y */
00052 
00053   if ((hx - hy) == 0)
00054     {
00055       *quo = qs ? -1 : 1;
00056       return zero * x;
00057     }
00058 
00059   x  = fabsf (x);
00060   y  = fabsf (y);
00061   cquo = 0;
00062 
00063   if (x >= 4 * y)
00064     {
00065       x -= 4 * y;
00066       cquo += 4;
00067     }
00068   if (x >= 2 * y)
00069     {
00070       x -= 2 * y;
00071       cquo += 2;
00072     }
00073 
00074   if (hy < 0x01000000)
00075     {
00076       if (x + x > y)
00077        {
00078          x -= y;
00079          ++cquo;
00080          if (x + x >= y)
00081            {
00082              x -= y;
00083              ++cquo;
00084            }
00085        }
00086     }
00087   else
00088     {
00089       float y_half = 0.5 * y;
00090       if (x > y_half)
00091        {
00092          x -= y;
00093          ++cquo;
00094          if (x >= y_half)
00095            {
00096              x -= y;
00097              ++cquo;
00098            }
00099        }
00100     }
00101 
00102   *quo = qs ? -cquo : cquo;
00103 
00104   if (sx)
00105     x = -x;
00106   return x;
00107 }
00108 weak_alias (__remquof, remquof)