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 /*
 280  * Keep in mind that carrying into high part of multiplication result
 281  * can not overflow, because it cannot be all-ones.
 282  */
 283 #if 0
 284 /* original macros are kept for reference purposes */
 285 #define mul_add_c(a,b,c0,c1,c2)         do {    \
 286         BN_ULONG ta = (a), tb = (b);            \
 287         BN_ULONG lo, hi;                        \
 288         BN_UMULT_LOHI(lo,hi,ta,tb);             \
 289         c0 += lo; hi += (c0<lo)?1:0;            \
 290         c1 += hi; c2 += (c1<hi)?1:0;            \
 291         } while(0)
 292
 293 #define mul_add_c2(a,b,c0,c1,c2)        do {    \
 294         BN_ULONG ta = (a), tb = (b);            \
 295         BN_ULONG lo, hi, tt;                    \
 296         BN_UMULT_LOHI(lo,hi,ta,tb);             \
 297         c0 += lo; tt = hi+((c0<lo)?1:0);        \
 298         c1 += tt; c2 += (c1<tt)?1:0;            \
 299         c0 += lo; hi += (c0<lo)?1:0;            \
 300         c1 += hi; c2 += (c1<hi)?1:0;            \
 301         } while(0)
 302
 303 #define sqr_add_c(a,i,c0,c1,c2)         do {    \
 304         BN_ULONG ta = (a)[i];                   \
 305         BN_ULONG lo, hi;                        \
 306         BN_UMULT_LOHI(lo,hi,ta,ta);             \
 307         c0 += lo; hi += (c0<lo)?1:0;            \
 308         c1 += hi; c2 += (c1<hi)?1:0;            \
 309         } while(0)
 310 #else
 311 #define mul_add_c(a,b,c0,c1,c2) do {    \
 312         BN_ULONG t1,t2;                 \
 313         asm ("mulq %3"                  \
 314                 : "=a"(t1),"=d"(t2)     \
 315                 : "a"(a),"m"(b)         \
 316                 : "cc");                \
 317         asm ("addq %3,%0; adcq %4,%1; adcq %5,%2"       \
 318                 : "+r"(c0),"+r"(c1),"+r"(c2)            \
 319                 : "r"(t1),"r"(t2),"g"(0)                \
 320                 : "cc");                                \
 321         } while (0)
 322
 323 #define sqr_add_c(a,i,c0,c1,c2) do {    \
 324         BN_ULONG t1,t2;                 \
 325         asm ("mulq %2"                  \
 326                 : "=a"(t1),"=d"(t2)     \
 327                 : "a"(a[i])             \
 328                 : "cc");                \
 329         asm ("addq %3,%0; adcq %4,%1; adcq %5,%2"       \
 330                 : "+r"(c0),"+r"(c1),"+r"(c2)            \
 331                 : "r"(t1),"r"(t2),"g"(0)                \
 332                 : "cc");                                \
 333         } while (0)
 334
 335 #define mul_add_c2(a,b,c0,c1,c2) do {   \
 336         BN_ULONG t1,t2;                 \
 337         asm ("mulq %3"                  \
 338                 : "=a"(t1),"=d"(t2)     \
 339                 : "a"(a),"m"(b)         \
 340                 : "cc");                \
 341         asm ("addq %3,%0; adcq %4,%1; adcq %5,%2"       \
 342                 : "+r"(c0),"+r"(c1),"+r"(c2)            \
 343                 : "r"(t1),"r"(t2),"g"(0)                \
 344                 : "cc");                                \
 345         asm ("addq %3,%0; adcq %4,%1; adcq %5,%2"       \
 346                 : "+r"(c0),"+r"(c1),"+r"(c2)            \
 347                 : "r"(t1),"r"(t2),"g"(0)                \
 348                 : "cc");                                \
 349         } while (0)
 350 #endif
 351
 352 #define sqr_add_c2(a,i,j,c0,c1,c2)      \
 353         mul_add_c2((a)[i],(a)[j],c0,c1,c2)
 354
 355 void bn_mul_comba8(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
 356         {
 357         BN_ULONG c1,c2,c3;
 358
 359         c1=0;
 360         c2=0;
 361         c3=0;
 362         mul_add_c(a[0],b[0],c1,c2,c3);
 363         r[0]=c1;
 364         c1=0;
 365         mul_add_c(a[0],b[1],c2,c3,c1);
 366         mul_add_c(a[1],b[0],c2,c3,c1);
 367         r[1]=c2;
 368         c2=0;
 369         mul_add_c(a[2],b[0],c3,c1,c2);
 370         mul_add_c(a[1],b[1],c3,c1,c2);
 371         mul_add_c(a[0],b[2],c3,c1,c2);
 372         r[2]=c3;
 373         c3=0;
 374         mul_add_c(a[0],b[3],c1,c2,c3);
 375         mul_add_c(a[1],b[2],c1,c2,c3);
 376         mul_add_c(a[2],b[1],c1,c2,c3);
 377         mul_add_c(a[3],b[0],c1,c2,c3);
 378         r[3]=c1;
 379         c1=0;
 380         mul_add_c(a[4],b[0],c2,c3,c1);
 381         mul_add_c(a[3],b[1],c2,c3,c1);
 382         mul_add_c(a[2],b[2],c2,c3,c1);
 383         mul_add_c(a[1],b[3],c2,c3,c1);
 384         mul_add_c(a[0],b[4],c2,c3,c1);
 385         r[4]=c2;
 386         c2=0;
 387         mul_add_c(a[0],b[5],c3,c1,c2);
 388         mul_add_c(a[1],b[4],c3,c1,c2);
 389         mul_add_c(a[2],b[3],c3,c1,c2);
 390         mul_add_c(a[3],b[2],c3,c1,c2);
 391         mul_add_c(a[4],b[1],c3,c1,c2);
 392         mul_add_c(a[5],b[0],c3,c1,c2);
 393         r[5]=c3;
 394         c3=0;
 395         mul_add_c(a[6],b[0],c1,c2,c3);
 396         mul_add_c(a[5],b[1],c1,c2,c3);
 397         mul_add_c(a[4],b[2],c1,c2,c3);
 398         mul_add_c(a[3],b[3],c1,c2,c3);
 399         mul_add_c(a[2],b[4],c1,c2,c3);
 400         mul_add_c(a[1],b[5],c1,c2,c3);
 401         mul_add_c(a[0],b[6],c1,c2,c3);
 402         r[6]=c1;
 403         c1=0;
 404         mul_add_c(a[0],b[7],c2,c3,c1);
 405         mul_add_c(a[1],b[6],c2,c3,c1);
 406         mul_add_c(a[2],b[5],c2,c3,c1);
 407         mul_add_c(a[3],b[4],c2,c3,c1);
 408         mul_add_c(a[4],b[3],c2,c3,c1);
 409         mul_add_c(a[5],b[2],c2,c3,c1);
 410         mul_add_c(a[6],b[1],c2,c3,c1);
 411         mul_add_c(a[7],b[0],c2,c3,c1);
 412         r[7]=c2;
 413         c2=0;
 414         mul_add_c(a[7],b[1],c3,c1,c2);
 415         mul_add_c(a[6],b[2],c3,c1,c2);
 416         mul_add_c(a[5],b[3],c3,c1,c2);
 417         mul_add_c(a[4],b[4],c3,c1,c2);
 418         mul_add_c(a[3],b[5],c3,c1,c2);
 419         mul_add_c(a[2],b[6],c3,c1,c2);
 420         mul_add_c(a[1],b[7],c3,c1,c2);
 421         r[8]=c3;
 422         c3=0;
 423         mul_add_c(a[2],b[7],c1,c2,c3);
 424         mul_add_c(a[3],b[6],c1,c2,c3);
 425         mul_add_c(a[4],b[5],c1,c2,c3);
 426         mul_add_c(a[5],b[4],c1,c2,c3);
 427         mul_add_c(a[6],b[3],c1,c2,c3);
 428         mul_add_c(a[7],b[2],c1,c2,c3);
 429         r[9]=c1;
 430         c1=0;
 431         mul_add_c(a[7],b[3],c2,c3,c1);
 432         mul_add_c(a[6],b[4],c2,c3,c1);
 433         mul_add_c(a[5],b[5],c2,c3,c1);
 434         mul_add_c(a[4],b[6],c2,c3,c1);
 435         mul_add_c(a[3],b[7],c2,c3,c1);
 436         r[10]=c2;
 437         c2=0;
 438         mul_add_c(a[4],b[7],c3,c1,c2);
 439         mul_add_c(a[5],b[6],c3,c1,c2);
 440         mul_add_c(a[6],b[5],c3,c1,c2);
 441         mul_add_c(a[7],b[4],c3,c1,c2);
 442         r[11]=c3;
 443         c3=0;
 444         mul_add_c(a[7],b[5],c1,c2,c3);
 445         mul_add_c(a[6],b[6],c1,c2,c3);
 446         mul_add_c(a[5],b[7],c1,c2,c3);
 447         r[12]=c1;
 448         c1=0;
 449         mul_add_c(a[6],b[7],c2,c3,c1);
 450         mul_add_c(a[7],b[6],c2,c3,c1);
 451         r[13]=c2;
 452         c2=0;
 453         mul_add_c(a[7],b[7],c3,c1,c2);
 454         r[14]=c3;
 455         r[15]=c1;
 456         }
 457
 458 void bn_mul_comba4(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
 459         {
 460         BN_ULONG c1,c2,c3;
 461
 462         c1=0;
 463         c2=0;
 464         c3=0;
 465         mul_add_c(a[0],b[0],c1,c2,c3);
 466         r[0]=c1;
 467         c1=0;
 468         mul_add_c(a[0],b[1],c2,c3,c1);
 469         mul_add_c(a[1],b[0],c2,c3,c1);
 470         r[1]=c2;
 471         c2=0;
 472         mul_add_c(a[2],b[0],c3,c1,c2);
 473         mul_add_c(a[1],b[1],c3,c1,c2);
 474         mul_add_c(a[0],b[2],c3,c1,c2);
 475         r[2]=c3;
 476         c3=0;
 477         mul_add_c(a[0],b[3],c1,c2,c3);
 478         mul_add_c(a[1],b[2],c1,c2,c3);
 479         mul_add_c(a[2],b[1],c1,c2,c3);
 480         mul_add_c(a[3],b[0],c1,c2,c3);
 481         r[3]=c1;
 482         c1=0;
 483         mul_add_c(a[3],b[1],c2,c3,c1);
 484         mul_add_c(a[2],b[2],c2,c3,c1);
 485         mul_add_c(a[1],b[3],c2,c3,c1);
 486         r[4]=c2;
 487         c2=0;
 488         mul_add_c(a[2],b[3],c3,c1,c2);
 489         mul_add_c(a[3],b[2],c3,c1,c2);
 490         r[5]=c3;
 491         c3=0;
 492         mul_add_c(a[3],b[3],c1,c2,c3);
 493         r[6]=c1;
 494         r[7]=c2;
 495         }
 496
 497 void bn_sqr_comba8(BN_ULONG *r, const BN_ULONG *a)
 498         {
 499         BN_ULONG c1,c2,c3;
 500
 501         c1=0;
 502         c2=0;
 503         c3=0;
 504         sqr_add_c(a,0,c1,c2,c3);
 505         r[0]=c1;
 506         c1=0;
 507         sqr_add_c2(a,1,0,c2,c3,c1);
 508         r[1]=c2;
 509         c2=0;
 510         sqr_add_c(a,1,c3,c1,c2);
 511         sqr_add_c2(a,2,0,c3,c1,c2);
 512         r[2]=c3;
 513         c3=0;
 514         sqr_add_c2(a,3,0,c1,c2,c3);
 515         sqr_add_c2(a,2,1,c1,c2,c3);
 516         r[3]=c1;
 517         c1=0;
 518         sqr_add_c(a,2,c2,c3,c1);
 519         sqr_add_c2(a,3,1,c2,c3,c1);
 520         sqr_add_c2(a,4,0,c2,c3,c1);
 521         r[4]=c2;
 522         c2=0;
 523         sqr_add_c2(a,5,0,c3,c1,c2);
 524         sqr_add_c2(a,4,1,c3,c1,c2);
 525         sqr_add_c2(a,3,2,c3,c1,c2);
 526         r[5]=c3;
 527         c3=0;
 528         sqr_add_c(a,3,c1,c2,c3);
 529         sqr_add_c2(a,4,2,c1,c2,c3);
 530         sqr_add_c2(a,5,1,c1,c2,c3);
 531         sqr_add_c2(a,6,0,c1,c2,c3);
 532         r[6]=c1;
 533         c1=0;
 534         sqr_add_c2(a,7,0,c2,c3,c1);
 535         sqr_add_c2(a,6,1,c2,c3,c1);
 536         sqr_add_c2(a,5,2,c2,c3,c1);
 537         sqr_add_c2(a,4,3,c2,c3,c1);
 538         r[7]=c2;
 539         c2=0;
 540         sqr_add_c(a,4,c3,c1,c2);
 541         sqr_add_c2(a,5,3,c3,c1,c2);
 542         sqr_add_c2(a,6,2,c3,c1,c2);
 543         sqr_add_c2(a,7,1,c3,c1,c2);
 544         r[8]=c3;
 545         c3=0;
 546         sqr_add_c2(a,7,2,c1,c2,c3);
 547         sqr_add_c2(a,6,3,c1,c2,c3);
 548         sqr_add_c2(a,5,4,c1,c2,c3);
 549         r[9]=c1;
 550         c1=0;
 551         sqr_add_c(a,5,c2,c3,c1);
 552         sqr_add_c2(a,6,4,c2,c3,c1);
 553         sqr_add_c2(a,7,3,c2,c3,c1);
 554         r[10]=c2;
 555         c2=0;
 556         sqr_add_c2(a,7,4,c3,c1,c2);
 557         sqr_add_c2(a,6,5,c3,c1,c2);
 558         r[11]=c3;
 559         c3=0;
 560         sqr_add_c(a,6,c1,c2,c3);
 561         sqr_add_c2(a,7,5,c1,c2,c3);
 562         r[12]=c1;
 563         c1=0;
 564         sqr_add_c2(a,7,6,c2,c3,c1);
 565         r[13]=c2;
 566         c2=0;
 567         sqr_add_c(a,7,c3,c1,c2);
 568         r[14]=c3;
 569         r[15]=c1;
 570         }
 571
 572 void bn_sqr_comba4(BN_ULONG *r, const BN_ULONG *a)
 573         {
 574         BN_ULONG c1,c2,c3;
 575
 576         c1=0;
 577         c2=0;
 578         c3=0;
 579         sqr_add_c(a,0,c1,c2,c3);
 580         r[0]=c1;
 581         c1=0;
 582         sqr_add_c2(a,1,0,c2,c3,c1);
 583         r[1]=c2;
 584         c2=0;
 585         sqr_add_c(a,1,c3,c1,c2);
 586         sqr_add_c2(a,2,0,c3,c1,c2);
 587         r[2]=c3;
 588         c3=0;
 589         sqr_add_c2(a,3,0,c1,c2,c3);
 590         sqr_add_c2(a,2,1,c1,c2,c3);
 591         r[3]=c1;
 592         c1=0;
 593         sqr_add_c(a,2,c2,c3,c1);
 594         sqr_add_c2(a,3,1,c2,c3,c1);
 595         r[4]=c2;
 596         c2=0;
 597         sqr_add_c2(a,3,2,c3,c1,c2);
 598         r[5]=c3;
 599         c3=0;
 600         sqr_add_c(a,3,c1,c2,c3);
 601         r[6]=c1;
 602         r[7]=c2;
 603         }
 604 #endif