Return to sqrt.c CVS log | Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / php / ext / bcmath / libbcmath / src |
1.1 misho 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: { 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 TSRMLS_CC); 75: bc_init_num(&guess1 TSRMLS_CC); 76: bc_init_num(&diff TSRMLS_CC); 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 TSRMLS_CC); 95: guess1->n_scale = 0; 96: bc_raise (guess, guess1, &guess, 0 TSRMLS_CC); 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 TSRMLS_CC); 108: bc_add (guess, guess1, &guess, 0); 109: bc_multiply (guess, point5, &guess, cscale TSRMLS_CC); 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 TSRMLS_CC); 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: