85191ecb0847b85ab00fbd210b9eb7deb5ba665f
[openssl.git] / crypto / bn / asm / x86_64-gcc.c
1 #ifdef __SUNPRO_C
2 # include "../bn_asm.c" /* kind of dirty hack for Sun Studio */
3 #else
4 #include <sys/types.h>
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 #define BN_ULONG unsigned long
59
60 /*
61  * "m"(a), "+m"(r)      is the way to favor DirectPath ยต-code;
62  * "g"(0)               let the compiler to decide where does it
63  *                      want to keep the value of zero;
64  */
65 #define mul_add(r,a,word,carry) do {    \
66         register BN_ULONG high,low;     \
67         asm ("mulq %3"                  \
68                 : "=a"(low),"=d"(high)  \
69                 : "a"(word),"m"(a)      \
70                 : "cc");                \
71         asm ("addq %2,%0; adcq %3,%1"   \
72                 : "+r"(carry),"+d"(high)\
73                 : "a"(low),"g"(0)       \
74                 : "cc");                \
75         asm ("addq %2,%0; adcq %3,%1"   \
76                 : "+m"(r),"+d"(high)    \
77                 : "r"(carry),"g"(0)     \
78                 : "cc");                \
79         carry=high;                     \
80         } while (0)
81
82 #define mul(r,a,word,carry) do {        \
83         register BN_ULONG high,low;     \
84         asm ("mulq %3"                  \
85                 : "=a"(low),"=d"(high)  \
86                 : "a"(word),"g"(a)      \
87                 : "cc");                \
88         asm ("addq %2,%0; adcq %3,%1"   \
89                 : "+r"(carry),"+d"(high)\
90                 : "a"(low),"g"(0)       \
91                 : "cc");                \
92         (r)=carry, carry=high;          \
93         } while (0)
94
95 #define sqr(r0,r1,a)                    \
96         asm ("mulq %2"                  \
97                 : "=a"(r0),"=d"(r1)     \
98                 : "a"(a)                \
99                 : "cc");
100
101 BN_ULONG bn_mul_add_words(BN_ULONG *rp, BN_ULONG *ap, int num, BN_ULONG w)
102         {
103         BN_ULONG c1=0;
104
105         if (num <= 0) return(c1);
106
107         while (num&~3)
108                 {
109                 mul_add(rp[0],ap[0],w,c1);
110                 mul_add(rp[1],ap[1],w,c1);
111                 mul_add(rp[2],ap[2],w,c1);
112                 mul_add(rp[3],ap[3],w,c1);
113                 ap+=4; rp+=4; num-=4;
114                 }
115         if (num)
116                 {
117                 mul_add(rp[0],ap[0],w,c1); if (--num==0) return c1;
118                 mul_add(rp[1],ap[1],w,c1); if (--num==0) return c1;
119                 mul_add(rp[2],ap[2],w,c1); return c1;
120                 }
121         
122         return(c1);
123         } 
124
125 BN_ULONG bn_mul_words(BN_ULONG *rp, BN_ULONG *ap, int num, BN_ULONG w)
126         {
127         BN_ULONG c1=0;
128
129         if (num <= 0) return(c1);
130
131         while (num&~3)
132                 {
133                 mul(rp[0],ap[0],w,c1);
134                 mul(rp[1],ap[1],w,c1);
135                 mul(rp[2],ap[2],w,c1);
136                 mul(rp[3],ap[3],w,c1);
137                 ap+=4; rp+=4; num-=4;
138                 }
139         if (num)
140                 {
141                 mul(rp[0],ap[0],w,c1); if (--num == 0) return c1;
142                 mul(rp[1],ap[1],w,c1); if (--num == 0) return c1;
143                 mul(rp[2],ap[2],w,c1);
144                 }
145         return(c1);
146         } 
147
148 void bn_sqr_words(BN_ULONG *r, BN_ULONG *a, int n)
149         {
150         if (n <= 0) return;
151
152         while (n&~3)
153                 {
154                 sqr(r[0],r[1],a[0]);
155                 sqr(r[2],r[3],a[1]);
156                 sqr(r[4],r[5],a[2]);
157                 sqr(r[6],r[7],a[3]);
158                 a+=4; r+=8; n-=4;
159                 }
160         if (n)
161                 {
162                 sqr(r[0],r[1],a[0]); if (--n == 0) return;
163                 sqr(r[2],r[3],a[1]); if (--n == 0) return;
164                 sqr(r[4],r[5],a[2]);
165                 }
166         }
167
168 BN_ULONG bn_div_words(BN_ULONG h, BN_ULONG l, BN_ULONG d)
169 {       BN_ULONG ret,waste;
170
171         asm ("divq      %4"
172                 : "=a"(ret),"=d"(waste)
173                 : "a"(l),"d"(h),"g"(d)
174                 : "cc");
175
176         return ret;
177 }
178
179 BN_ULONG bn_add_words (BN_ULONG *rp, BN_ULONG *ap, BN_ULONG *bp, size_t n)
180 { BN_ULONG ret=0,i=0;
181
182         if (n <= 0) return 0;
183
184         asm (
185         "       subq    %2,%2           \n"
186         ".align 16                      \n"
187         "1:     movq    (%4,%2,8),%0    \n"
188         "       adcq    (%5,%2,8),%0    \n"
189         "       movq    %0,(%3,%2,8)    \n"
190         "       leaq    1(%2),%2        \n"
191         "       loop    1b              \n"
192         "       sbbq    %0,%0           \n"
193                 : "=&a"(ret),"+c"(n),"=&r"(i)
194                 : "r"(rp),"r"(ap),"r"(bp)
195                 : "cc"
196         );
197
198   return ret&1;
199 }
200
201 #ifndef SIMICS
202 BN_ULONG bn_sub_words (BN_ULONG *rp, BN_ULONG *ap, BN_ULONG *bp,int n)
203 { BN_ULONG ret=0,i=0;
204
205         if (n <= 0) return 0;
206
207         asm (
208         "       subq    %2,%2           \n"
209         ".align 16                      \n"
210         "1:     movq    (%4,%2,8),%0    \n"
211         "       sbbq    (%5,%2,8),%0    \n"
212         "       movq    %0,(%3,%2,8)    \n"
213         "       leaq    1(%2),%2        \n"
214         "       loop    1b              \n"
215         "       sbbq    %0,%0           \n"
216                 : "=&a"(ret),"+c"(n),"=&r"(i)
217                 : "r"(rp),"r"(ap),"r"(bp)
218                 : "cc"
219         );
220
221   return ret&1;
222 }
223 #else
224 /* Simics 1.4<7 has buggy sbbq:-( */
225 #define BN_MASK2 0xffffffffffffffffL
226 BN_ULONG bn_sub_words(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b, int n)
227         {
228         BN_ULONG t1,t2;
229         int c=0;
230
231         if (n <= 0) return((BN_ULONG)0);
232
233         for (;;)
234                 {
235                 t1=a[0]; t2=b[0];
236                 r[0]=(t1-t2-c)&BN_MASK2;
237                 if (t1 != t2) c=(t1 < t2);
238                 if (--n <= 0) break;
239
240                 t1=a[1]; t2=b[1];
241                 r[1]=(t1-t2-c)&BN_MASK2;
242                 if (t1 != t2) c=(t1 < t2);
243                 if (--n <= 0) break;
244
245                 t1=a[2]; t2=b[2];
246                 r[2]=(t1-t2-c)&BN_MASK2;
247                 if (t1 != t2) c=(t1 < t2);
248                 if (--n <= 0) break;
249
250                 t1=a[3]; t2=b[3];
251                 r[3]=(t1-t2-c)&BN_MASK2;
252                 if (t1 != t2) c=(t1 < t2);
253                 if (--n <= 0) break;
254
255                 a+=4;
256                 b+=4;
257                 r+=4;
258                 }
259         return(c);
260         }
261 #endif
262
263 /* mul_add_c(a,b,c0,c1,c2)  -- c+=a*b for three word number c=(c2,c1,c0) */
264 /* mul_add_c2(a,b,c0,c1,c2) -- c+=2*a*b for three word number c=(c2,c1,c0) */
265 /* sqr_add_c(a,i,c0,c1,c2)  -- c+=a[i]^2 for three word number c=(c2,c1,c0) */
266 /* sqr_add_c2(a,i,c0,c1,c2) -- c+=2*a[i]*a[j] for three word number c=(c2,c1,c0) */
267
268 #if 0
269 /* original macros are kept for reference purposes */
270 #define mul_add_c(a,b,c0,c1,c2) {       \
271         BN_ULONG ta=(a),tb=(b);         \
272         t1 = ta * tb;                   \
273         t2 = BN_UMULT_HIGH(ta,tb);      \
274         c0 += t1; t2 += (c0<t1)?1:0;    \
275         c1 += t2; c2 += (c1<t2)?1:0;    \
276         }
277
278 #define mul_add_c2(a,b,c0,c1,c2) {      \
279         BN_ULONG ta=(a),tb=(b),t0;      \
280         t1 = BN_UMULT_HIGH(ta,tb);      \
281         t0 = ta * tb;                   \
282         t2 = t1+t1; c2 += (t2<t1)?1:0;  \
283         t1 = t0+t0; t2 += (t1<t0)?1:0;  \
284         c0 += t1; t2 += (c0<t1)?1:0;    \
285         c1 += t2; c2 += (c1<t2)?1:0;    \
286         }
287 #else
288 #define mul_add_c(a,b,c0,c1,c2) do {    \
289         asm ("mulq %3"                  \
290                 : "=a"(t1),"=d"(t2)     \
291                 : "a"(a),"m"(b)         \
292                 : "cc");                \
293         asm ("addq %2,%0; adcq %3,%1"   \
294                 : "+r"(c0),"+d"(t2)     \
295                 : "a"(t1),"g"(0)        \
296                 : "cc");                \
297         asm ("addq %2,%0; adcq %3,%1"   \
298                 : "+r"(c1),"+r"(c2)     \
299                 : "d"(t2),"g"(0)        \
300                 : "cc");                \
301         } while (0)
302
303 #define sqr_add_c(a,i,c0,c1,c2) do {    \
304         asm ("mulq %2"                  \
305                 : "=a"(t1),"=d"(t2)     \
306                 : "a"(a[i])             \
307                 : "cc");                \
308         asm ("addq %2,%0; adcq %3,%1"   \
309                 : "+r"(c0),"+d"(t2)     \
310                 : "a"(t1),"g"(0)        \
311                 : "cc");                \
312         asm ("addq %2,%0; adcq %3,%1"   \
313                 : "+r"(c1),"+r"(c2)     \
314                 : "d"(t2),"g"(0)        \
315                 : "cc");                \
316         } while (0)
317
318 #define mul_add_c2(a,b,c0,c1,c2) do {   \
319         asm ("mulq %3"                  \
320                 : "=a"(t1),"=d"(t2)     \
321                 : "a"(a),"m"(b)         \
322                 : "cc");                \
323         asm ("addq %0,%0; adcq %2,%1"   \
324                 : "+d"(t2),"+r"(c2)     \
325                 : "g"(0)                \
326                 : "cc");                \
327         asm ("addq %0,%0; adcq %2,%1"   \
328                 : "+a"(t1),"+d"(t2)     \
329                 : "g"(0)                \
330                 : "cc");                \
331         asm ("addq %2,%0; adcq %3,%1"   \
332                 : "+r"(c0),"+d"(t2)     \
333                 : "a"(t1),"g"(0)        \
334                 : "cc");                \
335         asm ("addq %2,%0; adcq %3,%1"   \
336                 : "+r"(c1),"+r"(c2)     \
337                 : "d"(t2),"g"(0)        \
338                 : "cc");                \
339         } while (0)
340 #endif
341
342 #define sqr_add_c2(a,i,j,c0,c1,c2)      \
343         mul_add_c2((a)[i],(a)[j],c0,c1,c2)
344
345 void bn_mul_comba8(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
346         {
347         BN_ULONG t1,t2;
348         BN_ULONG c1,c2,c3;
349
350         c1=0;
351         c2=0;
352         c3=0;
353         mul_add_c(a[0],b[0],c1,c2,c3);
354         r[0]=c1;
355         c1=0;
356         mul_add_c(a[0],b[1],c2,c3,c1);
357         mul_add_c(a[1],b[0],c2,c3,c1);
358         r[1]=c2;
359         c2=0;
360         mul_add_c(a[2],b[0],c3,c1,c2);
361         mul_add_c(a[1],b[1],c3,c1,c2);
362         mul_add_c(a[0],b[2],c3,c1,c2);
363         r[2]=c3;
364         c3=0;
365         mul_add_c(a[0],b[3],c1,c2,c3);
366         mul_add_c(a[1],b[2],c1,c2,c3);
367         mul_add_c(a[2],b[1],c1,c2,c3);
368         mul_add_c(a[3],b[0],c1,c2,c3);
369         r[3]=c1;
370         c1=0;
371         mul_add_c(a[4],b[0],c2,c3,c1);
372         mul_add_c(a[3],b[1],c2,c3,c1);
373         mul_add_c(a[2],b[2],c2,c3,c1);
374         mul_add_c(a[1],b[3],c2,c3,c1);
375         mul_add_c(a[0],b[4],c2,c3,c1);
376         r[4]=c2;
377         c2=0;
378         mul_add_c(a[0],b[5],c3,c1,c2);
379         mul_add_c(a[1],b[4],c3,c1,c2);
380         mul_add_c(a[2],b[3],c3,c1,c2);
381         mul_add_c(a[3],b[2],c3,c1,c2);
382         mul_add_c(a[4],b[1],c3,c1,c2);
383         mul_add_c(a[5],b[0],c3,c1,c2);
384         r[5]=c3;
385         c3=0;
386         mul_add_c(a[6],b[0],c1,c2,c3);
387         mul_add_c(a[5],b[1],c1,c2,c3);
388         mul_add_c(a[4],b[2],c1,c2,c3);
389         mul_add_c(a[3],b[3],c1,c2,c3);
390         mul_add_c(a[2],b[4],c1,c2,c3);
391         mul_add_c(a[1],b[5],c1,c2,c3);
392         mul_add_c(a[0],b[6],c1,c2,c3);
393         r[6]=c1;
394         c1=0;
395         mul_add_c(a[0],b[7],c2,c3,c1);
396         mul_add_c(a[1],b[6],c2,c3,c1);
397         mul_add_c(a[2],b[5],c2,c3,c1);
398         mul_add_c(a[3],b[4],c2,c3,c1);
399         mul_add_c(a[4],b[3],c2,c3,c1);
400         mul_add_c(a[5],b[2],c2,c3,c1);
401         mul_add_c(a[6],b[1],c2,c3,c1);
402         mul_add_c(a[7],b[0],c2,c3,c1);
403         r[7]=c2;
404         c2=0;
405         mul_add_c(a[7],b[1],c3,c1,c2);
406         mul_add_c(a[6],b[2],c3,c1,c2);
407         mul_add_c(a[5],b[3],c3,c1,c2);
408         mul_add_c(a[4],b[4],c3,c1,c2);
409         mul_add_c(a[3],b[5],c3,c1,c2);
410         mul_add_c(a[2],b[6],c3,c1,c2);
411         mul_add_c(a[1],b[7],c3,c1,c2);
412         r[8]=c3;
413         c3=0;
414         mul_add_c(a[2],b[7],c1,c2,c3);
415         mul_add_c(a[3],b[6],c1,c2,c3);
416         mul_add_c(a[4],b[5],c1,c2,c3);
417         mul_add_c(a[5],b[4],c1,c2,c3);
418         mul_add_c(a[6],b[3],c1,c2,c3);
419         mul_add_c(a[7],b[2],c1,c2,c3);
420         r[9]=c1;
421         c1=0;
422         mul_add_c(a[7],b[3],c2,c3,c1);
423         mul_add_c(a[6],b[4],c2,c3,c1);
424         mul_add_c(a[5],b[5],c2,c3,c1);
425         mul_add_c(a[4],b[6],c2,c3,c1);
426         mul_add_c(a[3],b[7],c2,c3,c1);
427         r[10]=c2;
428         c2=0;
429         mul_add_c(a[4],b[7],c3,c1,c2);
430         mul_add_c(a[5],b[6],c3,c1,c2);
431         mul_add_c(a[6],b[5],c3,c1,c2);
432         mul_add_c(a[7],b[4],c3,c1,c2);
433         r[11]=c3;
434         c3=0;
435         mul_add_c(a[7],b[5],c1,c2,c3);
436         mul_add_c(a[6],b[6],c1,c2,c3);
437         mul_add_c(a[5],b[7],c1,c2,c3);
438         r[12]=c1;
439         c1=0;
440         mul_add_c(a[6],b[7],c2,c3,c1);
441         mul_add_c(a[7],b[6],c2,c3,c1);
442         r[13]=c2;
443         c2=0;
444         mul_add_c(a[7],b[7],c3,c1,c2);
445         r[14]=c3;
446         r[15]=c1;
447         }
448
449 void bn_mul_comba4(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
450         {
451         BN_ULONG t1,t2;
452         BN_ULONG c1,c2,c3;
453
454         c1=0;
455         c2=0;
456         c3=0;
457         mul_add_c(a[0],b[0],c1,c2,c3);
458         r[0]=c1;
459         c1=0;
460         mul_add_c(a[0],b[1],c2,c3,c1);
461         mul_add_c(a[1],b[0],c2,c3,c1);
462         r[1]=c2;
463         c2=0;
464         mul_add_c(a[2],b[0],c3,c1,c2);
465         mul_add_c(a[1],b[1],c3,c1,c2);
466         mul_add_c(a[0],b[2],c3,c1,c2);
467         r[2]=c3;
468         c3=0;
469         mul_add_c(a[0],b[3],c1,c2,c3);
470         mul_add_c(a[1],b[2],c1,c2,c3);
471         mul_add_c(a[2],b[1],c1,c2,c3);
472         mul_add_c(a[3],b[0],c1,c2,c3);
473         r[3]=c1;
474         c1=0;
475         mul_add_c(a[3],b[1],c2,c3,c1);
476         mul_add_c(a[2],b[2],c2,c3,c1);
477         mul_add_c(a[1],b[3],c2,c3,c1);
478         r[4]=c2;
479         c2=0;
480         mul_add_c(a[2],b[3],c3,c1,c2);
481         mul_add_c(a[3],b[2],c3,c1,c2);
482         r[5]=c3;
483         c3=0;
484         mul_add_c(a[3],b[3],c1,c2,c3);
485         r[6]=c1;
486         r[7]=c2;
487         }
488
489 void bn_sqr_comba8(BN_ULONG *r, BN_ULONG *a)
490         {
491         BN_ULONG t1,t2;
492         BN_ULONG c1,c2,c3;
493
494         c1=0;
495         c2=0;
496         c3=0;
497         sqr_add_c(a,0,c1,c2,c3);
498         r[0]=c1;
499         c1=0;
500         sqr_add_c2(a,1,0,c2,c3,c1);
501         r[1]=c2;
502         c2=0;
503         sqr_add_c(a,1,c3,c1,c2);
504         sqr_add_c2(a,2,0,c3,c1,c2);
505         r[2]=c3;
506         c3=0;
507         sqr_add_c2(a,3,0,c1,c2,c3);
508         sqr_add_c2(a,2,1,c1,c2,c3);
509         r[3]=c1;
510         c1=0;
511         sqr_add_c(a,2,c2,c3,c1);
512         sqr_add_c2(a,3,1,c2,c3,c1);
513         sqr_add_c2(a,4,0,c2,c3,c1);
514         r[4]=c2;
515         c2=0;
516         sqr_add_c2(a,5,0,c3,c1,c2);
517         sqr_add_c2(a,4,1,c3,c1,c2);
518         sqr_add_c2(a,3,2,c3,c1,c2);
519         r[5]=c3;
520         c3=0;
521         sqr_add_c(a,3,c1,c2,c3);
522         sqr_add_c2(a,4,2,c1,c2,c3);
523         sqr_add_c2(a,5,1,c1,c2,c3);
524         sqr_add_c2(a,6,0,c1,c2,c3);
525         r[6]=c1;
526         c1=0;
527         sqr_add_c2(a,7,0,c2,c3,c1);
528         sqr_add_c2(a,6,1,c2,c3,c1);
529         sqr_add_c2(a,5,2,c2,c3,c1);
530         sqr_add_c2(a,4,3,c2,c3,c1);
531         r[7]=c2;
532         c2=0;
533         sqr_add_c(a,4,c3,c1,c2);
534         sqr_add_c2(a,5,3,c3,c1,c2);
535         sqr_add_c2(a,6,2,c3,c1,c2);
536         sqr_add_c2(a,7,1,c3,c1,c2);
537         r[8]=c3;
538         c3=0;
539         sqr_add_c2(a,7,2,c1,c2,c3);
540         sqr_add_c2(a,6,3,c1,c2,c3);
541         sqr_add_c2(a,5,4,c1,c2,c3);
542         r[9]=c1;
543         c1=0;
544         sqr_add_c(a,5,c2,c3,c1);
545         sqr_add_c2(a,6,4,c2,c3,c1);
546         sqr_add_c2(a,7,3,c2,c3,c1);
547         r[10]=c2;
548         c2=0;
549         sqr_add_c2(a,7,4,c3,c1,c2);
550         sqr_add_c2(a,6,5,c3,c1,c2);
551         r[11]=c3;
552         c3=0;
553         sqr_add_c(a,6,c1,c2,c3);
554         sqr_add_c2(a,7,5,c1,c2,c3);
555         r[12]=c1;
556         c1=0;
557         sqr_add_c2(a,7,6,c2,c3,c1);
558         r[13]=c2;
559         c2=0;
560         sqr_add_c(a,7,c3,c1,c2);
561         r[14]=c3;
562         r[15]=c1;
563         }
564
565 void bn_sqr_comba4(BN_ULONG *r, BN_ULONG *a)
566         {
567         BN_ULONG t1,t2;
568         BN_ULONG c1,c2,c3;
569
570         c1=0;
571         c2=0;
572         c3=0;
573         sqr_add_c(a,0,c1,c2,c3);
574         r[0]=c1;
575         c1=0;
576         sqr_add_c2(a,1,0,c2,c3,c1);
577         r[1]=c2;
578         c2=0;
579         sqr_add_c(a,1,c3,c1,c2);
580         sqr_add_c2(a,2,0,c3,c1,c2);
581         r[2]=c3;
582         c3=0;
583         sqr_add_c2(a,3,0,c1,c2,c3);
584         sqr_add_c2(a,2,1,c1,c2,c3);
585         r[3]=c1;
586         c1=0;
587         sqr_add_c(a,2,c2,c3,c1);
588         sqr_add_c2(a,3,1,c2,c3,c1);
589         r[4]=c2;
590         c2=0;
591         sqr_add_c2(a,3,2,c3,c1,c2);
592         r[5]=c3;
593         c3=0;
594         sqr_add_c(a,3,c1,c2,c3);
595         r[6]=c1;
596         r[7]=c2;
597         }
598 #endif