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 small number theoretic routines
37 * mrsmall.c
38 */
39
40 #include "miracl.h"
41
42 #ifdef MR_FP
43 #include <math.h>
44 #endif
45
46 #ifdef MR_WIN64
47 #include <intrin.h>
48 #endif
49
50
smul(mr_small x,mr_small y,mr_small n)51 mr_small smul(mr_small x,mr_small y,mr_small n)
52 { /* returns x*y mod n */
53 mr_small r;
54
55 #ifdef MR_ITANIUM
56 mr_small tm;
57 #endif
58 #ifdef MR_WIN64
59 mr_small tm;
60 #endif
61 #ifdef MR_FP
62 mr_small dres;
63 #endif
64 #ifndef MR_NOFULLWIDTH
65 if (n==0)
66 { /* Assume n=2^MIRACL */
67 muldvd(x,y,(mr_small)0,&r);
68 return r;
69 }
70 #endif
71 x=MR_REMAIN(x,n);
72 y=MR_REMAIN(y,n);
73 muldiv(x,y,(mr_small)0,n,&r);
74 return r;
75 }
76
invers(mr_small x,mr_small y)77 mr_small invers(mr_small x,mr_small y)
78 { /* returns inverse of x mod y */
79 mr_small r,s,q,t,p;
80 #ifdef MR_FP
81 mr_small dres;
82 #endif
83
84 BOOL pos;
85 if (y!=0) x=MR_REMAIN(x,y);
86 r=1;
87 s=0;
88 p=y;
89 pos=TRUE;
90 #ifndef MR_NOFULLWIDTH
91 if (p==0)
92 { /* if modulus is 0, assume its actually 2^MIRACL */
93 if (x==1) return (mr_small)1;
94 t=r; r=s; s=t;
95 p=x;
96 q=muldvm((mr_small)1,(mr_small)0,p,&t);
97 t=r+s*q; r=s; s=t;
98 t=0-p*q; x=p; p=t;
99 }
100 #endif
101 while (p!=0)
102 { /* main euclidean loop */
103 q=MR_DIV(x,p);
104 t=r+s*q; r=s; s=t;
105 t=x-p*q; x=p; p=t;
106 pos=!pos;
107 }
108 if (!pos) r=y-r;
109 return r;
110 }
111
jac(mr_small x,mr_small n)112 int jac(mr_small x,mr_small n)
113 { /* finds (x/n) as (-1)^m */
114 int m,k,n8,u4;
115 mr_small t;
116 #ifdef MR_FP
117 mr_small dres;
118 #endif
119 if (x==0)
120 {
121 if (n==1) return 1;
122 else return 0;
123 }
124 if (MR_REMAIN(n,2)==0) return 0;
125 x=MR_REMAIN(x,n);
126 m=0;
127 while(n>1)
128 { /* main loop */
129 if (x==0) return 0;
130
131 /* extract powers of 2 */
132 for (k=0;MR_REMAIN(x,2)==0;k++) x=MR_DIV(x,2);
133 n8=(int)MR_REMAIN(n,8);
134 if (k%2==1) m+=(n8*n8-1)/8;
135
136 /* quadratic reciprocity */
137 u4=(int)MR_REMAIN(x,4);
138 m+=(n8-1)*(u4-1)/4;
139 t=n; t=MR_REMAIN(t,x);
140 n=x; x=t;
141 m%=2;
142 }
143 if (m==0) return 1;
144 else return (-1);
145 }
146
147 #ifndef MR_STATIC
148
spmd(mr_small x,mr_small n,mr_small m)149 mr_small spmd(mr_small x,mr_small n,mr_small m)
150 { /* returns x^n mod m */
151 mr_small r,sx;
152 #ifdef MR_FP
153 mr_small dres;
154 #endif
155 x=MR_REMAIN(x,m);
156 r=0;
157 if (x==0) return r;
158 r=1;
159 if (n==0) return r;
160 sx=x;
161 forever
162 {
163 if (MR_REMAIN(n,2)!=0) muldiv(r,sx,(mr_small)0,m,&r);
164 n=MR_DIV(n,2);
165 if (n==0) return r;
166 muldiv(sx,sx,(mr_small)0,m,&sx);
167 }
168 }
169
sqrmp(mr_small x,mr_small m)170 mr_small sqrmp(mr_small x,mr_small m)
171 { /* square root mod a small prime by Shanks method *
172 * returns 0 if root does not exist or m not prime */
173 mr_small z,y,v,w,t,q;
174 #ifdef MR_FP
175 mr_small dres;
176 #endif
177 int i,e,n,r;
178 BOOL pp;
179 x=MR_REMAIN(x,m);
180 if (x==0) return 0;
181 if (x==1) return 1;
182 if (spmd(x,(mr_small)((m-1)/2),m)!=1) return 0; /* Legendre symbol not 1 */
183 if (MR_REMAIN(m,4)==3) return spmd(x,(mr_small)((m+1)/4),m); /* easy case for m=4.k+3 */
184 if (MR_REMAIN(m,8)==5)
185 { /* also relatively easy */
186 t=spmd(x,(mr_small)((m-1)/4),m);
187 if (t==1) return spmd(x,(mr_small)((m+3)/8),m);
188 if (t==(mr_small)(m-1))
189 {
190 muldiv((mr_small)4,x,(mr_small)0,m,&t);
191 t=spmd(t,(mr_small)((m+3)/8),m);
192 muldiv(t,(mr_small)((m+1)/2),(mr_small)0,m,&t);
193 return t;
194 }
195 return 0;
196 }
197 q=m-1;
198 e=0;
199 while (MR_REMAIN(q,2)==0)
200 {
201 q=MR_DIV(q,2);
202 e++;
203 }
204 if (e==0) return 0; /* even m */
205 for (r=2;;r++)
206 { /* find suitable z */
207 z=spmd((mr_small)r,q,m);
208 if (z==1) continue;
209 t=z;
210 pp=FALSE;
211 for (i=1;i<e;i++)
212 { /* check for composite m */
213 if (t==(m-1)) pp=TRUE;
214 muldiv(t,t,(mr_small)0,m,&t);
215 if (t==1 && !pp) return 0;
216 }
217 if (t==(m-1)) break;
218 if (!pp) return 0; /* m is not prime */
219 }
220 y=z;
221 r=e;
222 v=spmd(x,(mr_small)((q+1)/2),m);
223 w=spmd(x,q,m);
224 while (w!=1)
225 {
226 t=w;
227 for (n=0;t!=1;n++) muldiv(t,t,(mr_small)0,m,&t);
228 if (n>=r) return 0;
229 y=spmd(y,mr_shiftbits(1,r-n-1),m);
230 muldiv(v,y,(mr_small)0,m,&v);
231 muldiv(y,y,(mr_small)0,m,&y);
232 muldiv(w,y,(mr_small)0,m,&w);
233 r=n;
234 }
235 return v;
236 }
237
238 #endif
239