/home/ntakagi/work/tcl8.5.1/libtommath/bn_mp_sqrt.c Source File

00001 #include <tommath.h>
00002 
00003 #ifdef BN_MP_SQRT_C
00004 /* LibTomMath, multiple-precision integer library -- Tom St Denis
00005  *
00006  * LibTomMath is a library that provides multiple-precision
00007  * integer arithmetic as well as number theoretic functionality.
00008  *
00009  * The library was designed directly after the MPI library by
00010  * Michael Fromberger but has been written from scratch with
00011  * additional optimizations in place.
00012  *
00013  * The library is free for all purposes without any express
00014  * guarantee it works.
00015  *
00016  * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
00017  */
00018 
00019 #ifndef NO_FLOATING_POINT
00020 #include <math.h>
00021 #endif
00022 
00023 /* this function is less generic than mp_n_root, simpler and faster */
00024 int mp_sqrt(mp_int *arg, mp_int *ret) 
00025 {
00026   int res;
00027   mp_int t1,t2;
00028   int i, j, k;
00029 #ifndef NO_FLOATING_POINT
00030   double d;
00031   mp_digit dig;
00032 #endif
00033 
00034   /* must be positive */
00035   if (arg->sign == MP_NEG) {
00036     return MP_VAL;
00037   }
00038 
00039   /* easy out */
00040   if (mp_iszero(arg) == MP_YES) {
00041     mp_zero(ret);
00042     return MP_OKAY;
00043   }
00044   
00045   i = (arg->used / 2) - 1;
00046   j = 2 * i;
00047   if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) {
00048       return res;
00049   }
00050   
00051   if ((res = mp_init(&t2)) != MP_OKAY) {
00052     goto E2;
00053   }
00054 
00055   for (k = 0; k < i; ++k) {
00056       t1.dp[k] = (mp_digit) 0;
00057   }
00058       
00059 #ifndef NO_FLOATING_POINT
00060 
00061   /* Estimate the square root using the hardware floating point unit. */
00062 
00063   d = 0.0;
00064   for (k = arg->used-1; k >= j; --k) {
00065       d = ldexp(d, DIGIT_BIT) + (double) (arg->dp[k]);
00066   }
00067   d = sqrt(d);
00068   dig = (mp_digit) ldexp(d, -DIGIT_BIT);
00069   if (dig) {
00070       t1.used = i+2;
00071       d -= ldexp((double) dig, DIGIT_BIT);
00072       if (d != 0.0) {
00073           t1.dp[i+1] = dig;
00074           t1.dp[i] = ((mp_digit) d) - 1;
00075       } else {
00076           t1.dp[i+1] = dig-1;
00077           t1.dp[i] = MP_DIGIT_MAX;
00078       }
00079   } else {
00080       t1.used = i+1;
00081       t1.dp[i] = ((mp_digit) d) - 1;
00082   }
00083 
00084 #else
00085 
00086   /* Estimate the square root as having 1 in the most significant place. */
00087 
00088   t1.used = i + 2;
00089   t1.dp[i+1] = (mp_digit) 1;
00090   t1.dp[i] = (mp_digit) 0;
00091 
00092 #endif
00093 
00094   /* t1 > 0  */ 
00095   if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
00096     goto E1;
00097   }
00098   if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) {
00099     goto E1;
00100   }
00101   if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) {
00102     goto E1;
00103   }
00104   /* And now t1 > sqrt(arg) */
00105   do { 
00106     if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
00107       goto E1;
00108     }
00109     if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) {
00110       goto E1;
00111     }
00112     if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) {
00113       goto E1;
00114     }
00115     /* t1 >= sqrt(arg) >= t2 at this point */
00116   } while (mp_cmp_mag(&t1,&t2) == MP_GT);
00117 
00118   mp_exch(&t1,ret);
00119 
00120 E1: mp_clear(&t2);
00121 E2: mp_clear(&t1);
00122   return res;
00123 }
00124 
00125 #endif
00126 
00127 /* $Source: /cvsroot/tcl/libtommath/bn_mp_sqrt.c,v $ */
00128 /* Based on Tom's 1.3 */
00129 /* $Revision: 1.5 $ */
00130 /* $Date: 2006/12/01 05:48:23 $ */