c6d12f4f81a1613baa1a3199d9eaf62eddf966df
[openssl.git] / crypto / bn / asm / x86_64-gcc.c
1 #include "../bn_lcl.h"
2 #if !(defined(__GNUC__) && __GNUC__>=2)
3 # include "../bn_asm.c" /* kind of dirty hack for Sun Studio */
4 #else
5 /*-
6  * x86_64 BIGNUM accelerator version 0.1, December 2002.
7  *
8  * Implemented by Andy Polyakov <appro@fy.chalmers.se> for the OpenSSL
9  * project.
10  *
11  * Rights for redistribution and usage in source and binary forms are
12  * granted according to the OpenSSL license. Warranty of any kind is
13  * disclaimed.
14  *
15  * Q. Version 0.1? It doesn't sound like Andy, he used to assign real
16  *    versions, like 1.0...
17  * A. Well, that's because this code is basically a quick-n-dirty
18  *    proof-of-concept hack. As you can see it's implemented with
19  *    inline assembler, which means that you're bound to GCC and that
20  *    there might be enough room for further improvement.
21  *
22  * Q. Why inline assembler?
23  * A. x86_64 features own ABI which I'm not familiar with. This is
24  *    why I decided to let the compiler take care of subroutine
25  *    prologue/epilogue as well as register allocation. For reference.
26  *    Win64 implements different ABI for AMD64, different from Linux.
27  *
28  * Q. How much faster does it get?
29  * A. 'apps/openssl speed rsa dsa' output with no-asm:
30  *
31  *                        sign    verify    sign/s verify/s
32  *      rsa  512 bits   0.0006s   0.0001s   1683.8  18456.2
33  *      rsa 1024 bits   0.0028s   0.0002s    356.0   6407.0
34  *      rsa 2048 bits   0.0172s   0.0005s     58.0   1957.8
35  *      rsa 4096 bits   0.1155s   0.0018s      8.7    555.6
36  *                        sign    verify    sign/s verify/s
37  *      dsa  512 bits   0.0005s   0.0006s   2100.8   1768.3
38  *      dsa 1024 bits   0.0014s   0.0018s    692.3    559.2
39  *      dsa 2048 bits   0.0049s   0.0061s    204.7    165.0
40  *
41  *    'apps/openssl speed rsa dsa' output with this module:
42  *
43  *                        sign    verify    sign/s verify/s
44  *      rsa  512 bits   0.0004s   0.0000s   2767.1  33297.9
45  *      rsa 1024 bits   0.0012s   0.0001s    867.4  14674.7
46  *      rsa 2048 bits   0.0061s   0.0002s    164.0   5270.0
47  *      rsa 4096 bits   0.0384s   0.0006s     26.1   1650.8
48  *                        sign    verify    sign/s verify/s
49  *      dsa  512 bits   0.0002s   0.0003s   4442.2   3786.3
50  *      dsa 1024 bits   0.0005s   0.0007s   1835.1   1497.4
51  *      dsa 2048 bits   0.0016s   0.0020s    620.4    504.6
52  *
53  *    For the reference. IA-32 assembler implementation performs
54  *    very much like 64-bit code compiled with no-asm on the same
55  *    machine.
56  */
57
58 #if defined(_WIN64) || !defined(__LP64__)
59 #define BN_ULONG unsigned long long
60 #else
61 #define BN_ULONG unsigned long
62 #endif
63
64 #undef mul
65 #undef mul_add
66
67 /*-
68  * "m"(a), "+m"(r)      is the way to favor DirectPath ยต-code;
69  * "g"(0)               let the compiler to decide where does it
70  *                      want to keep the value of zero;
71  */
72 #define mul_add(r,a,word,carry) do {    \
73         register BN_ULONG high,low;     \
74         asm ("mulq %3"                  \
75                 : "=a"(low),"=d"(high)  \
76                 : "a"(word),"m"(a)      \
77                 : "cc");                \
78         asm ("addq %2,%0; adcq %3,%1"   \
79                 : "+r"(carry),"+d"(high)\
80                 : "a"(low),"g"(0)       \
81                 : "cc");                \
82         asm ("addq %2,%0; adcq %3,%1"   \
83                 : "+m"(r),"+d"(high)    \
84                 : "r"(carry),"g"(0)     \
85                 : "cc");                \
86         carry=high;                     \
87         } while (0)
88
89 #define mul(r,a,word,carry) do {        \
90         register BN_ULONG high,low;     \
91         asm ("mulq %3"                  \
92                 : "=a"(low),"=d"(high)  \
93                 : "a"(word),"g"(a)      \
94                 : "cc");                \
95         asm ("addq %2,%0; adcq %3,%1"   \
96                 : "+r"(carry),"+d"(high)\
97                 : "a"(low),"g"(0)       \
98                 : "cc");                \
99         (r)=carry, carry=high;          \
100         } while (0)
101 #undef sqr
102 #define sqr(r0,r1,a)                    \
103         asm ("mulq %2"                  \
104                 : "=a"(r0),"=d"(r1)     \
105                 : "a"(a)                \
106                 : "cc");
107
108 BN_ULONG bn_mul_add_words(BN_ULONG *rp, const BN_ULONG *ap, int num, BN_ULONG w)
109         {
110         BN_ULONG c1=0;
111
112         if (num <= 0) return(c1);
113
114         while (num&~3)
115                 {
116                 mul_add(rp[0],ap[0],w,c1);
117                 mul_add(rp[1],ap[1],w,c1);
118                 mul_add(rp[2],ap[2],w,c1);
119                 mul_add(rp[3],ap[3],w,c1);
120                 ap+=4; rp+=4; num-=4;
121                 }
122         if (num)
123                 {
124                 mul_add(rp[0],ap[0],w,c1); if (--num==0) return c1;
125                 mul_add(rp[1],ap[1],w,c1); if (--num==0) return c1;
126                 mul_add(rp[2],ap[2],w,c1); return c1;
127                 }
128         
129         return(c1);
130         } 
131
132 BN_ULONG bn_mul_words(BN_ULONG *rp, const BN_ULONG *ap, int num, BN_ULONG w)
133         {
134         BN_ULONG c1=0;
135
136         if (num <= 0) return(c1);
137
138         while (num&~3)
139                 {
140                 mul(rp[0],ap[0],w,c1);
141                 mul(rp[1],ap[1],w,c1);
142                 mul(rp[2],ap[2],w,c1);
143                 mul(rp[3],ap[3],w,c1);
144                 ap+=4; rp+=4; num-=4;
145                 }
146         if (num)
147                 {
148                 mul(rp[0],ap[0],w,c1); if (--num == 0) return c1;
149                 mul(rp[1],ap[1],w,c1); if (--num == 0) return c1;
150                 mul(rp[2],ap[2],w,c1);
151                 }
152         return(c1);
153         } 
154
155 void bn_sqr_words(BN_ULONG *r, const BN_ULONG *a, int n)
156         {
157         if (n <= 0) return;
158
159         while (n&~3)
160                 {
161                 sqr(r[0],r[1],a[0]);
162                 sqr(r[2],r[3],a[1]);
163                 sqr(r[4],r[5],a[2]);
164                 sqr(r[6],r[7],a[3]);
165                 a+=4; r+=8; n-=4;
166                 }
167         if (n)
168                 {
169                 sqr(r[0],r[1],a[0]); if (--n == 0) return;
170                 sqr(r[2],r[3],a[1]); if (--n == 0) return;
171                 sqr(r[4],r[5],a[2]);
172                 }
173         }
174
175 BN_ULONG bn_div_words(BN_ULONG h, BN_ULONG l, BN_ULONG d)
176 {       BN_ULONG ret,waste;
177
178         asm ("divq      %4"
179                 : "=a"(ret),"=d"(waste)
180                 : "a"(l),"d"(h),"g"(d)
181                 : "cc");
182
183         return ret;
184 }
185
186 BN_ULONG bn_add_words (BN_ULONG *rp, const BN_ULONG *ap, const BN_ULONG *bp,int n)
187 { BN_ULONG ret;
188   size_t i=0;
189
190         if (n <= 0) return 0;
191
192         asm volatile (
193         "       subq    %0,%0           \n"     /* clear carry */
194         "       jmp     1f              \n"
195         ".p2align 4                     \n"
196         "1:     movq    (%4,%2,8),%0    \n"
197         "       adcq    (%5,%2,8),%0    \n"
198         "       movq    %0,(%3,%2,8)    \n"
199         "       lea     1(%2),%2        \n"
200         "       loop    1b              \n"
201         "       sbbq    %0,%0           \n"
202                 : "=&r"(ret),"+c"(n),"+r"(i)
203                 : "r"(rp),"r"(ap),"r"(bp)
204                 : "cc", "memory"
205         );
206
207   return ret&1;
208 }
209
210 #ifndef SIMICS
211 BN_ULONG bn_sub_words (BN_ULONG *rp, const BN_ULONG *ap, const BN_ULONG *bp,int n)
212 { BN_ULONG ret;
213   size_t i=0;
214
215         if (n <= 0) return 0;
216
217         asm volatile (
218         "       subq    %0,%0           \n"     /* clear borrow */
219         "       jmp     1f              \n"
220         ".p2align 4                     \n"
221         "1:     movq    (%4,%2,8),%0    \n"
222         "       sbbq    (%5,%2,8),%0    \n"
223         "       movq    %0,(%3,%2,8)    \n"
224         "       lea     1(%2),%2        \n"
225         "       loop    1b              \n"
226         "       sbbq    %0,%0           \n"
227                 : "=&r"(ret),"+c"(n),"+r"(i)
228                 : "r"(rp),"r"(ap),"r"(bp)
229                 : "cc", "memory"
230         );
231
232   return ret&1;
233 }
234 #else
235 /* Simics 1.4<7 has buggy sbbq:-( */
236 #define BN_MASK2 0xffffffffffffffffL
237 BN_ULONG bn_sub_words(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n)
238         {
239         BN_ULONG t1,t2;
240         int c=0;
241
242         if (n <= 0) return((BN_ULONG)0);
243
244         for (;;)
245                 {
246                 t1=a[0]; t2=b[0];
247                 r[0]=(t1-t2-c)&BN_MASK2;
248                 if (t1 != t2) c=(t1 < t2);
249                 if (--n <= 0) break;
250
251                 t1=a[1]; t2=b[1];
252                 r[1]=(t1-t2-c)&BN_MASK2;
253                 if (t1 != t2) c=(t1 < t2);
254                 if (--n <= 0) break;
255
256                 t1=a[2]; t2=b[2];
257                 r[2]=(t1-t2-c)&BN_MASK2;
258                 if (t1 != t2) c=(t1 < t2);
259                 if (--n <= 0) break;
260
261                 t1=a[3]; t2=b[3];
262                 r[3]=(t1-t2-c)&BN_MASK2;
263                 if (t1 != t2) c=(t1 < t2);
264                 if (--n <= 0) break;
265
266                 a+=4;
267                 b+=4;
268                 r+=4;
269                 }
270         return(c);
271         }
272 #endif
273
274 /* mul_add_c(a,b,c0,c1,c2)  -- c+=a*b for three word number c=(c2,c1,c0) */
275 /* mul_add_c2(a,b,c0,c1,c2) -- c+=2*a*b for three word number c=(c2,c1,c0) */
276 /* sqr_add_c(a,i,c0,c1,c2)  -- c+=a[i]^2 for three word number c=(c2,c1,c0) */
277 /* sqr_add_c2(a,i,c0,c1,c2) -- c+=2*a[i]*a[j] for three word number c=(c2,c1,c0) */
278
279 #if 0
280 /* original macros are kept for reference purposes */
281 #define mul_add_c(a,b,c0,c1,c2) {       \
282         BN_ULONG ta=(a),tb=(b);         \
283         t1 = ta * tb;                   \
284         t2 = BN_UMULT_HIGH(ta,tb);      \
285         c0 += t1; t2 += (c0<t1)?1:0;    \
286         c1 += t2; c2 += (c1<t2)?1:0;    \
287         }
288
289 #define mul_add_c2(a,b,c0,c1,c2) {      \
290         BN_ULONG ta=(a),tb=(b),t0;      \
291         t1 = BN_UMULT_HIGH(ta,tb);      \
292         t0 = ta * tb;                   \
293         t2 = t1+t1; c2 += (t2<t1)?1:0;  \
294         t1 = t0+t0; t2 += (t1<t0)?1:0;  \
295         c0 += t1; t2 += (c0<t1)?1:0;    \
296         c1 += t2; c2 += (c1<t2)?1:0;    \
297         }
298 #else
299 #define mul_add_c(a,b,c0,c1,c2) do {    \
300         asm ("mulq %3"                  \
301                 : "=a"(t1),"=d"(t2)     \
302                 : "a"(a),"m"(b)         \
303                 : "cc");                \
304         asm ("addq %2,%0; adcq %3,%1"   \
305                 : "+r"(c0),"+d"(t2)     \
306                 : "a"(t1),"g"(0)        \
307                 : "cc");                \
308         asm ("addq %2,%0; adcq %3,%1"   \
309                 : "+r"(c1),"+r"(c2)     \
310                 : "d"(t2),"g"(0)        \
311                 : "cc");                \
312         } while (0)
313
314 #define sqr_add_c(a,i,c0,c1,c2) do {    \
315         asm ("mulq %2"                  \
316                 : "=a"(t1),"=d"(t2)     \
317                 : "a"(a[i])             \
318                 : "cc");                \
319         asm ("addq %2,%0; adcq %3,%1"   \
320                 : "+r"(c0),"+d"(t2)     \
321                 : "a"(t1),"g"(0)        \
322                 : "cc");                \
323         asm ("addq %2,%0; adcq %3,%1"   \
324                 : "+r"(c1),"+r"(c2)     \
325                 : "d"(t2),"g"(0)        \
326                 : "cc");                \
327         } while (0)
328
329 #define mul_add_c2(a,b,c0,c1,c2) do {   \
330         asm ("mulq %3"                  \
331                 : "=a"(t1),"=d"(t2)     \
332                 : "a"(a),"m"(b)         \
333                 : "cc");                \
334         asm ("addq %0,%0; adcq %2,%1"   \
335                 : "+d"(t2),"+r"(c2)     \
336                 : "g"(0)                \
337                 : "cc");                \
338         asm ("addq %0,%0; adcq %2,%1"   \
339                 : "+a"(t1),"+d"(t2)     \
340                 : "g"(0)                \
341                 : "cc");                \
342         asm ("addq %2,%0; adcq %3,%1"   \
343                 : "+r"(c0),"+d"(t2)     \
344                 : "a"(t1),"g"(0)        \
345                 : "cc");                \
346         asm ("addq %2,%0; adcq %3,%1"   \
347                 : "+r"(c1),"+r"(c2)     \
348                 : "d"(t2),"g"(0)        \
349                 : "cc");                \
350         } while (0)
351 #endif
352
353 #define sqr_add_c2(a,i,j,c0,c1,c2)      \
354         mul_add_c2((a)[i],(a)[j],c0,c1,c2)
355
356 void bn_mul_comba8(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
357         {
358         BN_ULONG t1,t2;
359         BN_ULONG c1,c2,c3;
360
361         c1=0;
362         c2=0;
363         c3=0;
364         mul_add_c(a[0],b[0],c1,c2,c3);
365         r[0]=c1;
366         c1=0;
367         mul_add_c(a[0],b[1],c2,c3,c1);
368         mul_add_c(a[1],b[0],c2,c3,c1);
369         r[1]=c2;
370         c2=0;
371         mul_add_c(a[2],b[0],c3,c1,c2);
372         mul_add_c(a[1],b[1],c3,c1,c2);
373         mul_add_c(a[0],b[2],c3,c1,c2);
374         r[2]=c3;
375         c3=0;
376         mul_add_c(a[0],b[3],c1,c2,c3);
377         mul_add_c(a[1],b[2],c1,c2,c3);
378         mul_add_c(a[2],b[1],c1,c2,c3);
379         mul_add_c(a[3],b[0],c1,c2,c3);
380         r[3]=c1;
381         c1=0;
382         mul_add_c(a[4],b[0],c2,c3,c1);
383         mul_add_c(a[3],b[1],c2,c3,c1);
384         mul_add_c(a[2],b[2],c2,c3,c1);
385         mul_add_c(a[1],b[3],c2,c3,c1);
386         mul_add_c(a[0],b[4],c2,c3,c1);
387         r[4]=c2;
388         c2=0;
389         mul_add_c(a[0],b[5],c3,c1,c2);
390         mul_add_c(a[1],b[4],c3,c1,c2);
391         mul_add_c(a[2],b[3],c3,c1,c2);
392         mul_add_c(a[3],b[2],c3,c1,c2);
393         mul_add_c(a[4],b[1],c3,c1,c2);
394         mul_add_c(a[5],b[0],c3,c1,c2);
395         r[5]=c3;
396         c3=0;
397         mul_add_c(a[6],b[0],c1,c2,c3);
398         mul_add_c(a[5],b[1],c1,c2,c3);
399         mul_add_c(a[4],b[2],c1,c2,c3);
400         mul_add_c(a[3],b[3],c1,c2,c3);
401         mul_add_c(a[2],b[4],c1,c2,c3);
402         mul_add_c(a[1],b[5],c1,c2,c3);
403         mul_add_c(a[0],b[6],c1,c2,c3);
404         r[6]=c1;
405         c1=0;
406         mul_add_c(a[0],b[7],c2,c3,c1);
407         mul_add_c(a[1],b[6],c2,c3,c1);
408         mul_add_c(a[2],b[5],c2,c3,c1);
409         mul_add_c(a[3],b[4],c2,c3,c1);
410         mul_add_c(a[4],b[3],c2,c3,c1);
411         mul_add_c(a[5],b[2],c2,c3,c1);
412         mul_add_c(a[6],b[1],c2,c3,c1);
413         mul_add_c(a[7],b[0],c2,c3,c1);
414         r[7]=c2;
415         c2=0;
416         mul_add_c(a[7],b[1],c3,c1,c2);
417         mul_add_c(a[6],b[2],c3,c1,c2);
418         mul_add_c(a[5],b[3],c3,c1,c2);
419         mul_add_c(a[4],b[4],c3,c1,c2);
420         mul_add_c(a[3],b[5],c3,c1,c2);
421         mul_add_c(a[2],b[6],c3,c1,c2);
422         mul_add_c(a[1],b[7],c3,c1,c2);
423         r[8]=c3;
424         c3=0;
425         mul_add_c(a[2],b[7],c1,c2,c3);
426         mul_add_c(a[3],b[6],c1,c2,c3);
427         mul_add_c(a[4],b[5],c1,c2,c3);
428         mul_add_c(a[5],b[4],c1,c2,c3);
429         mul_add_c(a[6],b[3],c1,c2,c3);
430         mul_add_c(a[7],b[2],c1,c2,c3);
431         r[9]=c1;
432         c1=0;
433         mul_add_c(a[7],b[3],c2,c3,c1);
434         mul_add_c(a[6],b[4],c2,c3,c1);
435         mul_add_c(a[5],b[5],c2,c3,c1);
436         mul_add_c(a[4],b[6],c2,c3,c1);
437         mul_add_c(a[3],b[7],c2,c3,c1);
438         r[10]=c2;
439         c2=0;
440         mul_add_c(a[4],b[7],c3,c1,c2);
441         mul_add_c(a[5],b[6],c3,c1,c2);
442         mul_add_c(a[6],b[5],c3,c1,c2);
443         mul_add_c(a[7],b[4],c3,c1,c2);
444         r[11]=c3;
445         c3=0;
446         mul_add_c(a[7],b[5],c1,c2,c3);
447         mul_add_c(a[6],b[6],c1,c2,c3);
448         mul_add_c(a[5],b[7],c1,c2,c3);
449         r[12]=c1;
450         c1=0;
451         mul_add_c(a[6],b[7],c2,c3,c1);
452         mul_add_c(a[7],b[6],c2,c3,c1);
453         r[13]=c2;
454         c2=0;
455         mul_add_c(a[7],b[7],c3,c1,c2);
456         r[14]=c3;
457         r[15]=c1;
458         }
459
460 void bn_mul_comba4(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
461         {
462         BN_ULONG t1,t2;
463         BN_ULONG c1,c2,c3;
464
465         c1=0;
466         c2=0;
467         c3=0;
468         mul_add_c(a[0],b[0],c1,c2,c3);
469         r[0]=c1;
470         c1=0;
471         mul_add_c(a[0],b[1],c2,c3,c1);
472         mul_add_c(a[1],b[0],c2,c3,c1);
473         r[1]=c2;
474         c2=0;
475         mul_add_c(a[2],b[0],c3,c1,c2);
476         mul_add_c(a[1],b[1],c3,c1,c2);
477         mul_add_c(a[0],b[2],c3,c1,c2);
478         r[2]=c3;
479         c3=0;
480         mul_add_c(a[0],b[3],c1,c2,c3);
481         mul_add_c(a[1],b[2],c1,c2,c3);
482         mul_add_c(a[2],b[1],c1,c2,c3);
483         mul_add_c(a[3],b[0],c1,c2,c3);
484         r[3]=c1;
485         c1=0;
486         mul_add_c(a[3],b[1],c2,c3,c1);
487         mul_add_c(a[2],b[2],c2,c3,c1);
488         mul_add_c(a[1],b[3],c2,c3,c1);
489         r[4]=c2;
490         c2=0;
491         mul_add_c(a[2],b[3],c3,c1,c2);
492         mul_add_c(a[3],b[2],c3,c1,c2);
493         r[5]=c3;
494         c3=0;
495         mul_add_c(a[3],b[3],c1,c2,c3);
496         r[6]=c1;
497         r[7]=c2;
498         }
499
500 void bn_sqr_comba8(BN_ULONG *r, const BN_ULONG *a)
501         {
502         BN_ULONG t1,t2;
503         BN_ULONG c1,c2,c3;
504
505         c1=0;
506         c2=0;
507         c3=0;
508         sqr_add_c(a,0,c1,c2,c3);
509         r[0]=c1;
510         c1=0;
511         sqr_add_c2(a,1,0,c2,c3,c1);
512         r[1]=c2;
513         c2=0;
514         sqr_add_c(a,1,c3,c1,c2);
515         sqr_add_c2(a,2,0,c3,c1,c2);
516         r[2]=c3;
517         c3=0;
518         sqr_add_c2(a,3,0,c1,c2,c3);
519         sqr_add_c2(a,2,1,c1,c2,c3);
520         r[3]=c1;
521         c1=0;
522         sqr_add_c(a,2,c2,c3,c1);
523         sqr_add_c2(a,3,1,c2,c3,c1);
524         sqr_add_c2(a,4,0,c2,c3,c1);
525         r[4]=c2;
526         c2=0;
527         sqr_add_c2(a,5,0,c3,c1,c2);
528         sqr_add_c2(a,4,1,c3,c1,c2);
529         sqr_add_c2(a,3,2,c3,c1,c2);
530         r[5]=c3;
531         c3=0;
532         sqr_add_c(a,3,c1,c2,c3);
533         sqr_add_c2(a,4,2,c1,c2,c3);
534         sqr_add_c2(a,5,1,c1,c2,c3);
535         sqr_add_c2(a,6,0,c1,c2,c3);
536         r[6]=c1;
537         c1=0;
538         sqr_add_c2(a,7,0,c2,c3,c1);
539         sqr_add_c2(a,6,1,c2,c3,c1);
540         sqr_add_c2(a,5,2,c2,c3,c1);
541         sqr_add_c2(a,4,3,c2,c3,c1);
542         r[7]=c2;
543         c2=0;
544         sqr_add_c(a,4,c3,c1,c2);
545         sqr_add_c2(a,5,3,c3,c1,c2);
546         sqr_add_c2(a,6,2,c3,c1,c2);
547         sqr_add_c2(a,7,1,c3,c1,c2);
548         r[8]=c3;
549         c3=0;
550         sqr_add_c2(a,7,2,c1,c2,c3);
551         sqr_add_c2(a,6,3,c1,c2,c3);
552         sqr_add_c2(a,5,4,c1,c2,c3);
553         r[9]=c1;
554         c1=0;
555         sqr_add_c(a,5,c2,c3,c1);
556         sqr_add_c2(a,6,4,c2,c3,c1);
557         sqr_add_c2(a,7,3,c2,c3,c1);
558         r[10]=c2;
559         c2=0;
560         sqr_add_c2(a,7,4,c3,c1,c2);
561         sqr_add_c2(a,6,5,c3,c1,c2);
562         r[11]=c3;
563         c3=0;
564         sqr_add_c(a,6,c1,c2,c3);
565         sqr_add_c2(a,7,5,c1,c2,c3);
566         r[12]=c1;
567         c1=0;
568         sqr_add_c2(a,7,6,c2,c3,c1);
569         r[13]=c2;
570         c2=0;
571         sqr_add_c(a,7,c3,c1,c2);
572         r[14]=c3;
573         r[15]=c1;
574         }
575
576 void bn_sqr_comba4(BN_ULONG *r, const BN_ULONG *a)
577         {
578         BN_ULONG t1,t2;
579         BN_ULONG c1,c2,c3;
580
581         c1=0;
582         c2=0;
583         c3=0;
584         sqr_add_c(a,0,c1,c2,c3);
585         r[0]=c1;
586         c1=0;
587         sqr_add_c2(a,1,0,c2,c3,c1);
588         r[1]=c2;
589         c2=0;
590         sqr_add_c(a,1,c3,c1,c2);
591         sqr_add_c2(a,2,0,c3,c1,c2);
592         r[2]=c3;
593         c3=0;
594         sqr_add_c2(a,3,0,c1,c2,c3);
595         sqr_add_c2(a,2,1,c1,c2,c3);
596         r[3]=c1;
597         c1=0;
598         sqr_add_c(a,2,c2,c3,c1);
599         sqr_add_c2(a,3,1,c2,c3,c1);
600         r[4]=c2;
601         c2=0;
602         sqr_add_c2(a,3,2,c3,c1,c2);
603         r[5]=c3;
604         c3=0;
605         sqr_add_c(a,3,c1,c2,c3);
606         r[6]=c1;
607         r[7]=c2;
608         }
609 #endif