Back to index

glibc  2.9
s_remquol.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 long double zero = 0.0;
00027 
00028 
00029 long double
00030 __remquol (long double x, long double p, int *quo)
00031 {
00032   int32_t ex,ep,hx,hp;
00033   u_int32_t sx,lx,lp;
00034   int cquo,qs;
00035 
00036   GET_LDOUBLE_WORDS (ex, hx, lx, x);
00037   GET_LDOUBLE_WORDS (ep, hp, lp, p);
00038   sx = ex & 0x8000;
00039   qs = (sx ^ (ep & 0x8000)) >> 15;
00040   ep &= 0x7fff;
00041   ex &= 0x7fff;
00042 
00043   /* Purge off exception values.  */
00044   if ((ep | hp | lp) == 0)
00045     return (x * p) / (x * p);                    /* p = 0 */
00046   if ((ex == 0x7fff)                      /* x not finite */
00047       || ((ep == 0x7fff)                  /* p is NaN */
00048          && ((hp | lp) != 0)))
00049     return (x * p) / (x * p);
00050 
00051   if (ep <= 0x7ffb)
00052     x = __ieee754_fmodl (x, 8 * p);              /* now x < 8p */
00053 
00054   if (((ex - ep) | (hx - hp) | (lx - lp)) == 0)
00055     {
00056       *quo = qs ? -1 : 1;
00057       return zero * x;
00058     }
00059 
00060   x  = fabsl (x);
00061   p  = fabsl (p);
00062   cquo = 0;
00063 
00064   if (x >= 4 * p)
00065     {
00066       x -= 4 * p;
00067       cquo += 4;
00068     }
00069   if (x >= 2 * p)
00070     {
00071       x -= 2 * p;
00072       cquo += 2;
00073     }
00074 
00075   if (ep < 0x0002)
00076     {
00077       if (x + x > p)
00078        {
00079          x -= p;
00080          ++cquo;
00081          if (x + x >= p)
00082            {
00083              x -= p;
00084              ++cquo;
00085            }
00086        }
00087     }
00088   else
00089     {
00090       long double p_half = 0.5 * p;
00091       if (x > p_half)
00092        {
00093          x -= p;
00094          ++cquo;
00095          if (x >= p_half)
00096            {
00097              x -= p;
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 (__remquol, remquol)