xref: /PHP-8.1/ext/bcmath/libbcmath/src/sqrt.c (revision 3d2819e4)
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(&guess1);
74   bc_init_num(&diff);
75   point5 = bc_new_num (1,1);
76   point5->n_value[1] = 5;
77 
78 
79   /* Calculate the initial guess. */
80   if (cmp_res < 0)
81     {
82       /* The number is between 0 and 1.  Guess should start at 1. */
83       guess = bc_copy_num (BCG(_one_));
84       cscale = (*num)->n_scale;
85     }
86   else
87     {
88       /* The number is greater than 1.  Guess should start at 10^(exp/2). */
89       bc_init_num(&guess);
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