crypto/openssl-0.9/crypto/bn/asm/x86_64-gcc.c

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