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 const bc_num local_num = *num;
42 /* Initial checks. */
43 if (bc_is_neg(local_num)) {
44 /* Cannot take the square root of a negative number */
45 return false;
46 }
47 /* Square root of 0 is 0 */
48 if (bc_is_zero(local_num)) {
49 bc_free_num (num);
50 *num = bc_copy_num(BCG(_zero_));
51 return true;
52 }
53
54 bcmath_compare_result num_cmp_one = bc_compare(local_num, BCG(_one_), local_num->n_scale);
55 /* Square root of 1 is 1 */
56 if (num_cmp_one == BCMATH_EQUAL) {
57 bc_free_num (num);
58 *num = bc_copy_num(BCG(_one_));
59 return true;
60 }
61
62 /* Initialize the variables. */
63 size_t cscale;
64 bc_num guess, guess1, point5, diff;
65 size_t rscale = MAX(scale, local_num->n_scale);
66
67 bc_init_num(&guess1);
68 bc_init_num(&diff);
69 point5 = bc_new_num (1, 1);
70 point5->n_value[1] = 5;
71
72
73 /* Calculate the initial guess. */
74 if (num_cmp_one == BCMATH_RIGHT_GREATER) {
75 /* The number is between 0 and 1. Guess should start at 1. */
76 guess = bc_copy_num(BCG(_one_));
77 cscale = local_num->n_scale;
78 } else {
79 /* The number is greater than 1. Guess should start at 10^(exp/2). */
80 bc_init_num(&guess);
81 bc_int2num(&guess, 10);
82
83 bc_int2num(&guess1, local_num->n_len);
84 bc_multiply_ex(guess1, point5, &guess1, 0);
85 guess1->n_scale = 0;
86 bc_raise_bc_exponent(guess, guess1, &guess, 0);
87 bc_free_num (&guess1);
88 cscale = 3;
89 }
90
91 /* Find the square root using Newton's algorithm. */
92 bool done = false;
93 while (!done) {
94 bc_free_num (&guess1);
95 guess1 = bc_copy_num(guess);
96 bc_divide(*num, guess, &guess, cscale);
97 bc_add_ex(guess, guess1, &guess, 0);
98 bc_multiply_ex(guess, point5, &guess, cscale);
99 bc_sub_ex(guess, guess1, &diff, cscale + 1);
100 if (bc_is_near_zero(diff, cscale)) {
101 if (cscale < rscale + 1) {
102 cscale = MIN (cscale * 3, rscale + 1);
103 } else {
104 done = true;
105 }
106 }
107 }
108
109 /* Assign the number and clean up. */
110 bc_free_num (num);
111 bc_divide(guess, BCG(_one_), num, rscale);
112 bc_free_num (&guess);
113 bc_free_num (&guess1);
114 bc_free_num (&point5);
115 bc_free_num (&diff);
116 return true;
117 }
118
