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. (COPYING.LIB)
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 <config.h>
33 #include <stdio.h>
34 #include <assert.h>
35 #include <stdlib.h>
36 #include <ctype.h>
37 #include <stdarg.h>
38 #include "bcmath.h"
39 #include "private.h"
40
41 /* Take the square root NUM and return it in NUM with SCALE digits
42 after the decimal place. */
43
44 int
45 bc_sqrt (bc_num *num, int scale)
46 {
47 int rscale, cmp_res, done;
48 int cscale;
49 bc_num guess, guess1, point5, diff;
50
51 /* Initial checks. */
52 cmp_res = bc_compare (*num, BCG(_zero_));
53 if (cmp_res < 0)
54 return 0; /* error */
55 else
56 {
57 if (cmp_res == 0)
58 {
59 bc_free_num (num);
60 *num = bc_copy_num (BCG(_zero_));
61 return 1;
62 }
63 }
64 cmp_res = bc_compare (*num, BCG(_one_));
65 if (cmp_res == 0)
66 {
67 bc_free_num (num);
68 *num = bc_copy_num (BCG(_one_));
69 return 1;
70 }
71
72 /* Initialize the variables. */
73 rscale = MAX (scale, (*num)->n_scale);
74 bc_init_num(&guess);
75 bc_init_num(&guess1);
76 bc_init_num(&diff);
77 point5 = bc_new_num (1,1);
78 point5->n_value[1] = 5;
79
80
81 /* Calculate the initial guess. */
82 if (cmp_res < 0)
83 {
84 /* The number is between 0 and 1. Guess should start at 1. */
85 guess = bc_copy_num (BCG(_one_));
86 cscale = (*num)->n_scale;
87 }
88 else
89 {
90 /* The number is greater than 1. Guess should start at 10^(exp/2). */
91 bc_int2num (&guess,10);
92
93 bc_int2num (&guess1,(*num)->n_len);
94 bc_multiply (guess1, point5, &guess1, 0);
95 guess1->n_scale = 0;
96 bc_raise (guess, guess1, &guess, 0);
97 bc_free_num (&guess1);
98 cscale = 3;
99 }
100
101 /* Find the square root using Newton's algorithm. */
102 done = FALSE;
103 while (!done)
104 {
105 bc_free_num (&guess1);
106 guess1 = bc_copy_num (guess);
107 bc_divide (*num, guess, &guess, cscale);
108 bc_add (guess, guess1, &guess, 0);
109 bc_multiply (guess, point5, &guess, cscale);
110 bc_sub (guess, guess1, &diff, cscale+1);
111 if (bc_is_near_zero (diff, cscale))
112 {
113 if (cscale < rscale+1)
114 cscale = MIN (cscale*3, rscale+1);
115 else
116 done = TRUE;
117 }
118 }
119
120 /* Assign the number and clean up. */
121 bc_free_num (num);
122 bc_divide (guess,BCG(_one_),num,rscale);
123 bc_free_num (&guess);
124 bc_free_num (&guess1);
125 bc_free_num (&point5);
126 bc_free_num (&diff);
127 return 1;
128 }
129