1 /* sqrt.c: bcmath library file. */
2 /*
3 Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
4 Copyright (C) 2000 Philip A. Nelson
5
6 This library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2 of the License, or (at your option) any later version.
10
11 This 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 Lesser General Public License for more details. (LICENSE)
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with this library; if not, write to:
18
19 The Free Software Foundation, Inc.
20 59 Temple Place, Suite 330
21 Boston, MA 02111-1307 USA.
22
23 You may contact the author by:
24 e-mail: philnelson@acm.org
25 us-mail: Philip A. Nelson
26 Computer Science Department, 9062
27 Western Washington University
28 Bellingham, WA 98226-9062
29
30 *************************************************************************/
31
32 #include "bcmath.h"
33 #include <stdbool.h>
34 #include <stddef.h>
35
36 /* Take the square root NUM and return it in NUM with SCALE digits
37 after the decimal place. */
38
bc_sqrt(bc_num * num,size_t scale)39 bool bc_sqrt(bc_num *num, size_t scale)
40 {
41 /* Initial checks. */
42 int cmp_res = bc_compare(*num, BCG(_zero_));
43 if (cmp_res < 0) {
44 return false; /* error */
45 } else {
46 if (cmp_res == 0) {
47 bc_free_num (num);
48 *num = bc_copy_num(BCG(_zero_));
49 return true;
50 }
51 }
52 cmp_res = bc_compare(*num, BCG(_one_));
53 if (cmp_res == 0) {
54 bc_free_num (num);
55 *num = bc_copy_num(BCG(_one_));
56 return true;
57 }
58
59 /* Initialize the variables. */
60 size_t rscale;
61 size_t cscale;
62 bc_num guess, guess1, point5, diff;
63
64 rscale = MAX (scale, (*num)->n_scale);
65 bc_init_num(&guess1);
66 bc_init_num(&diff);
67 point5 = bc_new_num (1, 1);
68 point5->n_value[1] = 5;
69
70
71 /* Calculate the initial guess. */
72 if (cmp_res < 0) {
73 /* The number is between 0 and 1. Guess should start at 1. */
74 guess = bc_copy_num(BCG(_one_));
75 cscale = (*num)->n_scale;
76 } else {
77 /* The number is greater than 1. Guess should start at 10^(exp/2). */
78 bc_init_num(&guess);
79 bc_int2num(&guess, 10);
80
81 bc_int2num(&guess1, (*num)->n_len);
82 bc_multiply(guess1, point5, &guess1, 0);
83 guess1->n_scale = 0;
84 bc_raise_bc_exponent(guess, guess1, &guess, 0);
85 bc_free_num (&guess1);
86 cscale = 3;
87 }
88
89 /* Find the square root using Newton's algorithm. */
90 bool done = false;
91 while (!done) {
92 bc_free_num (&guess1);
93 guess1 = bc_copy_num(guess);
94 bc_divide(*num, guess, &guess, cscale);
95 bc_add(guess, guess1, &guess, 0);
96 bc_multiply(guess, point5, &guess, cscale);
97 bc_sub(guess, guess1, &diff, cscale + 1);
98 if (bc_is_near_zero(diff, cscale)) {
99 if (cscale < rscale + 1) {
100 cscale = MIN (cscale * 3, rscale + 1);
101 } else {
102 done = true;
103 }
104 }
105 }
106
107 /* Assign the number and clean up. */
108 bc_free_num (num);
109 bc_divide(guess, BCG(_one_), num, rscale);
110 bc_free_num (&guess);
111 bc_free_num (&guess1);
112 bc_free_num (&point5);
113 bc_free_num (&diff);
114 return true;
115 }
116
