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