| |
1 /* This Source Code Form is subject to the terms of the Mozilla Public |
| |
2 * License, v. 2.0. If a copy of the MPL was not distributed with this |
| |
3 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ |
| |
4 |
| |
5 #ifdef SOLARIS |
| |
6 #define RF_INLINE_MACROS 1 |
| |
7 #endif |
| |
8 |
| |
9 static const double TwoTo16=65536.0; |
| |
10 static const double TwoToMinus16=1.0/65536.0; |
| |
11 static const double Zero=0.0; |
| |
12 static const double TwoTo32=65536.0*65536.0; |
| |
13 static const double TwoToMinus32=1.0/(65536.0*65536.0); |
| |
14 |
| |
15 #ifdef RF_INLINE_MACROS |
| |
16 |
| |
17 double upper32(double); |
| |
18 double lower32(double, double); |
| |
19 double mod(double, double, double); |
| |
20 |
| |
21 void i16_to_d16_and_d32x4(const double * /*1/(2^16)*/, |
| |
22 const double * /* 2^16*/, |
| |
23 const double * /* 0 */, |
| |
24 double * /*result16*/, |
| |
25 double * /* result32 */, |
| |
26 float * /*source - should be unsigned int* |
| |
27 converted to float* */); |
| |
28 |
| |
29 #else |
| |
30 #ifdef MP_USE_FLOOR |
| |
31 #include <math.h> |
| |
32 #else |
| |
33 #define floor(d) ((double)((unsigned long long)(d))) |
| |
34 #endif |
| |
35 |
| |
36 static double upper32(double x) |
| |
37 { |
| |
38 return floor(x*TwoToMinus32); |
| |
39 } |
| |
40 |
| |
41 static double lower32(double x, double y) |
| |
42 { |
| |
43 return x-TwoTo32*floor(x*TwoToMinus32); |
| |
44 } |
| |
45 |
| |
46 static double mod(double x, double oneoverm, double m) |
| |
47 { |
| |
48 return x-m*floor(x*oneoverm); |
| |
49 } |
| |
50 |
| |
51 #endif |
| |
52 |
| |
53 |
| |
54 static void cleanup(double *dt, int from, int tlen) |
| |
55 { |
| |
56 int i; |
| |
57 double tmp,tmp1,x,x1; |
| |
58 |
| |
59 tmp=tmp1=Zero; |
| |
60 /* original code ** |
| |
61 for(i=2*from;i<2*tlen-2;i++) |
| |
62 { |
| |
63 x=dt[i]; |
| |
64 dt[i]=lower32(x,Zero)+tmp1; |
| |
65 tmp1=tmp; |
| |
66 tmp=upper32(x); |
| |
67 } |
| |
68 dt[tlen-2]+=tmp1; |
| |
69 dt[tlen-1]+=tmp; |
| |
70 **end original code ***/ |
| |
71 /* new code ***/ |
| |
72 for(i=2*from;i<2*tlen;i+=2) |
| |
73 { |
| |
74 x=dt[i]; |
| |
75 x1=dt[i+1]; |
| |
76 dt[i]=lower32(x,Zero)+tmp; |
| |
77 dt[i+1]=lower32(x1,Zero)+tmp1; |
| |
78 tmp=upper32(x); |
| |
79 tmp1=upper32(x1); |
| |
80 } |
| |
81 /** end new code **/ |
| |
82 } |
| |
83 |
| |
84 |
| |
85 void conv_d16_to_i32(unsigned int *i32, double *d16, long long *tmp, int ilen) |
| |
86 { |
| |
87 int i; |
| |
88 long long t, t1, a, b, c, d; |
| |
89 |
| |
90 t1=0; |
| |
91 a=(long long)d16[0]; |
| |
92 b=(long long)d16[1]; |
| |
93 for(i=0; i<ilen-1; i++) |
| |
94 { |
| |
95 c=(long long)d16[2*i+2]; |
| |
96 t1+=(unsigned int)a; |
| |
97 t=(a>>32); |
| |
98 d=(long long)d16[2*i+3]; |
| |
99 t1+=(b&0xffff)<<16; |
| |
100 t+=(b>>16)+(t1>>32); |
| |
101 i32[i]=(unsigned int)t1; |
| |
102 t1=t; |
| |
103 a=c; |
| |
104 b=d; |
| |
105 } |
| |
106 t1+=(unsigned int)a; |
| |
107 t=(a>>32); |
| |
108 t1+=(b&0xffff)<<16; |
| |
109 i32[i]=(unsigned int)t1; |
| |
110 } |
| |
111 |
| |
112 void conv_i32_to_d32(double *d32, unsigned int *i32, int len) |
| |
113 { |
| |
114 int i; |
| |
115 |
| |
116 #pragma pipeloop(0) |
| |
117 for(i=0;i<len;i++) d32[i]=(double)(i32[i]); |
| |
118 } |
| |
119 |
| |
120 |
| |
121 void conv_i32_to_d16(double *d16, unsigned int *i32, int len) |
| |
122 { |
| |
123 int i; |
| |
124 unsigned int a; |
| |
125 |
| |
126 #pragma pipeloop(0) |
| |
127 for(i=0;i<len;i++) |
| |
128 { |
| |
129 a=i32[i]; |
| |
130 d16[2*i]=(double)(a&0xffff); |
| |
131 d16[2*i+1]=(double)(a>>16); |
| |
132 } |
| |
133 } |
| |
134 |
| |
135 |
| |
136 void conv_i32_to_d32_and_d16(double *d32, double *d16, |
| |
137 unsigned int *i32, int len) |
| |
138 { |
| |
139 int i = 0; |
| |
140 unsigned int a; |
| |
141 |
| |
142 #pragma pipeloop(0) |
| |
143 #ifdef RF_INLINE_MACROS |
| |
144 for(;i<len-3;i+=4) |
| |
145 { |
| |
146 i16_to_d16_and_d32x4(&TwoToMinus16, &TwoTo16, &Zero, |
| |
147 &(d16[2*i]), &(d32[i]), (float *)(&(i32[i]))); |
| |
148 } |
| |
149 #endif |
| |
150 for(;i<len;i++) |
| |
151 { |
| |
152 a=i32[i]; |
| |
153 d32[i]=(double)(i32[i]); |
| |
154 d16[2*i]=(double)(a&0xffff); |
| |
155 d16[2*i+1]=(double)(a>>16); |
| |
156 } |
| |
157 } |
| |
158 |
| |
159 |
| |
160 void adjust_montf_result(unsigned int *i32, unsigned int *nint, int len) |
| |
161 { |
| |
162 long long acc; |
| |
163 int i; |
| |
164 |
| |
165 if(i32[len]>0) i=-1; |
| |
166 else |
| |
167 { |
| |
168 for(i=len-1; i>=0; i--) |
| |
169 { |
| |
170 if(i32[i]!=nint[i]) break; |
| |
171 } |
| |
172 } |
| |
173 if((i<0)||(i32[i]>nint[i])) |
| |
174 { |
| |
175 acc=0; |
| |
176 for(i=0;i<len;i++) |
| |
177 { |
| |
178 acc=acc+(unsigned long long)(i32[i])-(unsigned long long)(nint[i]); |
| |
179 i32[i]=(unsigned int)acc; |
| |
180 acc=acc>>32; |
| |
181 } |
| |
182 } |
| |
183 } |
| |
184 |
| |
185 |
| |
186 |
| |
187 |
| |
188 /* |
| |
189 ** the lengths of the input arrays should be at least the following: |
| |
190 ** result[nlen+1], dm1[nlen], dm2[2*nlen+1], dt[4*nlen+2], dn[nlen], nint[nlen] |
| |
191 ** all of them should be different from one another |
| |
192 ** |
| |
193 */ |
| |
194 void mont_mulf_noconv(unsigned int *result, |
| |
195 double *dm1, double *dm2, double *dt, |
| |
196 double *dn, unsigned int *nint, |
| |
197 int nlen, double dn0) |
| |
198 { |
| |
199 int i, j, jj; |
| |
200 int tmp; |
| |
201 double digit, m2j, nextm2j, a, b; |
| |
202 double *dptmp, *pdm1, *pdm2, *pdn, *pdtj, pdn_0, pdm1_0; |
| |
203 |
| |
204 pdm1=&(dm1[0]); |
| |
205 pdm2=&(dm2[0]); |
| |
206 pdn=&(dn[0]); |
| |
207 pdm2[2*nlen]=Zero; |
| |
208 |
| |
209 if (nlen!=16) |
| |
210 { |
| |
211 for(i=0;i<4*nlen+2;i++) dt[i]=Zero; |
| |
212 |
| |
213 a=dt[0]=pdm1[0]*pdm2[0]; |
| |
214 digit=mod(lower32(a,Zero)*dn0,TwoToMinus16,TwoTo16); |
| |
215 |
| |
216 pdtj=&(dt[0]); |
| |
217 for(j=jj=0;j<2*nlen;j++,jj++,pdtj++) |
| |
218 { |
| |
219 m2j=pdm2[j]; |
| |
220 a=pdtj[0]+pdn[0]*digit; |
| |
221 b=pdtj[1]+pdm1[0]*pdm2[j+1]+a*TwoToMinus16; |
| |
222 pdtj[1]=b; |
| |
223 |
| |
224 #pragma pipeloop(0) |
| |
225 for(i=1;i<nlen;i++) |
| |
226 { |
| |
227 pdtj[2*i]+=pdm1[i]*m2j+pdn[i]*digit; |
| |
228 } |
| |
229 if((jj==30)) {cleanup(dt,j/2+1,2*nlen+1); jj=0;} |
| |
230 |
| |
231 digit=mod(lower32(b,Zero)*dn0,TwoToMinus16,TwoTo16); |
| |
232 } |
| |
233 } |
| |
234 else |
| |
235 { |
| |
236 a=dt[0]=pdm1[0]*pdm2[0]; |
| |
237 |
| |
238 dt[65]= dt[64]= dt[63]= dt[62]= dt[61]= dt[60]= |
| |
239 dt[59]= dt[58]= dt[57]= dt[56]= dt[55]= dt[54]= |
| |
240 dt[53]= dt[52]= dt[51]= dt[50]= dt[49]= dt[48]= |
| |
241 dt[47]= dt[46]= dt[45]= dt[44]= dt[43]= dt[42]= |
| |
242 dt[41]= dt[40]= dt[39]= dt[38]= dt[37]= dt[36]= |
| |
243 dt[35]= dt[34]= dt[33]= dt[32]= dt[31]= dt[30]= |
| |
244 dt[29]= dt[28]= dt[27]= dt[26]= dt[25]= dt[24]= |
| |
245 dt[23]= dt[22]= dt[21]= dt[20]= dt[19]= dt[18]= |
| |
246 dt[17]= dt[16]= dt[15]= dt[14]= dt[13]= dt[12]= |
| |
247 dt[11]= dt[10]= dt[ 9]= dt[ 8]= dt[ 7]= dt[ 6]= |
| |
248 dt[ 5]= dt[ 4]= dt[ 3]= dt[ 2]= dt[ 1]=Zero; |
| |
249 |
| |
250 pdn_0=pdn[0]; |
| |
251 pdm1_0=pdm1[0]; |
| |
252 |
| |
253 digit=mod(lower32(a,Zero)*dn0,TwoToMinus16,TwoTo16); |
| |
254 pdtj=&(dt[0]); |
| |
255 |
| |
256 for(j=0;j<32;j++,pdtj++) |
| |
257 { |
| |
258 |
| |
259 m2j=pdm2[j]; |
| |
260 a=pdtj[0]+pdn_0*digit; |
| |
261 b=pdtj[1]+pdm1_0*pdm2[j+1]+a*TwoToMinus16; |
| |
262 pdtj[1]=b; |
| |
263 |
| |
264 /**** this loop will be fully unrolled: |
| |
265 for(i=1;i<16;i++) |
| |
266 { |
| |
267 pdtj[2*i]+=pdm1[i]*m2j+pdn[i]*digit; |
| |
268 } |
| |
269 *************************************/ |
| |
270 pdtj[2]+=pdm1[1]*m2j+pdn[1]*digit; |
| |
271 pdtj[4]+=pdm1[2]*m2j+pdn[2]*digit; |
| |
272 pdtj[6]+=pdm1[3]*m2j+pdn[3]*digit; |
| |
273 pdtj[8]+=pdm1[4]*m2j+pdn[4]*digit; |
| |
274 pdtj[10]+=pdm1[5]*m2j+pdn[5]*digit; |
| |
275 pdtj[12]+=pdm1[6]*m2j+pdn[6]*digit; |
| |
276 pdtj[14]+=pdm1[7]*m2j+pdn[7]*digit; |
| |
277 pdtj[16]+=pdm1[8]*m2j+pdn[8]*digit; |
| |
278 pdtj[18]+=pdm1[9]*m2j+pdn[9]*digit; |
| |
279 pdtj[20]+=pdm1[10]*m2j+pdn[10]*digit; |
| |
280 pdtj[22]+=pdm1[11]*m2j+pdn[11]*digit; |
| |
281 pdtj[24]+=pdm1[12]*m2j+pdn[12]*digit; |
| |
282 pdtj[26]+=pdm1[13]*m2j+pdn[13]*digit; |
| |
283 pdtj[28]+=pdm1[14]*m2j+pdn[14]*digit; |
| |
284 pdtj[30]+=pdm1[15]*m2j+pdn[15]*digit; |
| |
285 /* no need for cleenup, cannot overflow */ |
| |
286 digit=mod(lower32(b,Zero)*dn0,TwoToMinus16,TwoTo16); |
| |
287 } |
| |
288 } |
| |
289 |
| |
290 conv_d16_to_i32(result,dt+2*nlen,(long long *)dt,nlen+1); |
| |
291 |
| |
292 adjust_montf_result(result,nint,nlen); |
| |
293 |
| |
294 } |
| |
295 |