summary refs log tree commit diff
path: root/sysdeps/ieee754/sqrt.c
blob: 7e350e0d91b664201031bdc1ae8fdb13a0ffc719 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
/*
 * Copyright (c) 1985 Regents of the University of California.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms are permitted provided
 * that: (1) source distributions retain this entire copyright notice and
 * comment, and (2) distributions including binaries display the following
 * acknowledgement:  ``This product includes software developed by the
 * University of California, Berkeley and its contributors'' in the
 * documentation or other materials provided with the distribution and in
 * all advertising materials mentioning features or use of this software.
 * Neither the name of the University nor the names of its contributors may
 * be used to endorse or promote products derived from this software without
 * specific prior written permission.
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 */

#include <ansidecl.h>
#include <errno.h>
#include <math.h>

/* Return the square root of X.  */
double
DEFUN (sqrt, (x), double x)
{
  double q, s, b, r, t;
  CONST double zero = 0.0;
  int m, n, i;

  /* sqrt (NaN) is NaN; sqrt (+-0) is +-0.  */
  if (__isnan (x) || x == zero)
    return x;

  if (x < zero)
    return zero / zero;

  /* sqrt (Inf) is Inf.  */
  if (__isinf (x))
    return x;

  /* Scale X to [1,4).  */
  n = __logb (x);
  x = __scalb (x, -n);
  m = __logb (x);
  if (m != 0)
    /* Subnormal number.  */
    x = __scalb (x, -m);

  m += n;
  n = m / 2;

  if ((n + n) != m)
    {
      x *= 2;
      --m;
      n = m / 2;
    }

  /* Generate sqrt (X) bit by bit (accumulating in Q).  */
  q = 1.0;
  s = 4.0;
  x -= 1.0;
  r = 1;
  for (i = 1; i <= 51; i++)
    {
      t = s + 1;
      x *= 4;
      r /= 2;
      if (t <= x)
	{
	  s = t + t + 2, x -= t;
	  q += r;
	}
      else
	s *= 2;
    }

  /* Generate the last bit and determine the final rounding.  */
  r /= 2;
  x *= 4;
  if (x == zero)
    goto end;
  (void) (100 + r);		/* Trigger inexact flag.  */
  if (s < x)
    {
      q += r;
      x -= s;
      s += 2;
      s *= 2;
      x *= 4;
      t = (x - s) - 5;
      b = 1.0 + 3 * r / 4;
      if (b == 1.0)
	goto end;		/* B == 1: Round to zero.  */
      b = 1.0 + r / 4;
      if (b > 1.0)
	t = 1;			/* B > 1: Round to +Inf.  */
      if (t >= 0)
	q += r;
    }				/* Else round to nearest.  */
  else
    {
      s *= 2;
      x *= 4;
      t = (x - s) - 1;
      b = 1.0 + 3 * r / 4;
      if (b == 1.0)
	goto end;
      b = 1.0 + r / 4;
      if (b > 1.0)
	t = 1;
      if (t >= 0)
	q += r;
    }

end:
  return __scalb (q, n);
}