Back to index

glibc  2.9
s_remquo.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 double zero = 0.0;
00027 
00028 
00029 double
00030 __remquo (double x, double y, int *quo)
00031 {
00032   int32_t hx,hy;
00033   u_int32_t sx,lx,ly;
00034   int cquo, qs;
00035 
00036   EXTRACT_WORDS (hx, lx, x);
00037   EXTRACT_WORDS (hy, ly, 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 | ly) == 0)
00045     return (x * y) / (x * y);                    /* y = 0 */
00046   if ((hx >= 0x7ff00000)                  /* x not finite */
00047       || ((hy >= 0x7ff00000)                     /* p is NaN */
00048          && (((hy - 0x7ff00000) | ly) != 0)))
00049     return (x * y) / (x * y);
00050 
00051   if (hy <= 0x7fbfffff)
00052     x = __ieee754_fmod (x, 8 * y);        /* now x < 8y */
00053 
00054   if (((hx - hy) | (lx - ly)) == 0)
00055     {
00056       *quo = qs ? -1 : 1;
00057       return zero * x;
00058     }
00059 
00060   x  = fabs (x);
00061   y  = fabs (y);
00062   cquo = 0;
00063 
00064   if (x >= 4 * y)
00065     {
00066       x -= 4 * y;
00067       cquo += 4;
00068     }
00069   if (x >= 2 * y)
00070     {
00071       x -= 2 * y;
00072       cquo += 2;
00073     }
00074 
00075   if (hy < 0x00200000)
00076     {
00077       if (x + x > y)
00078        {
00079          x -= y;
00080          ++cquo;
00081          if (x + x >= y)
00082            {
00083              x -= y;
00084              ++cquo;
00085            }
00086        }
00087     }
00088   else
00089     {
00090       double y_half = 0.5 * y;
00091       if (x > y_half)
00092        {
00093          x -= y;
00094          ++cquo;
00095          if (x >= y_half)
00096            {
00097              x -= y;
00098              ++cquo;
00099            }
00100        }
00101     }
00102 
00103   *quo = qs ? -cquo : cquo;
00104 
00105   if (sx)
00106     x = -x;
00107   return x;
00108 }
00109 weak_alias (__remquo, remquo)
00110 #ifdef NO_LONG_DOUBLE
00111 strong_alias (__remquo, __remquol)
00112 weak_alias (__remquo, remquol)
00113 #endif