Back to index

php5  5.3.10
sqrt.c
Go to the documentation of this file.
00001 /* sqrt.c: bcmath library file. */
00002 /*
00003     Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
00004     Copyright (C) 2000 Philip A. Nelson
00005 
00006     This library is free software; you can redistribute it and/or
00007     modify it under the terms of the GNU Lesser General Public
00008     License as published by the Free Software Foundation; either
00009     version 2 of the License, or (at your option) any later version.
00010 
00011     This library is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014     Lesser General Public License for more details.  (COPYING.LIB)
00015 
00016     You should have received a copy of the GNU Lesser General Public
00017     License along with this library; if not, write to:
00018 
00019       The Free Software Foundation, Inc.
00020       59 Temple Place, Suite 330
00021       Boston, MA 02111-1307 USA.
00022 
00023     You may contact the author by:
00024        e-mail:  philnelson@acm.org
00025       us-mail:  Philip A. Nelson
00026                 Computer Science Department, 9062
00027                 Western Washington University
00028                 Bellingham, WA 98226-9062
00029        
00030 *************************************************************************/
00031 
00032 #include <config.h>
00033 #include <stdio.h>
00034 #include <assert.h>
00035 #include <stdlib.h>
00036 #include <ctype.h>
00037 #include <stdarg.h>
00038 #include "bcmath.h"
00039 #include "private.h"
00040 
00041 /* Take the square root NUM and return it in NUM with SCALE digits
00042    after the decimal place. */
00043 
00044 int
00045 bc_sqrt (bc_num *num, int scale TSRMLS_DC)
00046 {
00047   int rscale, cmp_res, done;
00048   int cscale;
00049   bc_num guess, guess1, point5, diff;
00050 
00051   /* Initial checks. */
00052   cmp_res = bc_compare (*num, BCG(_zero_));
00053   if (cmp_res < 0)
00054     return 0;        /* error */
00055   else
00056     {
00057       if (cmp_res == 0)
00058        {
00059          bc_free_num (num);
00060          *num = bc_copy_num (BCG(_zero_));
00061          return 1;
00062        }
00063     }
00064   cmp_res = bc_compare (*num, BCG(_one_));
00065   if (cmp_res == 0)
00066     {
00067       bc_free_num (num);
00068       *num = bc_copy_num (BCG(_one_));
00069       return 1;
00070     }
00071 
00072   /* Initialize the variables. */
00073   rscale = MAX (scale, (*num)->n_scale);
00074   bc_init_num(&guess TSRMLS_CC);
00075   bc_init_num(&guess1 TSRMLS_CC);
00076   bc_init_num(&diff TSRMLS_CC);
00077   point5 = bc_new_num (1,1);
00078   point5->n_value[1] = 5;
00079 
00080 
00081   /* Calculate the initial guess. */
00082   if (cmp_res < 0)
00083     {
00084       /* The number is between 0 and 1.  Guess should start at 1. */
00085       guess = bc_copy_num (BCG(_one_));
00086       cscale = (*num)->n_scale;
00087     }
00088   else
00089     {
00090       /* The number is greater than 1.  Guess should start at 10^(exp/2). */
00091       bc_int2num (&guess,10);
00092 
00093       bc_int2num (&guess1,(*num)->n_len);
00094       bc_multiply (guess1, point5, &guess1, 0 TSRMLS_CC);
00095       guess1->n_scale = 0;
00096       bc_raise (guess, guess1, &guess, 0 TSRMLS_CC);
00097       bc_free_num (&guess1);
00098       cscale = 3;
00099     }
00100 
00101   /* Find the square root using Newton's algorithm. */
00102   done = FALSE;
00103   while (!done)
00104     {
00105       bc_free_num (&guess1);
00106       guess1 = bc_copy_num (guess);
00107       bc_divide (*num, guess, &guess, cscale TSRMLS_CC);
00108       bc_add (guess, guess1, &guess, 0);
00109       bc_multiply (guess, point5, &guess, cscale TSRMLS_CC);
00110       bc_sub (guess, guess1, &diff, cscale+1);
00111       if (bc_is_near_zero (diff, cscale))
00112        {
00113          if (cscale < rscale+1)
00114            cscale = MIN (cscale*3, rscale+1);
00115          else
00116            done = TRUE;
00117        }
00118     }
00119 
00120   /* Assign the number and clean up. */
00121   bc_free_num (num);
00122   bc_divide (guess,BCG(_one_),num,rscale TSRMLS_CC);
00123   bc_free_num (&guess);
00124   bc_free_num (&guess1);
00125   bc_free_num (&point5);
00126   bc_free_num (&diff);
00127   return 1;
00128 }
00129