1
2 /***************************************************************************
3 *
4 Copyright 2012 CertiVox IOM Ltd. *
5 *
6 This file is part of CertiVox MIRACL Crypto SDK. *
7 *
8 The CertiVox MIRACL Crypto SDK provides developers with an *
9 extensive and efficient set of cryptographic functions. *
10 For further information about its features and functionalities please *
11 refer to http://www.certivox.com *
12 *
13 * The CertiVox MIRACL Crypto SDK is free software: you can *
14 redistribute it and/or modify it under the terms of the *
15 GNU Affero General Public License as published by the *
16 Free Software Foundation, either version 3 of the License, *
17 or (at your option) any later version. *
18 *
19 * The CertiVox MIRACL Crypto SDK is distributed in the hope *
20 that it will be useful, but WITHOUT ANY WARRANTY; without even the *
21 implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
22 See the GNU Affero General Public License for more details. *
23 *
24 * You should have received a copy of the GNU Affero General Public *
25 License along with CertiVox MIRACL Crypto SDK. *
26 If not, see <http://www.gnu.org/licenses/>. *
27 *
28 You can be released from the requirements of the license by purchasing *
29 a commercial license. Buying such a license is mandatory as soon as you *
30 develop commercial activities involving the CertiVox MIRACL Crypto SDK *
31 without disclosing the source code of your own applications, or shipping *
32 the CertiVox MIRACL Crypto SDK with a closed source product. *
33 *
34 ***************************************************************************/
35 /*
36 * MIRACL arithmetic routines 3 - simple powers and roots
37 * mrarth3.c
38 */
39
40 #include <stdlib.h>
41 #include "miracl.h"
42
expint(_MIPD_ int b,int n,big x)43 void expint(_MIPD_ int b,int n,big x)
44 { /* sets x=b^n */
45 unsigned int bit,un;
46 #ifdef MR_OS_THREADS
47 miracl *mr_mip=get_mip();
48 #endif
49 if (mr_mip->ERNUM) return;
50 convert(_MIPP_ 1,x);
51 if (n==0) return;
52
53 MR_IN(50)
54
55 if (n<0)
56 {
57 mr_berror(_MIPP_ MR_ERR_NEG_POWER);
58 MR_OUT
59 return;
60 }
61 if (b==2) expb2(_MIPP_ n,x);
62 else
63 {
64 bit=1;
65 un=(unsigned int)n;
66 while (un>=bit) bit<<=1;
67 bit>>=1;
68 while (bit>0)
69 { /* ltr method */
70 multiply(_MIPP_ x,x,x);
71 if ((bit&un)!=0) premult(_MIPP_ x,b,x);
72 bit>>=1;
73 }
74 }
75 MR_OUT
76 }
77
power(_MIPD_ big x,long n,big z,big w)78 void power(_MIPD_ big x,long n,big z,big w)
79 { /* raise big number to int power w=x^n *
80 * (mod z if z and w distinct) */
81 mr_small norm;
82 #ifdef MR_OS_THREADS
83 miracl *mr_mip=get_mip();
84 #endif
85
86 copy(x,mr_mip->w5);
87 zero(w);
88 if(mr_mip->ERNUM || size(mr_mip->w5)==0) return;
89 convert(_MIPP_ 1,w);
90 if (n==0L) return;
91
92 MR_IN(17)
93
94 if (n<0L)
95 {
96 mr_berror(_MIPP_ MR_ERR_NEG_POWER);
97 MR_OUT
98 return;
99 }
100
101 if (w==z) forever
102 { /* "Russian peasant" exponentiation */
103 if (n%2!=0L)
104 multiply(_MIPP_ w,mr_mip->w5,w);
105 n/=2L;
106 if (mr_mip->ERNUM || n==0L) break;
107 multiply(_MIPP_ mr_mip->w5,mr_mip->w5,mr_mip->w5);
108 }
109 else
110 {
111 norm=normalise(_MIPP_ z,z);
112 divide(_MIPP_ mr_mip->w5,z,z);
113 forever
114 {
115 if (mr_mip->user!=NULL) (*mr_mip->user)();
116
117 if (n%2!=0L) mad(_MIPP_ w,mr_mip->w5,mr_mip->w5,z,z,w);
118 n/=2L;
119 if (mr_mip->ERNUM || n==0L) break;
120 mad(_MIPP_ mr_mip->w5,mr_mip->w5,mr_mip->w5,z,z,mr_mip->w5);
121 }
122 if (norm!=1)
123 {
124 #ifdef MR_FP_ROUNDING
125 mr_sdiv(_MIPP_ z,norm,mr_invert(norm),z);
126 #else
127 mr_sdiv(_MIPP_ z,norm,z);
128 #endif
129 divide(_MIPP_ w,z,z);
130 }
131 }
132
133 MR_OUT
134 }
135
nroot(_MIPD_ big x,int n,big w)136 BOOL nroot(_MIPD_ big x,int n,big w)
137 { /* extract lower approximation to nth root *
138 * w=x^(1/n) returns TRUE for exact root *
139 * uses Newtons method */
140 int sx,dif,s,p,d,lg2,lgx,rem;
141 BOOL full;
142 #ifdef MR_OS_THREADS
143 miracl *mr_mip=get_mip();
144 #endif
145 if (mr_mip->ERNUM) return FALSE;
146 if (size(x)==0 || n==1)
147 {
148 copy(x,w);
149 return TRUE;
150 }
151
152 MR_IN(16)
153
154 if (n<1) mr_berror(_MIPP_ MR_ERR_BAD_ROOT);
155 sx=exsign(x);
156 if (n%2==0 && sx==MINUS) mr_berror(_MIPP_ MR_ERR_NEG_ROOT);
157 if (mr_mip->ERNUM)
158 {
159 MR_OUT
160 return FALSE;
161 }
162 insign(PLUS,x);
163 lgx=logb2(_MIPP_ x);
164 if (n>=lgx)
165 { /* root must be 1 */
166 insign(sx,x);
167 convert(_MIPP_ sx,w);
168 MR_OUT
169 if (lgx==1) return TRUE;
170 else return FALSE;
171 }
172 expb2(_MIPP_ 1+(lgx-1)/n,mr_mip->w2); /* guess root as 2^(log2(x)/n) */
173 s=(-(((int)x->len-1)/n)*n);
174 mr_shift(_MIPP_ mr_mip->w2,s/n,mr_mip->w2);
175 lg2=logb2(_MIPP_ mr_mip->w2)-1;
176 full=FALSE;
177 if (s==0) full=TRUE;
178 d=0;
179 p=1;
180 while (!mr_mip->ERNUM)
181 { /* Newtons method */
182 copy(mr_mip->w2,mr_mip->w3);
183 mr_shift(_MIPP_ x,s,mr_mip->w4);
184 mr_mip->check=OFF;
185 power(_MIPP_ mr_mip->w2,n-1,mr_mip->w6,mr_mip->w6);
186 mr_mip->check=ON;
187 divide(_MIPP_ mr_mip->w4,mr_mip->w6,mr_mip->w2);
188 rem=size(mr_mip->w4);
189 subtract(_MIPP_ mr_mip->w2,mr_mip->w3,mr_mip->w2);
190 dif=size(mr_mip->w2);
191 subdiv(_MIPP_ mr_mip->w2,n,mr_mip->w2);
192 add(_MIPP_ mr_mip->w2,mr_mip->w3,mr_mip->w2);
193 p*=2;
194 if(p<lg2+d*mr_mip->lg2b) continue;
195 if (full && mr_abs(dif)<n)
196 { /* test for finished */
197 while (dif<0)
198 {
199 rem=0;
200 decr(_MIPP_ mr_mip->w2,1,mr_mip->w2);
201 mr_mip->check=OFF;
202 power(_MIPP_ mr_mip->w2,n,mr_mip->w6,mr_mip->w6);
203 mr_mip->check=ON;
204 dif=mr_compare(x,mr_mip->w6);
205 }
206 copy(mr_mip->w2,w);
207 insign(sx,w);
208 insign(sx,x);
209 MR_OUT
210 if (rem==0 && dif==0) return TRUE;
211 else return FALSE;
212 }
213 else
214 { /* adjust precision */
215 d*=2;
216 if (d==0) d=1;
217 s+=d*n;
218 if (s>=0)
219 {
220 d-=s/n;
221 s=0;
222 full=TRUE;
223 }
224 mr_shift(_MIPP_ mr_mip->w2,d,mr_mip->w2);
225 }
226 p/=2;
227 }
228 MR_OUT
229 return FALSE;
230 }
231
232