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