(__ieee754_sqrt): Remove unused variables b and n.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-96 / s_remquol.c
1 /* Compute remainder and a congruent to the quotient.
2    Copyright (C) 1997 Free Software Foundation, Inc.
3    This file is part of the GNU C Library.
4    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
5
6    The GNU C Library is free software; you can redistribute it and/or
7    modify it under the terms of the GNU Library General Public License as
8    published by the Free Software Foundation; either version 2 of the
9    License, or (at your option) any later version.
10
11    The GNU C Library is distributed in the hope that it will be useful,
12    but WITHOUT ANY WARRANTY; without even the implied warranty of
13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14    Library General Public License for more details.
15
16    You should have received a copy of the GNU Library General Public
17    License along with the GNU C Library; see the file COPYING.LIB.  If not,
18    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
19    Boston, MA 02111-1307, USA.  */
20
21 #include <math.h>
22
23 #include "math_private.h"
24
25
26 static const long double zero = 0.0;
27
28
29 long double
30 __remquol (long double x, long double p, int *quo)
31 {
32   int32_t ex,ep,hx,hp;
33   u_int32_t sx,lx,lp;
34   int cquo,qs;
35
36   GET_LDOUBLE_WORDS (ex, hx, lx, x);
37   GET_LDOUBLE_WORDS (ep, hp, lp, p);
38   sx = ex & 0x8000;
39   qs = (sx ^ (ep & 0x8000)) >> 15;
40   ep &= 0x7fff;
41   ex &= 0x7fff;
42
43   /* Purge off exception values.  */
44   if ((ep | hp | lp) == 0)
45     return (x * p) / (x * p);                   /* p = 0 */
46   if ((ex == 0x7fff)                            /* x not finite */
47       || ((ep == 0x7fff)                        /* p is NaN */
48           && ((hp | lp) != 0)))
49     return (x * p) / (x * p);
50
51   if (ep <= 0x7ffb)
52     x = __ieee754_fmodl (x, 8 * p);             /* now x < 8p */
53
54   if (((ex - ep) | (hx - hp) | (lx - lp)) == 0)
55     {
56       *quo = qs ? -1 : 1;
57       return zero * x;
58     }
59
60   x  = fabsl (x);
61   p  = fabsl (p);
62   cquo = 0;
63
64   if (x >= 4 * p)
65     {
66       x -= 4 * p;
67       cquo += 4;
68     }
69   if (x >= 2 * p)
70     {
71       x -= 2 * p;
72       cquo += 2;
73     }
74
75   if (ep < 0x0002)
76     {
77       if (x + x > p)
78         {
79           x -= p;
80           ++cquo;
81           if (x + x >= p)
82             {
83               x -= p;
84               ++cquo;
85             }
86         }
87     }
88   else
89     {
90       long double p_half = 0.5 * p;
91       if (x > p_half)
92         {
93           x -= p;
94           ++cquo;
95           if (x >= p_half)
96             {
97               x -= p;
98               ++cquo;
99             }
100         }
101     }
102
103   *quo = qs ? -cquo : cquo;
104
105   if (sx)
106     x = -x;
107   return x;
108 }
109 weak_alias (__remquol, remquol)