- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# integer.cpp

    00001 // integer.cpp - written and placed in the public domain by Wei Dai
    00002 // contains public domain code contributed by Alister Lee and Leonard Janke
    00003
    00004 #include "pch.h"
    00005
    00006 #ifndef CRYPTOPP_IMPORTS
    00007
    00008 #include "integer.h"
    00009 #include "modarith.h"
    00010 #include "nbtheory.h"
    00011 #include "asn.h"
    00012 #include "oids.h"
    00013 #include "words.h"
    00014 #include "algparam.h"
    00015 #include "pubkey.h"             // for P1363_KDF2
    00016 #include "sha.h"
    00017 #include "cpu.h"
    00018
    00019 #include <iostream>
    00020
    00021 #if _MSC_VER >= 1400
    00022         #include <intrin.h>
    00023 #endif
    00024
    00025 #ifdef __DECCXX
    00026         #include <c_asm.h>
    00027 #endif
    00028
    00029 #ifdef CRYPTOPP_MSVC6_NO_PP
    00030         #pragma message("You do not seem to have the Visual C++ Processor Pack installed, so use of SSE2 instructions will be disabled.")
    00031 #endif
    00032
    00033 #define CRYPTOPP_INTEGER_SSE2 (CRYPTOPP_BOOL_SSE2_ASM_AVAILABLE && CRYPTOPP_BOOL_X86)
    00034
    00035 NAMESPACE_BEGIN(CryptoPP)
    00036
    00037 bool AssignIntToInteger(const std::type_info &valueType, void *pInteger, const void *pInt)
    00038 {
    00039         if (valueType != typeid(Integer))
    00040                 return false;
    00041         *reinterpret_cast<Integer *>(pInteger) = *reinterpret_cast<const int *>(pInt);
    00042         return true;
    00043 }
    00044
    00045 inline static int Compare(const word *A, const word *B, size_t N)
    00046 {
    00047         while (N--)
    00048                 if (A[N] > B[N])
    00049                         return 1;
    00050                 else if (A[N] < B[N])
    00051                         return -1;
    00052
    00053         return 0;
    00054 }
    00055
    00056 inline static int Increment(word *A, size_t N, word B=1)
    00057 {
    00058         assert(N);
    00059         word t = A[0];
    00060         A[0] = t+B;
    00061         if (A[0] >= t)
    00062                 return 0;
    00063         for (unsigned i=1; i<N; i++)
    00064                 if (++A[i])
    00065                         return 0;
    00066         return 1;
    00067 }
    00068
    00069 inline static int Decrement(word *A, size_t N, word B=1)
    00070 {
    00071         assert(N);
    00072         word t = A[0];
    00073         A[0] = t-B;
    00074         if (A[0] <= t)
    00075                 return 0;
    00076         for (unsigned i=1; i<N; i++)
    00077                 if (A[i]--)
    00078                         return 0;
    00079         return 1;
    00080 }
    00081
    00082 static void TwosComplement(word *A, size_t N)
    00083 {
    00084         Decrement(A, N);
    00085         for (unsigned i=0; i<N; i++)
    00086                 A[i] = ~A[i];
    00087 }
    00088
    00089 static word AtomicInverseModPower2(word A)
    00090 {
    00091         assert(A%2==1);
    00092
    00093         word R=A%8;
    00094
    00095         for (unsigned i=3; i<WORD_BITS; i*=2)
    00096                 R = R*(2-R*A);
    00097
    00098         assert(R*A==1);
    00099         return R;
    00100 }
    00101
    00102 // ********************************************************
    00103
    00104 #if !defined(CRYPTOPP_NATIVE_DWORD_AVAILABLE) || CRYPTOPP_BOOL_X64
    00105         #define Declare2Words(x)                        word x##0, x##1;
    00106         #define AssignWord(a, b)                        a##0 = b; a##1 = 0;
    00107         #define Add2WordsBy1(a, b, c)           a##0 = b##0 + c; a##1 = b##1 + (a##0 < c);
    00108         #define LowWord(a)                                      a##0
    00109         #define HighWord(a)                                     a##1
    00110         #ifdef _MSC_VER
    00111                 #define MultiplyWords(p, a, b)  p##0 = _umul128(a, b, &p##1);
    00112                 #define Double3Words(c, d)              d##1 = __shiftleft128(d##0, d##1, 1); d##0 = __shiftleft128(c, d##0, 1); c *= 2;
    00113         #elif defined(__DECCXX)
    00114                 #define MultiplyWords(p, a, b)  p##0 = a*b; p##1 = asm("umulh %a0, %a1, %v0", a, b);
    00115         #elif CRYPTOPP_BOOL_X64
    00116                 #define MultiplyWords(p, a, b)  asm ("mulq %3" : "=a"(p##0), "=d"(p##1) : "a"(a), "g"(b) : "cc");
    00117                 #define MulAcc(c, d, a, b)              asm ("mulq %6; addq %3, %0; adcq %4, %1; adcq $0, %2;" : "+r"(c), "+r"(d##0), "+r"(d##1), "=a"(p0), "=d"(p1) : "a"(a), "g"(b) : "cc");
    00118                 #define Double3Words(c, d)              asm ("addq %0, %0; adcq %1, %1; adcq %2, %2;" : "+r"(c), "+r"(d##0), "+r"(d##1) : : "cc");
    00119                 #define Acc2WordsBy1(a, b)              asm ("addq %2, %0; adcq $0, %1;" : "+r"(a##0), "+r"(a##1) : "r"(b) : "cc");
    00120                 #define Acc2WordsBy2(a, b)              asm ("addq %2, %0; adcq %3, %1;" : "+r"(a##0), "+r"(a##1) : "r"(b##0), "r"(b##1) : "cc");
    00121                 #define Acc3WordsBy2(c, d, e)   asm ("addq %5, %0; adcq %6, %1; adcq $0, %2;" : "+r"(c), "=r"(e##0), "=r"(e##1) : "1"(d##0), "2"(d##1), "r"(e##0), "r"(e##1) : "cc");
    00122         #endif
    00123         #ifndef Double3Words
    00124                 #define Double3Words(c, d)              d##1 = 2*d##1 + (d##0>>(WORD_BITS-1)); d##0 = 2*d##0 + (c>>(WORD_BITS-1)); c *= 2;
    00125         #endif
    00126         #ifndef Acc2WordsBy2
    00127                 #define Acc2WordsBy2(a, b)              a##0 += b##0; a##1 += a##0 < b##0; a##1 += b##1;
    00128         #endif
    00129         #define AddWithCarry(u, a, b)           {word t = a+b; u##0 = t + u##1; u##1 = (t<a) + (u##0<t);}
    00130         #define SubtractWithBorrow(u, a, b)     {word t = a-b; u##0 = t - u##1; u##1 = (t>a) + (u##0>t);}
    00131         #define GetCarry(u)                                     u##1
    00132         #define GetBorrow(u)                            u##1
    00133 #else
    00134         #define Declare2Words(x)                        dword x;
    00135         #if _MSC_VER >= 1400 && !defined(__INTEL_COMPILER)
    00136                 #define MultiplyWords(p, a, b)          p = __emulu(a, b);
    00137         #else
    00138                 #define MultiplyWords(p, a, b)          p = (dword)a*b;
    00139         #endif
    00140         #define AssignWord(a, b)                        a = b;
    00141         #define Add2WordsBy1(a, b, c)           a = b + c;
    00142         #define Acc2WordsBy2(a, b)                      a += b;
    00143         #define LowWord(a)                                      word(a)
    00144         #define HighWord(a)                                     word(a>>WORD_BITS)
    00145         #define Double3Words(c, d)                      d = 2*d + (c>>(WORD_BITS-1)); c *= 2;
    00146         #define AddWithCarry(u, a, b)           u = dword(a) + b + GetCarry(u);
    00147         #define SubtractWithBorrow(u, a, b)     u = dword(a) - b - GetBorrow(u);
    00148         #define GetCarry(u)                                     HighWord(u)
    00149         #define GetBorrow(u)                            word(u>>(WORD_BITS*2-1))
    00150 #endif
    00151 #ifndef MulAcc
    00152         #define MulAcc(c, d, a, b)                      MultiplyWords(p, a, b); Acc2WordsBy1(p, c); c = LowWord(p); Acc2WordsBy1(d, HighWord(p));
    00153 #endif
    00154 #ifndef Acc2WordsBy1
    00155         #define Acc2WordsBy1(a, b)                      Add2WordsBy1(a, a, b)
    00156 #endif
    00157 #ifndef Acc3WordsBy2
    00158         #define Acc3WordsBy2(c, d, e)           Acc2WordsBy1(e, c); c = LowWord(e); Add2WordsBy1(e, d, HighWord(e));
    00159 #endif
    00160
    00161 class DWord
    00162 {
    00163 public:
    00164         DWord() {}
    00165
    00166 #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00167         explicit DWord(word low)
    00168         {
    00169                 m_whole = low;
    00170         }
    00171 #else
    00172         explicit DWord(word low)
    00173         {
    00174                 m_halfs.low = low;
    00175                 m_halfs.high = 0;
    00176         }
    00177 #endif
    00178
    00179         DWord(word low, word high)
    00180         {
    00181                 m_halfs.low = low;
    00182                 m_halfs.high = high;
    00183         }
    00184
    00185         static DWord Multiply(word a, word b)
    00186         {
    00187                 DWord r;
    00188                 #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00189                         r.m_whole = (dword)a * b;
    00190                 #else
    00191                         r.m_halfs.low = _umul128(a, b, &r.m_halfs.high);
    00192                 #endif
    00193                 return r;
    00194         }
    00195
    00196         static DWord MultiplyAndAdd(word a, word b, word c)
    00197         {
    00198                 DWord r = Multiply(a, b);
    00199                 return r += c;
    00200         }
    00201
    00202         DWord & operator+=(word a)
    00203         {
    00204                 #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00205                         m_whole = m_whole + a;
    00206                 #else
    00207                         m_halfs.low += a;
    00208                         m_halfs.high += (m_halfs.low < a);
    00209                 #endif
    00210                 return *this;
    00211         }
    00212
    00213         DWord operator+(word a)
    00214         {
    00215                 DWord r;
    00216                 #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00217                         r.m_whole = m_whole + a;
    00218                 #else
    00219                         r.m_halfs.low = m_halfs.low + a;
    00220                         r.m_halfs.high = m_halfs.high + (r.m_halfs.low < a);
    00221                 #endif
    00222                 return r;
    00223         }
    00224
    00225         DWord operator-(DWord a)
    00226         {
    00227                 DWord r;
    00228                 #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00229                         r.m_whole = m_whole - a.m_whole;
    00230                 #else
    00231                         r.m_halfs.low = m_halfs.low - a.m_halfs.low;
    00232                         r.m_halfs.high = m_halfs.high - a.m_halfs.high - (r.m_halfs.low > m_halfs.low);
    00233                 #endif
    00234                 return r;
    00235         }
    00236
    00237         DWord operator-(word a)
    00238         {
    00239                 DWord r;
    00240                 #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00241                         r.m_whole = m_whole - a;
    00242                 #else
    00243                         r.m_halfs.low = m_halfs.low - a;
    00244                         r.m_halfs.high = m_halfs.high - (r.m_halfs.low > m_halfs.low);
    00245                 #endif
    00246                 return r;
    00247         }
    00248
    00249         // returns quotient, which must fit in a word
    00250         word operator/(word divisor);
    00251
    00252         word operator%(word a);
    00253
    00254         bool operator!() const
    00255         {
    00256         #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00257                 return !m_whole;
    00258         #else
    00259                 return !m_halfs.high && !m_halfs.low;
    00260         #endif
    00261         }
    00262
    00263         word GetLowHalf() const {return m_halfs.low;}
    00264         word GetHighHalf() const {return m_halfs.high;}
    00265         word GetHighHalfAsBorrow() const {return 0-m_halfs.high;}
    00266
    00267 private:
    00268         union
    00269         {
    00270         #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00271                 dword m_whole;
    00272         #endif
    00273                 struct
    00274                 {
    00275                 #ifdef IS_LITTLE_ENDIAN
    00276                         word low;
    00277                         word high;
    00278                 #else
    00279                         word high;
    00280                         word low;
    00281                 #endif
    00282                 } m_halfs;
    00283         };
    00284 };
    00285
    00286 class Word
    00287 {
    00288 public:
    00289         Word() {}
    00290
    00291         Word(word value)
    00292         {
    00293                 m_whole = value;
    00294         }
    00295
    00296         Word(hword low, hword high)
    00297         {
    00298                 m_whole = low | (word(high) << (WORD_BITS/2));
    00299         }
    00300
    00301         static Word Multiply(hword a, hword b)
    00302         {
    00303                 Word r;
    00304                 r.m_whole = (word)a * b;
    00305                 return r;
    00306         }
    00307
    00308         Word operator-(Word a)
    00309         {
    00310                 Word r;
    00311                 r.m_whole = m_whole - a.m_whole;
    00312                 return r;
    00313         }
    00314
    00315         Word operator-(hword a)
    00316         {
    00317                 Word r;
    00318                 r.m_whole = m_whole - a;
    00319                 return r;
    00320         }
    00321
    00322         // returns quotient, which must fit in a word
    00323         hword operator/(hword divisor)
    00324         {
    00325                 return hword(m_whole / divisor);
    00326         }
    00327
    00328         bool operator!() const
    00329         {
    00330                 return !m_whole;
    00331         }
    00332
    00333         word GetWhole() const {return m_whole;}
    00334         hword GetLowHalf() const {return hword(m_whole);}
    00335         hword GetHighHalf() const {return hword(m_whole>>(WORD_BITS/2));}
    00336         hword GetHighHalfAsBorrow() const {return 0-hword(m_whole>>(WORD_BITS/2));}
    00337
    00338 private:
    00339         word m_whole;
    00340 };
    00341
    00342 // do a 3 word by 2 word divide, returns quotient and leaves remainder in A
    00343 template <class S, class D>
    00344 S DivideThreeWordsByTwo(S *A, S B0, S B1, D *dummy=NULL)
    00345 {
    00346         // assert {A[2],A[1]} < {B1,B0}, so quotient can fit in a S
    00347         assert(A[2] < B1 || (A[2]==B1 && A[1] < B0));
    00348
    00349         // estimate the quotient: do a 2 S by 1 S divide
    00350         S Q;
    00351         if (S(B1+1) == 0)
    00352                 Q = A[2];
    00353         else
    00354                 Q = D(A[1], A[2]) / S(B1+1);
    00355
    00356         // now subtract Q*B from A
    00357         D p = D::Multiply(B0, Q);
    00358         D u = (D) A[0] - p.GetLowHalf();
    00359         A[0] = u.GetLowHalf();
    00360         u = (D) A[1] - p.GetHighHalf() - u.GetHighHalfAsBorrow() - D::Multiply(B1, Q);
    00361         A[1] = u.GetLowHalf();
    00362         A[2] += u.GetHighHalf();
    00363
    00364         // Q <= actual quotient, so fix it
    00365         while (A[2] || A[1] > B1 || (A[1]==B1 && A[0]>=B0))
    00366         {
    00367                 u = (D) A[0] - B0;
    00368                 A[0] = u.GetLowHalf();
    00369                 u = (D) A[1] - B1 - u.GetHighHalfAsBorrow();
    00370                 A[1] = u.GetLowHalf();
    00371                 A[2] += u.GetHighHalf();
    00372                 Q++;
    00373                 assert(Q);      // shouldn't overflow
    00374         }
    00375
    00376         return Q;
    00377 }
    00378
    00379 // do a 4 word by 2 word divide, returns 2 word quotient in Q0 and Q1
    00380 template <class S, class D>
    00381 inline D DivideFourWordsByTwo(S *T, const D &Al, const D &Ah, const D &B)
    00382 {
    00383         if (!B) // if divisor is 0, we assume divisor==2**(2*WORD_BITS)
    00384                 return D(Ah.GetLowHalf(), Ah.GetHighHalf());
    00385         else
    00386         {
    00387                 S Q[2];
    00388                 T[0] = Al.GetLowHalf();
    00389                 T[1] = Al.GetHighHalf();
    00390                 T[2] = Ah.GetLowHalf();
    00391                 T[3] = Ah.GetHighHalf();
    00392                 Q[1] = DivideThreeWordsByTwo<S, D>(T+1, B.GetLowHalf(), B.GetHighHalf());
    00393                 Q[0] = DivideThreeWordsByTwo<S, D>(T, B.GetLowHalf(), B.GetHighHalf());
    00394                 return D(Q[0], Q[1]);
    00395         }
    00396 }
    00397
    00398 // returns quotient, which must fit in a word
    00399 inline word DWord::operator/(word a)
    00400 {
    00401         #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00402                 return word(m_whole / a);
    00403         #else
    00404                 hword r[4];
    00405                 return DivideFourWordsByTwo<hword, Word>(r, m_halfs.low, m_halfs.high, a).GetWhole();
    00406         #endif
    00407 }
    00408
    00409 inline word DWord::operator%(word a)
    00410 {
    00411         #ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
    00412                 return word(m_whole % a);
    00413         #else
    00414                 if (a < (word(1) << (WORD_BITS/2)))
    00415                 {
    00416                         hword h = hword(a);
    00417                         word r = m_halfs.high % h;
    00418                         r = ((m_halfs.low >> (WORD_BITS/2)) + (r << (WORD_BITS/2))) % h;
    00419                         return hword((hword(m_halfs.low) + (r << (WORD_BITS/2))) % h);
    00420                 }
    00421                 else
    00422                 {
    00423                         hword r[4];
    00424                         DivideFourWordsByTwo<hword, Word>(r, m_halfs.low, m_halfs.high, a);
    00425                         return Word(r[0], r[1]).GetWhole();
    00426                 }
    00427         #endif
    00428 }
    00429
    00430 // ********************************************************
    00431
    00432 // use some tricks to share assembly code between MSVC and GCC
    00433 #if defined(__GNUC__)
    00434         #define AddPrologue \
    00435                 int result;     \
    00436                 __asm__ __volatile__ \
    00437                 ( \
    00438                         ".intel_syntax noprefix;"
    00439         #define AddEpilogue \
    00440                         ".att_syntax prefix;" \
    00441                                         : "=a" (result)\
    00442                                         : "d" (C), "a" (A), "D" (B), "c" (N) \
    00443                                         : "%esi", "memory", "cc" \
    00444                 );\
    00445                 return result;
    00446         #define MulPrologue \
    00447                 __asm__ __volatile__ \
    00448                 ( \
    00449                         ".intel_syntax noprefix;" \
    00450                         AS1(    push    ebx) \
    00451                         AS2(    mov             ebx, edx)
    00452         #define MulEpilogue \
    00453                         AS1(    pop             ebx) \
    00454                         ".att_syntax prefix;" \
    00455                         : \
    00456                         : "d" (s_maskLow16), "c" (C), "a" (A), "D" (B) \
    00457                         : "%esi", "memory", "cc" \
    00458                 );
    00459         #define SquPrologue             MulPrologue
    00460         #define SquEpilogue     \
    00461                         AS1(    pop             ebx) \
    00462                         ".att_syntax prefix;" \
    00463                         : \
    00464                         : "d" (s_maskLow16), "c" (C), "a" (A) \
    00465                         : "%esi", "%edi", "memory", "cc" \
    00466                 );
    00467         #define TopPrologue             MulPrologue
    00468         #define TopEpilogue     \
    00469                         AS1(    pop             ebx) \
    00470                         ".att_syntax prefix;" \
    00471                         : \
    00472                         : "d" (s_maskLow16), "c" (C), "a" (A), "D" (B), "S" (L) \
    00473                         : "memory", "cc" \
    00474                 );
    00475 #else
    00476         #define AddPrologue \
    00477                 __asm   push edi \
    00478                 __asm   push esi \
    00479                 __asm   mov             eax, [esp+12] \
    00480                 __asm   mov             edi, [esp+16]
    00481         #define AddEpilogue \
    00482                 __asm   pop esi \
    00483                 __asm   pop edi \
    00484                 __asm   ret 8
    00485 #if _MSC_VER < 1300
    00486         #define SaveEBX         __asm push ebx
    00487         #define RestoreEBX      __asm pop ebx
    00488 #else
    00489         #define SaveEBX
    00490         #define RestoreEBX
    00491 #endif
    00492         #define SquPrologue                                     \
    00493                 AS2(    mov             eax, A)                 \
    00494                 AS2(    mov             ecx, C)                 \
    00495                 SaveEBX                                                 \
    00496                 AS2(    lea             ebx, s_maskLow16)
    00497         #define MulPrologue                                     \
    00498                 AS2(    mov             eax, A)                 \
    00499                 AS2(    mov             edi, B)                 \
    00500                 AS2(    mov             ecx, C)                 \
    00501                 SaveEBX                                                 \
    00502                 AS2(    lea             ebx, s_maskLow16)
    00503         #define TopPrologue                                     \
    00504                 AS2(    mov             eax, A)                 \
    00505                 AS2(    mov             edi, B)                 \
    00506                 AS2(    mov             ecx, C)                 \
    00507                 AS2(    mov             esi, L)                 \
    00508                 SaveEBX                                                 \
    00509                 AS2(    lea             ebx, s_maskLow16)
    00510         #define SquEpilogue             RestoreEBX
    00511         #define MulEpilogue             RestoreEBX
    00512         #define TopEpilogue             RestoreEBX
    00513 #endif
    00514
    00515 #ifdef CRYPTOPP_X64_MASM_AVAILABLE
    00516 extern "C" {
    00517 int Baseline_Add(size_t N, word *C, const word *A, const word *B);
    00518 int Baseline_Sub(size_t N, word *C, const word *A, const word *B);
    00519 }
    00520 #elif defined(CRYPTOPP_X64_ASM_AVAILABLE) && defined(__GNUC__)
    00521 int Baseline_Add(size_t N, word *C, const word *A, const word *B)
    00522 {
    00523         word result;
    00524         __asm__ __volatile__
    00525         (
    00526         ".intel_syntax;"
    00527         AS1(    neg             %1)
    00528         ASJ(    jz,             1, f)
    00529         AS2(    mov             %0,[%3+8*%1])
    00530         AS2(    add             %0,[%4+8*%1])
    00531         AS2(    mov             [%2+8*%1],%0)
    00532         ASL(0)
    00533         AS2(    mov             %0,[%3+8*%1+8])
    00534         AS2(    adc             %0,[%4+8*%1+8])
    00535         AS2(    mov             [%2+8*%1+8],%0)
    00536         AS2(    lea             %1,[%1+2])
    00537         ASJ(    jrcxz,  1, f)
    00538         AS2(    mov             %0,[%3+8*%1])
    00539         AS2(    adc             %0,[%4+8*%1])
    00540         AS2(    mov             [%2+8*%1],%0)
    00541         ASJ(    jmp,    0, b)
    00542         ASL(1)
    00543         AS2(    mov             %0, 0)
    00544         AS2(    adc             %0, %0)
    00545         ".att_syntax;"
    00546         : "=&r" (result)
    00547         : "c" (N), "r" (C+N), "r" (A+N), "r" (B+N)
    00548         : "memory", "cc"
    00549         );
    00550         return (int)result;
    00551 }
    00552
    00553 int Baseline_Sub(size_t N, word *C, const word *A, const word *B)
    00554 {
    00555         word result;
    00556         __asm__ __volatile__
    00557         (
    00558         ".intel_syntax;"
    00559         AS1(    neg             %1)
    00560         ASJ(    jz,             1, f)
    00561         AS2(    mov             %0,[%3+8*%1])
    00562         AS2(    sub             %0,[%4+8*%1])
    00563         AS2(    mov             [%2+8*%1],%0)
    00564         ASL(0)
    00565         AS2(    mov             %0,[%3+8*%1+8])
    00566         AS2(    sbb             %0,[%4+8*%1+8])
    00567         AS2(    mov             [%2+8*%1+8],%0)
    00568         AS2(    lea             %1,[%1+2])
    00569         ASJ(    jrcxz,  1, f)
    00570         AS2(    mov             %0,[%3+8*%1])
    00571         AS2(    sbb             %0,[%4+8*%1])
    00572         AS2(    mov             [%2+8*%1],%0)
    00573         ASJ(    jmp,    0, b)
    00574         ASL(1)
    00575         AS2(    mov             %0, 0)
    00576         AS2(    adc             %0, %0)
    00577         ".att_syntax;"
    00578         : "=&r" (result)
    00579         : "c" (N), "r" (C+N), "r" (A+N), "r" (B+N)
    00580         : "memory", "cc"
    00581         );
    00582         return (int)result;
    00583 }
    00584 #elif defined(CRYPTOPP_X86_ASM_AVAILABLE) && CRYPTOPP_BOOL_X86
    00585 CRYPTOPP_NAKED int CRYPTOPP_FASTCALL Baseline_Add(size_t N, word *C, const word *A, const word *B)
    00586 {
    00587         AddPrologue
    00588
    00589         // now: eax = A, edi = B, edx = C, ecx = N
    00590         AS2(    lea             eax, [eax+4*ecx])
    00591         AS2(    lea             edi, [edi+4*ecx])
    00592         AS2(    lea             edx, [edx+4*ecx])
    00593
    00594         AS1(    neg             ecx)                            // ecx is negative index
    00595         AS2(    test    ecx, 2)                         // this clears carry flag
    00596         ASJ(    jz,             0, f)
    00597         AS2(    sub             ecx, 2)
    00598         ASJ(    jmp,    1, f)
    00599
    00600         ASL(0)
    00601         ASJ(    jecxz,  2, f)                           // loop until ecx overflows and becomes zero
    00602         AS2(    mov             esi,[eax+4*ecx])
    00603         AS2(    adc             esi,[edi+4*ecx])
    00604         AS2(    mov             [edx+4*ecx],esi)
    00605         AS2(    mov             esi,[eax+4*ecx+4])
    00606         AS2(    adc             esi,[edi+4*ecx+4])
    00607         AS2(    mov             [edx+4*ecx+4],esi)
    00608         ASL(1)
    00609         AS2(    mov             esi,[eax+4*ecx+8])
    00610         AS2(    adc             esi,[edi+4*ecx+8])
    00611         AS2(    mov             [edx+4*ecx+8],esi)
    00612         AS2(    mov             esi,[eax+4*ecx+12])
    00613         AS2(    adc             esi,[edi+4*ecx+12])
    00614         AS2(    mov             [edx+4*ecx+12],esi)
    00615
    00616         AS2(    lea             ecx,[ecx+4])            // advance index, avoid inc which causes slowdown on Intel Core 2
    00617         ASJ(    jmp,    0, b)
    00618
    00619         ASL(2)
    00620         AS2(    mov             eax, 0)
    00621         AS1(    setc    al)                                     // store carry into eax (return result register)
    00622
    00623         AddEpilogue
    00624 }
    00625
    00626 CRYPTOPP_NAKED int CRYPTOPP_FASTCALL Baseline_Sub(size_t N, word *C, const word *A, const word *B)
    00627 {
    00628         AddPrologue
    00629
    00630         // now: eax = A, edi = B, edx = C, ecx = N
    00631         AS2(    lea             eax, [eax+4*ecx])
    00632         AS2(    lea             edi, [edi+4*ecx])
    00633         AS2(    lea             edx, [edx+4*ecx])
    00634
    00635         AS1(    neg             ecx)                            // ecx is negative index
    00636         AS2(    test    ecx, 2)                         // this clears carry flag
    00637         ASJ(    jz,             0, f)
    00638         AS2(    sub             ecx, 2)
    00639         ASJ(    jmp,    1, f)
    00640
    00641         ASL(0)
    00642         ASJ(    jecxz,  2, f)                           // loop until ecx overflows and becomes zero
    00643         AS2(    mov             esi,[eax+4*ecx])
    00644         AS2(    sbb             esi,[edi+4*ecx])
    00645         AS2(    mov             [edx+4*ecx],esi)
    00646         AS2(    mov             esi,[eax+4*ecx+4])
    00647         AS2(    sbb             esi,[edi+4*ecx+4])
    00648         AS2(    mov             [edx+4*ecx+4],esi)
    00649         ASL(1)
    00650         AS2(    mov             esi,[eax+4*ecx+8])
    00651         AS2(    sbb             esi,[edi+4*ecx+8])
    00652         AS2(    mov             [edx+4*ecx+8],esi)
    00653         AS2(    mov             esi,[eax+4*ecx+12])
    00654         AS2(    sbb             esi,[edi+4*ecx+12])
    00655         AS2(    mov             [edx+4*ecx+12],esi)
    00656
    00657         AS2(    lea             ecx,[ecx+4])            // advance index, avoid inc which causes slowdown on Intel Core 2
    00658         ASJ(    jmp,    0, b)
    00659
    00660         ASL(2)
    00661         AS2(    mov             eax, 0)
    00662         AS1(    setc    al)                                     // store carry into eax (return result register)
    00663
    00664         AddEpilogue
    00665 }
    00666
    00667 #if CRYPTOPP_INTEGER_SSE2
    00668 CRYPTOPP_NAKED int CRYPTOPP_FASTCALL SSE2_Add(size_t N, word *C, const word *A, const word *B)
    00669 {
    00670         AddPrologue
    00671
    00672         // now: eax = A, edi = B, edx = C, ecx = N
    00673         AS2(    lea             eax, [eax+4*ecx])
    00674         AS2(    lea             edi, [edi+4*ecx])
    00675         AS2(    lea             edx, [edx+4*ecx])
    00676
    00677         AS1(    neg             ecx)                            // ecx is negative index
    00678         AS2(    pxor    mm2, mm2)
    00679         ASJ(    jz,             2, f)
    00680         AS2(    test    ecx, 2)                         // this clears carry flag
    00681         ASJ(    jz,             0, f)
    00682         AS2(    sub             ecx, 2)
    00683         ASJ(    jmp,    1, f)
    00684
    00685         ASL(0)
    00686         AS2(    movd     mm0, DWORD PTR [eax+4*ecx])
    00687         AS2(    movd     mm1, DWORD PTR [edi+4*ecx])
    00688         AS2(    paddq    mm0, mm1)
    00689         AS2(    paddq    mm2, mm0)
    00690         AS2(    movd     DWORD PTR [edx+4*ecx], mm2)
    00691         AS2(    psrlq    mm2, 32)
    00692
    00693         AS2(    movd     mm0, DWORD PTR [eax+4*ecx+4])
    00694         AS2(    movd     mm1, DWORD PTR [edi+4*ecx+4])
    00695         AS2(    paddq    mm0, mm1)
    00696         AS2(    paddq    mm2, mm0)
    00697         AS2(    movd     DWORD PTR [edx+4*ecx+4], mm2)
    00698         AS2(    psrlq    mm2, 32)
    00699
    00700         ASL(1)
    00701         AS2(    movd     mm0, DWORD PTR [eax+4*ecx+8])
    00702         AS2(    movd     mm1, DWORD PTR [edi+4*ecx+8])
    00703         AS2(    paddq    mm0, mm1)
    00704         AS2(    paddq    mm2, mm0)
    00705         AS2(    movd     DWORD PTR [edx+4*ecx+8], mm2)
    00706         AS2(    psrlq    mm2, 32)
    00707
    00708         AS2(    movd     mm0, DWORD PTR [eax+4*ecx+12])
    00709         AS2(    movd     mm1, DWORD PTR [edi+4*ecx+12])
    00710         AS2(    paddq    mm0, mm1)
    00711         AS2(    paddq    mm2, mm0)
    00712         AS2(    movd     DWORD PTR [edx+4*ecx+12], mm2)
    00713         AS2(    psrlq    mm2, 32)
    00714
    00715         AS2(    add             ecx, 4)
    00716         ASJ(    jnz,    0, b)
    00717
    00718         ASL(2)
    00719         AS2(    movd    eax, mm2)
    00720         AS1(    emms)
    00721
    00722         AddEpilogue
    00723 }
    00724 CRYPTOPP_NAKED int CRYPTOPP_FASTCALL SSE2_Sub(size_t N, word *C, const word *A, const word *B)
    00725 {
    00726         AddPrologue
    00727
    00728         // now: eax = A, edi = B, edx = C, ecx = N
    00729         AS2(    lea             eax, [eax+4*ecx])
    00730         AS2(    lea             edi, [edi+4*ecx])
    00731         AS2(    lea             edx, [edx+4*ecx])
    00732
    00733         AS1(    neg             ecx)                            // ecx is negative index
    00734         AS2(    pxor    mm2, mm2)
    00735         ASJ(    jz,             2, f)
    00736         AS2(    test    ecx, 2)                         // this clears carry flag
    00737         ASJ(    jz,             0, f)
    00738         AS2(    sub             ecx, 2)
    00739         ASJ(    jmp,    1, f)
    00740
    00741         ASL(0)
    00742         AS2(    movd     mm0, DWORD PTR [eax+4*ecx])
    00743         AS2(    movd     mm1, DWORD PTR [edi+4*ecx])
    00744         AS2(    psubq    mm0, mm1)
    00745         AS2(    psubq    mm0, mm2)
    00746         AS2(    movd     DWORD PTR [edx+4*ecx], mm0)
    00747         AS2(    psrlq    mm0, 63)
    00748
    00749         AS2(    movd     mm2, DWORD PTR [eax+4*ecx+4])
    00750         AS2(    movd     mm1, DWORD PTR [edi+4*ecx+4])
    00751         AS2(    psubq    mm2, mm1)
    00752         AS2(    psubq    mm2, mm0)
    00753         AS2(    movd     DWORD PTR [edx+4*ecx+4], mm2)
    00754         AS2(    psrlq    mm2, 63)
    00755
    00756         ASL(1)
    00757         AS2(    movd     mm0, DWORD PTR [eax+4*ecx+8])
    00758         AS2(    movd     mm1, DWORD PTR [edi+4*ecx+8])
    00759         AS2(    psubq    mm0, mm1)
    00760         AS2(    psubq    mm0, mm2)
    00761         AS2(    movd     DWORD PTR [edx+4*ecx+8], mm0)
    00762         AS2(    psrlq    mm0, 63)
    00763
    00764         AS2(    movd     mm2, DWORD PTR [eax+4*ecx+12])
    00765         AS2(    movd     mm1, DWORD PTR [edi+4*ecx+12])
    00766         AS2(    psubq    mm2, mm1)
    00767         AS2(    psubq    mm2, mm0)
    00768         AS2(    movd     DWORD PTR [edx+4*ecx+12], mm2)
    00769         AS2(    psrlq    mm2, 63)
    00770
    00771         AS2(    add             ecx, 4)
    00772         ASJ(    jnz,    0, b)
    00773
    00774         ASL(2)
    00775         AS2(    movd    eax, mm2)
    00776         AS1(    emms)
    00777
    00778         AddEpilogue
    00779 }
    00780 #endif  // #if CRYPTOPP_BOOL_SSE2_ASM_AVAILABLE
    00781 #else
    00782 int CRYPTOPP_FASTCALL Baseline_Add(size_t N, word *C, const word *A, const word *B)
    00783 {
    00784         assert (N%2 == 0);
    00785
    00786         Declare2Words(u);
    00787         AssignWord(u, 0);
    00788         for (size_t i=0; i<N; i+=2)
    00789         {
    00790                 AddWithCarry(u, A[i], B[i]);
    00791                 C[i] = LowWord(u);
    00792                 AddWithCarry(u, A[i+1], B[i+1]);
    00793                 C[i+1] = LowWord(u);
    00794         }
    00795         return int(GetCarry(u));
    00796 }
    00797
    00798 int CRYPTOPP_FASTCALL Baseline_Sub(size_t N, word *C, const word *A, const word *B)
    00799 {
    00800         assert (N%2 == 0);
    00801
    00802         Declare2Words(u);
    00803         AssignWord(u, 0);
    00804         for (size_t i=0; i<N; i+=2)
    00805         {
    00806                 SubtractWithBorrow(u, A[i], B[i]);
    00807                 C[i] = LowWord(u);
    00808                 SubtractWithBorrow(u, A[i+1], B[i+1]);
    00809                 C[i+1] = LowWord(u);
    00810         }
    00811         return int(GetBorrow(u));
    00812 }
    00813 #endif
    00814
    00815 static word LinearMultiply(word *C, const word *A, word B, size_t N)
    00816 {
    00817         word carry=0;
    00818         for(unsigned i=0; i<N; i++)
    00819         {
    00820                 Declare2Words(p);
    00821                 MultiplyWords(p, A[i], B);
    00822                 Acc2WordsBy1(p, carry);
    00823                 C[i] = LowWord(p);
    00824                 carry = HighWord(p);
    00825         }
    00826         return carry;
    00827 }
    00828
    00829 #ifndef CRYPTOPP_DOXYGEN_PROCESSING
    00830
    00831 #define Mul_2 \
    00832         Mul_Begin(2) \
    00833         Mul_SaveAcc(0, 0, 1) Mul_Acc(1, 0) \
    00834         Mul_End(1, 1)
    00835
    00836 #define Mul_4 \
    00837         Mul_Begin(4) \
    00838         Mul_SaveAcc(0, 0, 1) Mul_Acc(1, 0) \
    00839         Mul_SaveAcc(1, 0, 2) Mul_Acc(1, 1) Mul_Acc(2, 0)  \
    00840         Mul_SaveAcc(2, 0, 3) Mul_Acc(1, 2) Mul_Acc(2, 1) Mul_Acc(3, 0)  \
    00841         Mul_SaveAcc(3, 1, 3) Mul_Acc(2, 2) Mul_Acc(3, 1)  \
    00842         Mul_SaveAcc(4, 2, 3) Mul_Acc(3, 2) \
    00843         Mul_End(5, 3)
    00844
    00845 #define Mul_8 \
    00846         Mul_Begin(8) \
    00847         Mul_SaveAcc(0, 0, 1) Mul_Acc(1, 0) \
    00848         Mul_SaveAcc(1, 0, 2) Mul_Acc(1, 1) Mul_Acc(2, 0)  \
    00849         Mul_SaveAcc(2, 0, 3) Mul_Acc(1, 2) Mul_Acc(2, 1) Mul_Acc(3, 0)  \
    00850         Mul_SaveAcc(3, 0, 4) Mul_Acc(1, 3) Mul_Acc(2, 2) Mul_Acc(3, 1) Mul_Acc(4, 0) \
    00851         Mul_SaveAcc(4, 0, 5) Mul_Acc(1, 4) Mul_Acc(2, 3) Mul_Acc(3, 2) Mul_Acc(4, 1) Mul_Acc(5, 0) \
    00852         Mul_SaveAcc(5, 0, 6) Mul_Acc(1, 5) Mul_Acc(2, 4) Mul_Acc(3, 3) Mul_Acc(4, 2) Mul_Acc(5, 1) Mul_Acc(6, 0) \
    00853         Mul_SaveAcc(6, 0, 7) Mul_Acc(1, 6) Mul_Acc(2, 5) Mul_Acc(3, 4) Mul_Acc(4, 3) Mul_Acc(5, 2) Mul_Acc(6, 1) Mul_Acc(7, 0) \
    00854         Mul_SaveAcc(7, 1, 7) Mul_Acc(2, 6) Mul_Acc(3, 5) Mul_Acc(4, 4) Mul_Acc(5, 3) Mul_Acc(6, 2) Mul_Acc(7, 1) \
    00855         Mul_SaveAcc(8, 2, 7) Mul_Acc(3, 6) Mul_Acc(4, 5) Mul_Acc(5, 4) Mul_Acc(6, 3) Mul_Acc(7, 2) \
    00856         Mul_SaveAcc(9, 3, 7) Mul_Acc(4, 6) Mul_Acc(5, 5) Mul_Acc(6, 4) Mul_Acc(7, 3) \
    00857         Mul_SaveAcc(10, 4, 7) Mul_Acc(5, 6) Mul_Acc(6, 5) Mul_Acc(7, 4) \
    00858         Mul_SaveAcc(11, 5, 7) Mul_Acc(6, 6) Mul_Acc(7, 5) \
    00859         Mul_SaveAcc(12, 6, 7) Mul_Acc(7, 6) \
    00860         Mul_End(13, 7)
    00861
    00862 #define Mul_16 \
    00863         Mul_Begin(16) \
    00864         Mul_SaveAcc(0, 0, 1) Mul_Acc(1, 0) \
    00865         Mul_SaveAcc(1, 0, 2) Mul_Acc(1, 1) Mul_Acc(2, 0) \
    00866         Mul_SaveAcc(2, 0, 3) Mul_Acc(1, 2) Mul_Acc(2, 1) Mul_Acc(3, 0) \
    00867         Mul_SaveAcc(3, 0, 4) Mul_Acc(1, 3) Mul_Acc(2, 2) Mul_Acc(3, 1) Mul_Acc(4, 0) \
    00868         Mul_SaveAcc(4, 0, 5) Mul_Acc(1, 4) Mul_Acc(2, 3) Mul_Acc(3, 2) Mul_Acc(4, 1) Mul_Acc(5, 0) \
    00869         Mul_SaveAcc(5, 0, 6) Mul_Acc(1, 5) Mul_Acc(2, 4) Mul_Acc(3, 3) Mul_Acc(4, 2) Mul_Acc(5, 1) Mul_Acc(6, 0) \
    00870         Mul_SaveAcc(6, 0, 7) Mul_Acc(1, 6) Mul_Acc(2, 5) Mul_Acc(3, 4) Mul_Acc(4, 3) Mul_Acc(5, 2) Mul_Acc(6, 1) Mul_Acc(7, 0) \
    00871         Mul_SaveAcc(7, 0, 8) Mul_Acc(1, 7) Mul_Acc(2, 6) Mul_Acc(3, 5) Mul_Acc(4, 4) Mul_Acc(5, 3) Mul_Acc(6, 2) Mul_Acc(7, 1) Mul_Acc(8, 0) \
    00872         Mul_SaveAcc(8, 0, 9) Mul_Acc(1, 8) Mul_Acc(2, 7) Mul_Acc(3, 6) Mul_Acc(4, 5) Mul_Acc(5, 4) Mul_Acc(6, 3) Mul_Acc(7, 2) Mul_Acc(8, 1) Mul_Acc(9, 0) \
    00873         Mul_SaveAcc(9, 0, 10) Mul_Acc(1, 9) Mul_Acc(2, 8) Mul_Acc(3, 7) Mul_Acc(4, 6) Mul_Acc(5, 5) Mul_Acc(6, 4) Mul_Acc(7, 3) Mul_Acc(8, 2) Mul_Acc(9, 1) Mul_Acc(10, 0) \
    00874         Mul_SaveAcc(10, 0, 11) Mul_Acc(1, 10) Mul_Acc(2, 9) Mul_Acc(3, 8) Mul_Acc(4, 7) Mul_Acc(5, 6) Mul_Acc(6, 5) Mul_Acc(7, 4) Mul_Acc(8, 3) Mul_Acc(9, 2) Mul_Acc(10, 1) Mul_Acc(11, 0) \
    00875         Mul_SaveAcc(11, 0, 12) Mul_Acc(1, 11) Mul_Acc(2, 10) Mul_Acc(3, 9) Mul_Acc(4, 8) Mul_Acc(5, 7) Mul_Acc(6, 6) Mul_Acc(7, 5) Mul_Acc(8, 4) Mul_Acc(9, 3) Mul_Acc(10, 2) Mul_Acc(11, 1) Mul_Acc(12, 0) \
    00876         Mul_SaveAcc(12, 0, 13) Mul_Acc(1, 12) Mul_Acc(2, 11) Mul_Acc(3, 10) Mul_Acc(4, 9) Mul_Acc(5, 8) Mul_Acc(6, 7) Mul_Acc(7, 6) Mul_Acc(8, 5) Mul_Acc(9, 4) Mul_Acc(10, 3) Mul_Acc(11, 2) Mul_Acc(12, 1) Mul_Acc(13, 0) \
    00877         Mul_SaveAcc(13, 0, 14) Mul_Acc(1, 13) Mul_Acc(2, 12) Mul_Acc(3, 11) Mul_Acc(4, 10) Mul_Acc(5, 9) Mul_Acc(6, 8) Mul_Acc(7, 7) Mul_Acc(8, 6) Mul_Acc(9, 5) Mul_Acc(10, 4) Mul_Acc(11, 3) Mul_Acc(12, 2) Mul_Acc(13, 1) Mul_Acc(14, 0) \
    00878         Mul_SaveAcc(14, 0, 15) Mul_Acc(1, 14) Mul_Acc(2, 13) Mul_Acc(3, 12) Mul_Acc(4, 11) Mul_Acc(5, 10) Mul_Acc(6, 9) Mul_Acc(7, 8) Mul_Acc(8, 7) Mul_Acc(9, 6) Mul_Acc(10, 5) Mul_Acc(11, 4) Mul_Acc(12, 3) Mul_Acc(13, 2) Mul_Acc(14, 1) Mul_Acc(15, 0) \
    00879         Mul_SaveAcc(15, 1, 15) Mul_Acc(2, 14) Mul_Acc(3, 13) Mul_Acc(4, 12) Mul_Acc(5, 11) Mul_Acc(6, 10) Mul_Acc(7, 9) Mul_Acc(8, 8) Mul_Acc(9, 7) Mul_Acc(10, 6) Mul_Acc(11, 5) Mul_Acc(12, 4) Mul_Acc(13, 3) Mul_Acc(14, 2) Mul_Acc(15, 1) \
    00880         Mul_SaveAcc(16, 2, 15) Mul_Acc(3, 14) Mul_Acc(4, 13) Mul_Acc(5, 12) Mul_Acc(6, 11) Mul_Acc(7, 10) Mul_Acc(8, 9) Mul_Acc(9, 8) Mul_Acc(10, 7) Mul_Acc(11, 6) Mul_Acc(12, 5) Mul_Acc(13, 4) Mul_Acc(14, 3) Mul_Acc(15, 2) \
    00881         Mul_SaveAcc(17, 3, 15) Mul_Acc(4, 14) Mul_Acc(5, 13) Mul_Acc(6, 12) Mul_Acc(7, 11) Mul_Acc(8, 10) Mul_Acc(9, 9) Mul_Acc(10, 8) Mul_Acc(11, 7) Mul_Acc(12, 6) Mul_Acc(13, 5) Mul_Acc(14, 4) Mul_Acc(15, 3) \
    00882         Mul_SaveAcc(18, 4, 15) Mul_Acc(5, 14) Mul_Acc(6, 13) Mul_Acc(7, 12) Mul_Acc(8, 11) Mul_Acc(9, 10) Mul_Acc(10, 9) Mul_Acc(11, 8) Mul_Acc(12, 7) Mul_Acc(13, 6) Mul_Acc(14, 5) Mul_Acc(15, 4) \
    00883         Mul_SaveAcc(19, 5, 15) Mul_Acc(6, 14) Mul_Acc(7, 13) Mul_Acc(8, 12) Mul_Acc(9, 11) Mul_Acc(10, 10) Mul_Acc(11, 9) Mul_Acc(12, 8) Mul_Acc(13, 7) Mul_Acc(14, 6) Mul_Acc(15, 5) \
    00884         Mul_SaveAcc(20, 6, 15) Mul_Acc(7, 14) Mul_Acc(8, 13) Mul_Acc(9, 12) Mul_Acc(10, 11) Mul_Acc(11, 10) Mul_Acc(12, 9) Mul_Acc(13, 8) Mul_Acc(14, 7) Mul_Acc(15, 6) \
    00885         Mul_SaveAcc(21, 7, 15) Mul_Acc(8, 14) Mul_Acc(9, 13) Mul_Acc(10, 12) Mul_Acc(11, 11) Mul_Acc(12, 10) Mul_Acc(13, 9) Mul_Acc(14, 8) Mul_Acc(15, 7) \
    00886         Mul_SaveAcc(22, 8, 15) Mul_Acc(9, 14) Mul_Acc(10, 13) Mul_Acc(11, 12) Mul_Acc(12, 11) Mul_Acc(13, 10) Mul_Acc(14, 9) Mul_Acc(15, 8) \
    00887         Mul_SaveAcc(23, 9, 15) Mul_Acc(10, 14) Mul_Acc(11, 13) Mul_Acc(12, 12) Mul_Acc(13, 11) Mul_Acc(14, 10) Mul_Acc(15, 9) \
    00888         Mul_SaveAcc(24, 10, 15) Mul_Acc(11, 14) Mul_Acc(12, 13) Mul_Acc(13, 12) Mul_Acc(14, 11) Mul_Acc(15, 10) \
    00889         Mul_SaveAcc(25, 11, 15) Mul_Acc(12, 14) Mul_Acc(13, 13) Mul_Acc(14, 12) Mul_Acc(15, 11) \
    00890         Mul_SaveAcc(26, 12, 15) Mul_Acc(13, 14) Mul_Acc(14, 13) Mul_Acc(15, 12) \
    00891         Mul_SaveAcc(27, 13, 15) Mul_Acc(14, 14) Mul_Acc(15, 13) \
    00892         Mul_SaveAcc(28, 14, 15) Mul_Acc(15, 14) \
    00893         Mul_End(29, 15)
    00894
    00895 #define Squ_2 \
    00896         Squ_Begin(2) \
    00897         Squ_End(2)
    00898
    00899 #define Squ_4 \
    00900         Squ_Begin(4) \
    00901         Squ_SaveAcc(1, 0, 2) Squ_Diag(1) \
    00902         Squ_SaveAcc(2, 0, 3) Squ_Acc(1, 2) Squ_NonDiag \
    00903         Squ_SaveAcc(3, 1, 3) Squ_Diag(2) \
    00904         Squ_SaveAcc(4, 2, 3) Squ_NonDiag \
    00905         Squ_End(4)
    00906
    00907 #define Squ_8 \
    00908         Squ_Begin(8) \
    00909         Squ_SaveAcc(1, 0, 2) Squ_Diag(1) \
    00910         Squ_SaveAcc(2, 0, 3) Squ_Acc(1, 2) Squ_NonDiag \
    00911         Squ_SaveAcc(3, 0, 4) Squ_Acc(1, 3) Squ_Diag(2) \
    00912         Squ_SaveAcc(4, 0, 5) Squ_Acc(1, 4) Squ_Acc(2, 3) Squ_NonDiag \
    00913         Squ_SaveAcc(5, 0, 6) Squ_Acc(1, 5) Squ_Acc(2, 4) Squ_Diag(3) \
    00914         Squ_SaveAcc(6, 0, 7) Squ_Acc(1, 6) Squ_Acc(2, 5) Squ_Acc(3, 4) Squ_NonDiag \
    00915         Squ_SaveAcc(7, 1, 7) Squ_Acc(2, 6) Squ_Acc(3, 5) Squ_Diag(4) \
    00916         Squ_SaveAcc(8, 2, 7) Squ_Acc(3, 6) Squ_Acc(4, 5)  Squ_NonDiag \
    00917         Squ_SaveAcc(9, 3, 7) Squ_Acc(4, 6) Squ_Diag(5) \
    00918         Squ_SaveAcc(10, 4, 7) Squ_Acc(5, 6) Squ_NonDiag \
    00919         Squ_SaveAcc(11, 5, 7) Squ_Diag(6) \
    00920         Squ_SaveAcc(12, 6, 7) Squ_NonDiag \
    00921         Squ_End(8)
    00922
    00923 #define Squ_16 \
    00924         Squ_Begin(16) \
    00925         Squ_SaveAcc(1, 0, 2) Squ_Diag(1) \
    00926         Squ_SaveAcc(2, 0, 3) Squ_Acc(1, 2) Squ_NonDiag \
    00927         Squ_SaveAcc(3, 0, 4) Squ_Acc(1, 3) Squ_Diag(2) \
    00928         Squ_SaveAcc(4, 0, 5) Squ_Acc(1, 4) Squ_Acc(2, 3) Squ_NonDiag \
    00929         Squ_SaveAcc(5, 0, 6) Squ_Acc(1, 5) Squ_Acc(2, 4) Squ_Diag(3) \
    00930         Squ_SaveAcc(6, 0, 7) Squ_Acc(1, 6) Squ_Acc(2, 5) Squ_Acc(3, 4) Squ_NonDiag \
    00931         Squ_SaveAcc(7, 0, 8) Squ_Acc(1, 7) Squ_Acc(2, 6) Squ_Acc(3, 5) Squ_Diag(4) \
    00932         Squ_SaveAcc(8, 0, 9) Squ_Acc(1, 8) Squ_Acc(2, 7) Squ_Acc(3, 6) Squ_Acc(4, 5) Squ_NonDiag \
    00933         Squ_SaveAcc(9, 0, 10) Squ_Acc(1, 9) Squ_Acc(2, 8) Squ_Acc(3, 7) Squ_Acc(4, 6) Squ_Diag(5) \
    00934         Squ_SaveAcc(10, 0, 11) Squ_Acc(1, 10) Squ_Acc(2, 9) Squ_Acc(3, 8) Squ_Acc(4, 7) Squ_Acc(5, 6) Squ_NonDiag \
    00935         Squ_SaveAcc(11, 0, 12) Squ_Acc(1, 11) Squ_Acc(2, 10) Squ_Acc(3, 9) Squ_Acc(4, 8) Squ_Acc(5, 7) Squ_Diag(6) \
    00936         Squ_SaveAcc(12, 0, 13) Squ_Acc(1, 12) Squ_Acc(2, 11) Squ_Acc(3, 10) Squ_Acc(4, 9) Squ_Acc(5, 8) Squ_Acc(6, 7) Squ_NonDiag \
    00937         Squ_SaveAcc(13, 0, 14) Squ_Acc(1, 13) Squ_Acc(2, 12) Squ_Acc(3, 11) Squ_Acc(4, 10) Squ_Acc(5, 9) Squ_Acc(6, 8) Squ_Diag(7) \
    00938         Squ_SaveAcc(14, 0, 15) Squ_Acc(1, 14) Squ_Acc(2, 13) Squ_Acc(3, 12) Squ_Acc(4, 11) Squ_Acc(5, 10) Squ_Acc(6, 9) Squ_Acc(7, 8) Squ_NonDiag \
    00939         Squ_SaveAcc(15, 1, 15) Squ_Acc(2, 14) Squ_Acc(3, 13) Squ_Acc(4, 12) Squ_Acc(5, 11) Squ_Acc(6, 10) Squ_Acc(7, 9) Squ_Diag(8) \
    00940         Squ_SaveAcc(16, 2, 15) Squ_Acc(3, 14) Squ_Acc(4, 13) Squ_Acc(5, 12) Squ_Acc(6, 11) Squ_Acc(7, 10) Squ_Acc(8, 9) Squ_NonDiag \
    00941         Squ_SaveAcc(17, 3, 15) Squ_Acc(4, 14) Squ_Acc(5, 13) Squ_Acc(6, 12) Squ_Acc(7, 11) Squ_Acc(8, 10) Squ_Diag(9) \
    00942         Squ_SaveAcc(18, 4, 15) Squ_Acc(5, 14) Squ_Acc(6, 13) Squ_Acc(7, 12) Squ_Acc(8, 11) Squ_Acc(9, 10) Squ_NonDiag \
    00943         Squ_SaveAcc(19, 5, 15) Squ_Acc(6, 14) Squ_Acc(7, 13) Squ_Acc(8, 12) Squ_Acc(9, 11) Squ_Diag(10) \
    00944         Squ_SaveAcc(20, 6, 15) Squ_Acc(7, 14) Squ_Acc(8, 13) Squ_Acc(9, 12) Squ_Acc(10, 11) Squ_NonDiag \
    00945         Squ_SaveAcc(21, 7, 15) Squ_Acc(8, 14) Squ_Acc(9, 13) Squ_Acc(10, 12) Squ_Diag(11) \
    00946         Squ_SaveAcc(22, 8, 15) Squ_Acc(9, 14) Squ_Acc(10, 13) Squ_Acc(11, 12) Squ_NonDiag \
    00947         Squ_SaveAcc(23, 9, 15) Squ_Acc(10, 14) Squ_Acc(11, 13) Squ_Diag(12) \
    00948         Squ_SaveAcc(24, 10, 15) Squ_Acc(11, 14) Squ_Acc(12, 13) Squ_NonDiag \
    00949         Squ_SaveAcc(25, 11, 15) Squ_Acc(12, 14) Squ_Diag(13) \
    00950         Squ_SaveAcc(26, 12, 15) Squ_Acc(13, 14) Squ_NonDiag \
    00951         Squ_SaveAcc(27, 13, 15) Squ_Diag(14) \
    00952         Squ_SaveAcc(28, 14, 15) Squ_NonDiag \
    00953         Squ_End(16)
    00954
    00955 #define Bot_2 \
    00956         Mul_Begin(2) \
    00957         Bot_SaveAcc(0, 0, 1) Bot_Acc(1, 0) \
    00958         Bot_End(2)
    00959
    00960 #define Bot_4 \
    00961         Mul_Begin(4) \
    00962         Mul_SaveAcc(0, 0, 1) Mul_Acc(1, 0) \
    00963         Mul_SaveAcc(1, 2, 0) Mul_Acc(1, 1) Mul_Acc(0, 2)  \
    00964         Bot_SaveAcc(2, 0, 3) Bot_Acc(1, 2) Bot_Acc(2, 1) Bot_Acc(3, 0)  \
    00965         Bot_End(4)
    00966
    00967 #define Bot_8 \
    00968         Mul_Begin(8) \
    00969         Mul_SaveAcc(0, 0, 1) Mul_Acc(1, 0) \
    00970         Mul_SaveAcc(1, 0, 2) Mul_Acc(1, 1) Mul_Acc(2, 0)  \
    00971         Mul_SaveAcc(2, 0, 3) Mul_Acc(1, 2) Mul_Acc(2, 1) Mul_Acc(3, 0)  \
    00972         Mul_SaveAcc(3, 0, 4) Mul_Acc(1, 3) Mul_Acc(2, 2) Mul_Acc(3, 1) Mul_Acc(4, 0) \
    00973         Mul_SaveAcc(4, 0, 5) Mul_Acc(1, 4) Mul_Acc(2, 3) Mul_Acc(3, 2) Mul_Acc(4, 1) Mul_Acc(5, 0) \
    00974         Mul_SaveAcc(5, 0, 6) Mul_Acc(1, 5) Mul_Acc(2, 4) Mul_Acc(3, 3) Mul_Acc(4, 2) Mul_Acc(5, 1) Mul_Acc(6, 0) \
    00975         Bot_SaveAcc(6, 0, 7) Bot_Acc(1, 6) Bot_Acc(2, 5) Bot_Acc(3, 4) Bot_Acc(4, 3) Bot_Acc(5, 2) Bot_Acc(6, 1) Bot_Acc(7, 0) \
    00976         Bot_End(8)
    00977
    00978 #define Bot_16 \
    00979         Mul_Begin(16) \
    00980         Mul_SaveAcc(0, 0, 1) Mul_Acc(1, 0) \
    00981         Mul_SaveAcc(1, 0, 2) Mul_Acc(1, 1) Mul_Acc(2, 0) \
    00982         Mul_SaveAcc(2, 0, 3) Mul_Acc(1, 2) Mul_Acc(2, 1) Mul_Acc(3, 0) \
    00983         Mul_SaveAcc(3, 0, 4) Mul_Acc(1, 3) Mul_Acc(2, 2) Mul_Acc(3, 1) Mul_Acc(4, 0) \
    00984         Mul_SaveAcc(4, 0, 5) Mul_Acc(1, 4) Mul_Acc(2, 3) Mul_Acc(3, 2) Mul_Acc(4, 1) Mul_Acc(5, 0) \
    00985         Mul_SaveAcc(5, 0, 6) Mul_Acc(1, 5) Mul_Acc(2, 4) Mul_Acc(3, 3) Mul_Acc(4, 2) Mul_Acc(5, 1) Mul_Acc(6, 0) \
    00986         Mul_SaveAcc(6, 0, 7) Mul_Acc(1, 6) Mul_Acc(2, 5) Mul_Acc(3, 4) Mul_Acc(4, 3) Mul_Acc(5, 2) Mul_Acc(6, 1) Mul_Acc(7, 0) \
    00987         Mul_SaveAcc(7, 0, 8) Mul_Acc(1, 7) Mul_Acc(2, 6) Mul_Acc(3, 5) Mul_Acc(4, 4) Mul_Acc(5, 3) Mul_Acc(6, 2) Mul_Acc(7, 1) Mul_Acc(8, 0) \
    00988         Mul_SaveAcc(8, 0, 9) Mul_Acc(1, 8) Mul_Acc(2, 7) Mul_Acc(3, 6) Mul_Acc(4, 5) Mul_Acc(5, 4) Mul_Acc(6, 3) Mul_Acc(7, 2) Mul_Acc(8, 1) Mul_Acc(9, 0) \
    00989         Mul_SaveAcc(9, 0, 10) Mul_Acc(1, 9) Mul_Acc(2, 8) Mul_Acc(3, 7) Mul_Acc(4, 6) Mul_Acc(5, 5) Mul_Acc(6, 4) Mul_Acc(7, 3) Mul_Acc(8, 2) Mul_Acc(9, 1) Mul_Acc(10, 0) \
    00990         Mul_SaveAcc(10, 0, 11) Mul_Acc(1, 10) Mul_Acc(2, 9) Mul_Acc(3, 8) Mul_Acc(4, 7) Mul_Acc(5, 6) Mul_Acc(6, 5) Mul_Acc(7, 4) Mul_Acc(8, 3) Mul_Acc(9, 2) Mul_Acc(10, 1) Mul_Acc(11, 0) \
    00991         Mul_SaveAcc(11, 0, 12) Mul_Acc(1, 11) Mul_Acc(2, 10) Mul_Acc(3, 9) Mul_Acc(4, 8) Mul_Acc(5, 7) Mul_Acc(6, 6) Mul_Acc(7, 5) Mul_Acc(8, 4) Mul_Acc(9, 3) Mul_Acc(10, 2) Mul_Acc(11, 1) Mul_Acc(12, 0) \
    00992         Mul_SaveAcc(12, 0, 13) Mul_Acc(1, 12) Mul_Acc(2, 11) Mul_Acc(3, 10) Mul_Acc(4, 9) Mul_Acc(5, 8) Mul_Acc(6, 7) Mul_Acc(7, 6) Mul_Acc(8, 5) Mul_Acc(9, 4) Mul_Acc(10, 3) Mul_Acc(11, 2) Mul_Acc(12, 1) Mul_Acc(13, 0) \
    00993         Mul_SaveAcc(13, 0, 14) Mul_Acc(1, 13) Mul_Acc(2, 12) Mul_Acc(3, 11) Mul_Acc(4, 10) Mul_Acc(5, 9) Mul_Acc(6, 8) Mul_Acc(7, 7) Mul_Acc(8, 6) Mul_Acc(9, 5) Mul_Acc(10, 4) Mul_Acc(11, 3) Mul_Acc(12, 2) Mul_Acc(13, 1) Mul_Acc(14, 0) \
    00994         Bot_SaveAcc(14, 0, 15) Bot_Acc(1, 14) Bot_Acc(2, 13) Bot_Acc(3, 12) Bot_Acc(4, 11) Bot_Acc(5, 10) Bot_Acc(6, 9) Bot_Acc(7, 8) Bot_Acc(8, 7) Bot_Acc(9, 6) Bot_Acc(10, 5) Bot_Acc(11, 4) Bot_Acc(12, 3) Bot_Acc(13, 2) Bot_Acc(14, 1) Bot_Acc(15, 0) \
    00995         Bot_End(16)
    00996
    00997 #endif
    00998
    00999 #if 0
    01000 #define Mul_Begin(n)                            \
    01001         Declare2Words(p)                                \
    01002         Declare2Words(c)                                \
    01003         Declare2Words(d)                                \
    01004         MultiplyWords(p, A[0], B[0])    \
    01005         AssignWord(c, LowWord(p))               \
    01006         AssignWord(d, HighWord(p))
    01007
    01008 #define Mul_Acc(i, j)                           \
    01009         MultiplyWords(p, A[i], B[j])    \
    01010         Acc2WordsBy1(c, LowWord(p))             \
    01011         Acc2WordsBy1(d, HighWord(p))
    01012
    01013 #define Mul_SaveAcc(k, i, j)            \
    01014         R[k] = LowWord(c);                              \
    01015         Add2WordsBy1(c, d, HighWord(c)) \
    01016         MultiplyWords(p, A[i], B[j])    \
    01017         AssignWord(d, HighWord(p))              \
    01018         Acc2WordsBy1(c, LowWord(p))
    01019
    01020 #define Mul_End(n)                                      \
    01021         R[2*n-3] = LowWord(c);                  \
    01022         Acc2WordsBy1(d, HighWord(c))    \
    01023         MultiplyWords(p, A[n-1], B[n-1])\
    01024         Acc2WordsBy2(d, p)                              \
    01025         R[2*n-2] = LowWord(d);                  \
    01026         R[2*n-1] = HighWord(d);
    01027
    01028 #define Bot_SaveAcc(k, i, j)            \
    01029         R[k] = LowWord(c);                              \
    01030         word e = LowWord(d) + HighWord(c);      \
    01031         e += A[i] * B[j];
    01032
    01033 #define Bot_Acc(i, j)   \
    01034         e += A[i] * B[j];
    01035
    01036 #define Bot_End(n)              \
    01037         R[n-1] = e;
    01038 #else
    01039 #define Mul_Begin(n)                            \
    01040         Declare2Words(p)                                \
    01041         word c; \
    01042         Declare2Words(d)                                \
    01043         MultiplyWords(p, A[0], B[0])    \
    01044         c = LowWord(p);         \
    01045         AssignWord(d, HighWord(p))
    01046
    01047 #define Mul_Acc(i, j)                           \
    01048         MulAcc(c, d, A[i], B[j])
    01049
    01050 #define Mul_SaveAcc(k, i, j)            \
    01051         R[k] = c;                               \
    01052         c = LowWord(d); \
    01053         AssignWord(d, HighWord(d))      \
    01054         MulAcc(c, d, A[i], B[j])
    01055
    01056 #define Mul_End(k, i)                                   \
    01057         R[k] = c;                       \
    01058         MultiplyWords(p, A[i], B[i])    \
    01059         Acc2WordsBy2(p, d)                              \
    01060         R[k+1] = LowWord(p);                    \
    01061         R[k+2] = HighWord(p);
    01062
    01063 #define Bot_SaveAcc(k, i, j)            \
    01064         R[k] = c;                               \
    01065         c = LowWord(d); \
    01066         c += A[i] * B[j];
    01067
    01068 #define Bot_Acc(i, j)   \
    01069         c += A[i] * B[j];
    01070
    01071 #define Bot_End(n)              \
    01072         R[n-1] = c;
    01073 #endif
    01074
    01075 #define Squ_Begin(n)                            \
    01076         Declare2Words(p)                                \
    01077         word c;                         \
    01078         Declare2Words(d)                                \
    01079         Declare2Words(e)                                \
    01080         MultiplyWords(p, A[0], A[0])    \
    01081         R[0] = LowWord(p);                              \
    01082         AssignWord(e, HighWord(p))              \
    01083         MultiplyWords(p, A[0], A[1])    \
    01084         c = LowWord(p);         \
    01085         AssignWord(d, HighWord(p))              \
    01086         Squ_NonDiag                                             \
    01087
    01088 #define Squ_NonDiag                             \
    01089         Double3Words(c, d)
    01090
    01091 #define Squ_SaveAcc(k, i, j)            \
    01092         Acc3WordsBy2(c, d, e)                   \
    01093         R[k] = c;                               \
    01094         MultiplyWords(p, A[i], A[j])    \
    01095         c = LowWord(p);         \
    01096         AssignWord(d, HighWord(p))              \
    01097
    01098 #define Squ_Acc(i, j)                           \
    01099         MulAcc(c, d, A[i], A[j])
    01100
    01101 #define Squ_Diag(i)                                     \
    01102         Squ_NonDiag                                             \
    01103         MulAcc(c, d, A[i], A[i])
    01104
    01105 #define Squ_End(n)                                      \
    01106         Acc3WordsBy2(c, d, e)                   \
    01107         R[2*n-3] = c;                   \
    01108         MultiplyWords(p, A[n-1], A[n-1])\
    01109         Acc2WordsBy2(p, e)                              \
    01110         R[2*n-2] = LowWord(p);                  \
    01111         R[2*n-1] = HighWord(p);
    01112
    01113 void Baseline_Multiply2(word *R, const word *A, const word *B)
    01114 {
    01115         Mul_2
    01116 }
    01117
    01118 void Baseline_Multiply4(word *R, const word *A, const word *B)
    01119 {
    01120         Mul_4
    01121 }
    01122
    01123 void Baseline_Multiply8(word *R, const word *A, const word *B)
    01124 {
    01125         Mul_8
    01126 }
    01127
    01128 void Baseline_Square2(word *R, const word *A)
    01129 {
    01130         Squ_2
    01131 }
    01132
    01133 void Baseline_Square4(word *R, const word *A)
    01134 {
    01135         Squ_4
    01136 }
    01137
    01138 void Baseline_Square8(word *R, const word *A)
    01139 {
    01140         Squ_8
    01141 }
    01142
    01143 void Baseline_MultiplyBottom2(word *R, const word *A, const word *B)
    01144 {
    01145         Bot_2
    01146 }
    01147
    01148 void Baseline_MultiplyBottom4(word *R, const word *A, const word *B)
    01149 {
    01150         Bot_4
    01151 }
    01152
    01153 void Baseline_MultiplyBottom8(word *R, const word *A, const word *B)
    01154 {
    01155         Bot_8
    01156 }
    01157
    01158 #define Top_Begin(n)                            \
    01159         Declare2Words(p)                                \
    01160         word c; \
    01161         Declare2Words(d)                                \
    01162         MultiplyWords(p, A[0], B[n-2]);\
    01163         AssignWord(d, HighWord(p));
    01164
    01165 #define Top_Acc(i, j)   \
    01166         MultiplyWords(p, A[i], B[j]);\
    01167         Acc2WordsBy1(d, HighWord(p));
    01168
    01169 #define Top_SaveAcc0(i, j)              \
    01170         c = LowWord(d); \
    01171         AssignWord(d, HighWord(d))      \
    01172         MulAcc(c, d, A[i], B[j])
    01173
    01174 #define Top_SaveAcc1(i, j)              \
    01175         c = L<c; \
    01176         Acc2WordsBy1(d, c);     \
    01177         c = LowWord(d); \
    01178         AssignWord(d, HighWord(d))      \
    01179         MulAcc(c, d, A[i], B[j])
    01180
    01181 void Baseline_MultiplyTop2(word *R, const word *A, const word *B, word L)
    01182 {
    01183         word T[4];
    01184         Baseline_Multiply2(T, A, B);
    01185         R[0] = T[2];
    01186         R[1] = T[3];
    01187 }
    01188
    01189 void Baseline_MultiplyTop4(word *R, const word *A, const word *B, word L)
    01190 {
    01191         Top_Begin(4)
    01192         Top_Acc(1, 1) Top_Acc(2, 0)  \
    01193         Top_SaveAcc0(0, 3) Mul_Acc(1, 2) Mul_Acc(2, 1) Mul_Acc(3, 0)  \
    01194         Top_SaveAcc1(1, 3) Mul_Acc(2, 2) Mul_Acc(3, 1)  \
    01195         Mul_SaveAcc(0, 2, 3) Mul_Acc(3, 2) \
    01196         Mul_End(1, 3)
    01197 }
    01198
    01199 void Baseline_MultiplyTop8(word *R, const word *A, const word *B, word L)
    01200 {
    01201         Top_Begin(8)
    01202         Top_Acc(1, 5) Top_Acc(2, 4) Top_Acc(3, 3) Top_Acc(4, 2) Top_Acc(5, 1) Top_Acc(6, 0) \
    01203         Top_SaveAcc0(0, 7) Mul_Acc(1, 6) Mul_Acc(2, 5) Mul_Acc(3, 4) Mul_Acc(4, 3) Mul_Acc(5, 2) Mul_Acc(6, 1) Mul_Acc(7, 0) \
    01204         Top_SaveAcc1(1, 7) Mul_Acc(2, 6) Mul_Acc(3, 5) Mul_Acc(4, 4) Mul_Acc(5, 3) Mul_Acc(6, 2) Mul_Acc(7, 1) \
    01205         Mul_SaveAcc(0, 2, 7) Mul_Acc(3, 6) Mul_Acc(4, 5) Mul_Acc(5, 4) Mul_Acc(6, 3) Mul_Acc(7, 2) \
    01206         Mul_SaveAcc(1, 3, 7) Mul_Acc(4, 6) Mul_Acc(5, 5) Mul_Acc(6, 4) Mul_Acc(7, 3) \
    01207         Mul_SaveAcc(2, 4, 7) Mul_Acc(5, 6) Mul_Acc(6, 5) Mul_Acc(7, 4) \
    01208         Mul_SaveAcc(3, 5, 7) Mul_Acc(6, 6) Mul_Acc(7, 5) \
    01209         Mul_SaveAcc(4, 6, 7) Mul_Acc(7, 6) \
    01210         Mul_End(5, 7)
    01211 }
    01212
    01213 #if !CRYPTOPP_INTEGER_SSE2      // save memory by not compiling these functions when SSE2 is available
    01214 void Baseline_Multiply16(word *R, const word *A, const word *B)
    01215 {
    01216         Mul_16
    01217 }
    01218
    01219 void Baseline_Square16(word *R, const word *A)
    01220 {
    01221         Squ_16
    01222 }
    01223
    01224 void Baseline_MultiplyBottom16(word *R, const word *A, const word *B)
    01225 {
    01226         Bot_16
    01227 }
    01228
    01229 void Baseline_MultiplyTop16(word *R, const word *A, const word *B, word L)
    01230 {
    01231         Top_Begin(16)
    01232         Top_Acc(1, 13) Top_Acc(2, 12) Top_Acc(3, 11) Top_Acc(4, 10) Top_Acc(5, 9) Top_Acc(6, 8) Top_Acc(7, 7) Top_Acc(8, 6) Top_Acc(9, 5) Top_Acc(10, 4) Top_Acc(11, 3) Top_Acc(12, 2) Top_Acc(13, 1) Top_Acc(14, 0) \
    01233         Top_SaveAcc0(0, 15) Mul_Acc(1, 14) Mul_Acc(2, 13) Mul_Acc(3, 12) Mul_Acc(4, 11) Mul_Acc(5, 10) Mul_Acc(6, 9) Mul_Acc(7, 8) Mul_Acc(8, 7) Mul_Acc(9, 6) Mul_Acc(10, 5) Mul_Acc(11, 4) Mul_Acc(12, 3) Mul_Acc(13, 2) Mul_Acc(14, 1) Mul_Acc(15, 0) \
    01234         Top_SaveAcc1(1, 15) Mul_Acc(2, 14) Mul_Acc(3, 13) Mul_Acc(4, 12) Mul_Acc(5, 11) Mul_Acc(6, 10) Mul_Acc(7, 9) Mul_Acc(8, 8) Mul_Acc(9, 7) Mul_Acc(10, 6) Mul_Acc(11, 5) Mul_Acc(12, 4) Mul_Acc(13, 3) Mul_Acc(14, 2) Mul_Acc(15, 1) \
    01235         Mul_SaveAcc(0, 2, 15) Mul_Acc(3, 14) Mul_Acc(4, 13) Mul_Acc(5, 12) Mul_Acc(6, 11) Mul_Acc(7, 10) Mul_Acc(8, 9) Mul_Acc(9, 8) Mul_Acc(10, 7) Mul_Acc(11, 6) Mul_Acc(12, 5) Mul_Acc(13, 4) Mul_Acc(14, 3) Mul_Acc(15, 2) \
    01236         Mul_SaveAcc(1, 3, 15) Mul_Acc(4, 14) Mul_Acc(5, 13) Mul_Acc(6, 12) Mul_Acc(7, 11) Mul_Acc(8, 10) Mul_Acc(9, 9) Mul_Acc(10, 8) Mul_Acc(11, 7) Mul_Acc(12, 6) Mul_Acc(13, 5) Mul_Acc(14, 4) Mul_Acc(15, 3) \
    01237         Mul_SaveAcc(2, 4, 15) Mul_Acc(5, 14) Mul_Acc(6, 13) Mul_Acc(7, 12) Mul_Acc(8, 11) Mul_Acc(9, 10) Mul_Acc(10, 9) Mul_Acc(11, 8) Mul_Acc(12, 7) Mul_Acc(13, 6) Mul_Acc(14, 5) Mul_Acc(15, 4) \
    01238         Mul_SaveAcc(3, 5, 15) Mul_Acc(6, 14) Mul_Acc(7, 13) Mul_Acc(8, 12) Mul_Acc(9, 11) Mul_Acc(10, 10) Mul_Acc(11, 9) Mul_Acc(12, 8) Mul_Acc(13, 7) Mul_Acc(14, 6) Mul_Acc(15, 5) \
    01239         Mul_SaveAcc(4, 6, 15) Mul_Acc(7, 14) Mul_Acc(8, 13) Mul_Acc(9, 12) Mul_Acc(10, 11) Mul_Acc(11, 10) Mul_Acc(12, 9) Mul_Acc(13, 8) Mul_Acc(14, 7) Mul_Acc(15, 6) \
    01240         Mul_SaveAcc(5, 7, 15) Mul_Acc(8, 14) Mul_Acc(9, 13) Mul_Acc(10, 12) Mul_Acc(11, 11) Mul_Acc(12, 10) Mul_Acc(13, 9) Mul_Acc(14, 8) Mul_Acc(15, 7) \
    01241         Mul_SaveAcc(6, 8, 15) Mul_Acc(9, 14) Mul_Acc(10, 13) Mul_Acc(11, 12) Mul_Acc(12, 11) Mul_Acc(13, 10) Mul_Acc(14, 9) Mul_Acc(15, 8) \
    01242         Mul_SaveAcc(7, 9, 15) Mul_Acc(10, 14) Mul_Acc(11, 13) Mul_Acc(12, 12) Mul_Acc(13, 11) Mul_Acc(14, 10) Mul_Acc(15, 9) \
    01243         Mul_SaveAcc(8, 10, 15) Mul_Acc(11, 14) Mul_Acc(12, 13) Mul_Acc(13, 12) Mul_Acc(14, 11) Mul_Acc(15, 10) \
    01244         Mul_SaveAcc(9, 11, 15) Mul_Acc(12, 14) Mul_Acc(13, 13) Mul_Acc(14, 12) Mul_Acc(15, 11) \
    01245         Mul_SaveAcc(10, 12, 15) Mul_Acc(13, 14) Mul_Acc(14, 13) Mul_Acc(15, 12) \
    01246         Mul_SaveAcc(11, 13, 15) Mul_Acc(14, 14) Mul_Acc(15, 13) \
    01247         Mul_SaveAcc(12, 14, 15) Mul_Acc(15, 14) \
    01248         Mul_End(13, 15)
    01249 }
    01250 #endif
    01251
    01252 // ********************************************************
    01253
    01254 #if CRYPTOPP_INTEGER_SSE2
    01255
    01256 CRYPTOPP_ALIGN_DATA(16) static const word32 s_maskLow16[4] CRYPTOPP_SECTION_ALIGN16 = {0xffff,0xffff,0xffff,0xffff};
    01257
    01258 #undef Mul_Begin
    01259 #undef Mul_Acc
    01260 #undef Top_Begin
    01261 #undef Top_Acc
    01262 #undef Squ_Acc
    01263 #undef Squ_NonDiag
    01264 #undef Squ_Diag
    01265 #undef Squ_SaveAcc
    01266 #undef Squ_Begin
    01267 #undef Mul_SaveAcc
    01268 #undef Bot_Acc
    01269 #undef Bot_SaveAcc
    01270 #undef Bot_End
    01271 #undef Squ_End
    01272 #undef Mul_End
    01273
    01274 #define SSE2_FinalSave(k)                       \
    01275         AS2(    psllq           xmm5, 16)       \
    01276         AS2(    paddq           xmm4, xmm5)     \
    01277         AS2(    movq            QWORD PTR [ecx+8*(k)], xmm4)
    01278
    01279 #define SSE2_SaveShift(k)                       \
    01280         AS2(    movq            xmm0, xmm6)     \
    01281         AS2(    punpckhqdq      xmm6, xmm0)     \
    01282         AS2(    movq            xmm1, xmm7)     \
    01283         AS2(    punpckhqdq      xmm7, xmm1)     \
    01284         AS2(    paddd           xmm6, xmm0)     \
    01285         AS2(    pslldq          xmm6, 4)        \
    01286         AS2(    paddd           xmm7, xmm1)     \
    01287         AS2(    paddd           xmm4, xmm6)     \
    01288         AS2(    pslldq          xmm7, 4)        \
    01289         AS2(    movq            xmm6, xmm4)     \
    01290         AS2(    paddd           xmm5, xmm7)     \
    01291         AS2(    movq            xmm7, xmm5)     \
    01292         AS2(    movd            DWORD PTR [ecx+8*(k)], xmm4)    \
    01293         AS2(    psrlq           xmm6, 16)       \
    01294         AS2(    paddq           xmm6, xmm7)     \
    01295         AS2(    punpckhqdq      xmm4, xmm0)     \
    01296         AS2(    punpckhqdq      xmm5, xmm0)     \
    01297         AS2(    movq            QWORD PTR [ecx+8*(k)+2], xmm6)  \
    01298         AS2(    psrlq           xmm6, 3*16)     \
    01299         AS2(    paddd           xmm4, xmm6)     \
    01300
    01301 #define Squ_SSE2_SaveShift(k)                   \
    01302         AS2(    movq            xmm0, xmm6)     \
    01303         AS2(    punpckhqdq      xmm6, xmm0)     \
    01304         AS2(    movq            xmm1, xmm7)     \
    01305         AS2(    punpckhqdq      xmm7, xmm1)     \
    01306         AS2(    paddd           xmm6, xmm0)     \
    01307         AS2(    pslldq          xmm6, 4)        \
    01308         AS2(    paddd           xmm7, xmm1)     \
    01309         AS2(    paddd           xmm4, xmm6)     \
    01310         AS2(    pslldq          xmm7, 4)        \
    01311         AS2(    movhlps         xmm6, xmm4)     \
    01312         AS2(    movd            DWORD PTR [ecx+8*(k)], xmm4)    \
    01313         AS2(    paddd           xmm5, xmm7)     \
    01314         AS2(    movhps          QWORD PTR [esp+12], xmm5)\
    01315         AS2(    psrlq           xmm4, 16)       \
    01316         AS2(    paddq           xmm4, xmm5)     \
    01317         AS2(    movq            QWORD PTR [ecx+8*(k)+2], xmm4)  \
    01318         AS2(    psrlq           xmm4, 3*16)     \
    01319         AS2(    paddd           xmm4, xmm6)     \
    01320         AS2(    movq            QWORD PTR [esp+4], xmm4)\
    01321
    01322 #define SSE2_FirstMultiply(i)                           \
    01323         AS2(    movdqa          xmm7, [esi+(i)*16])\
    01324         AS2(    movdqa          xmm5, [edi-(i)*16])\
    01325         AS2(    pmuludq         xmm5, xmm7)             \
    01326         AS2(    movdqa          xmm4, [ebx])\
    01327         AS2(    movdqa          xmm6, xmm4)             \
    01328         AS2(    pand            xmm4, xmm5)             \
    01329         AS2(    psrld           xmm5, 16)               \
    01330         AS2(    pmuludq         xmm7, [edx-(i)*16])\
    01331         AS2(    pand            xmm6, xmm7)             \
    01332         AS2(    psrld           xmm7, 16)
    01333
    01334 #define Squ_Begin(n)                                                    \
    01335         SquPrologue                                                                     \
    01336         AS2(    mov             esi, esp)\
    01337         AS2(    and             esp, 0xfffffff0)\
    01338         AS2(    lea             edi, [esp-32*n])\
    01339         AS2(    sub             esp, 32*n+16)\
    01340         AS1(    push    esi)\
    01341         AS2(    mov             esi, edi)                                       \
    01342         AS2(    xor             edx, edx)                                       \
    01343         ASL(1)                                                                          \
    01344         ASS(    pshufd  xmm0, [eax+edx], 3,1,2,0)       \
    01345         ASS(    pshufd  xmm1, [eax+edx], 2,0,3,1)       \
    01346         AS2(    movdqa  [edi+2*edx], xmm0)              \
    01347         AS2(    psrlq   xmm0, 32)                                       \
    01348         AS2(    movdqa  [edi+2*edx+16], xmm0)   \
    01349         AS2(    movdqa  [edi+16*n+2*edx], xmm1)         \
    01350         AS2(    psrlq   xmm1, 32)                                       \
    01351         AS2(    movdqa  [edi+16*n+2*edx+16], xmm1)      \
    01352         AS2(    add             edx, 16)                                        \
    01353         AS2(    cmp             edx, 8*(n))                                     \
    01354         ASJ(    jne,    1, b)                                           \
    01355         AS2(    lea             edx, [edi+16*n])\
    01356         SSE2_FirstMultiply(0)                                                   \
    01357
    01358 #define Squ_Acc(i)                                                              \
    01359         ASL(LSqu##i)                                                            \
    01360         AS2(    movdqa          xmm1, [esi+(i)*16])     \
    01361         AS2(    movdqa          xmm0, [edi-(i)*16])     \
    01362         AS2(    movdqa          xmm2, [ebx])    \
    01363         AS2(    pmuludq         xmm0, xmm1)                             \
    01364         AS2(    pmuludq         xmm1, [edx-(i)*16])     \
    01365         AS2(    movdqa          xmm3, xmm2)                     \
    01366         AS2(    pand            xmm2, xmm0)                     \
    01367         AS2(    psrld           xmm0, 16)                       \
    01368         AS2(    paddd           xmm4, xmm2)                     \
    01369         AS2(    paddd           xmm5, xmm0)                     \
    01370         AS2(    pand            xmm3, xmm1)                     \
    01371         AS2(    psrld           xmm1, 16)                       \
    01372         AS2(    paddd           xmm6, xmm3)                     \
    01373         AS2(    paddd           xmm7, xmm1)             \
    01374
    01375 #define Squ_Acc1(i)
    01376 #define Squ_Acc2(i)             ASC(call, LSqu##i)
    01377 #define Squ_Acc3(i)             Squ_Acc2(i)
    01378 #define Squ_Acc4(i)             Squ_Acc2(i)
    01379 #define Squ_Acc5(i)             Squ_Acc2(i)
    01380 #define Squ_Acc6(i)             Squ_Acc2(i)
    01381 #define Squ_Acc7(i)             Squ_Acc2(i)
    01382 #define Squ_Acc8(i)             Squ_Acc2(i)
    01383
    01384 #define SSE2_End(E, n)                                  \
    01385         SSE2_SaveShift(2*(n)-3)                 \
    01386         AS2(    movdqa          xmm7, [esi+16]) \
    01387         AS2(    movdqa          xmm0, [edi])    \
    01388         AS2(    pmuludq         xmm0, xmm7)                             \
    01389         AS2(    movdqa          xmm2, [ebx])            \
    01390         AS2(    pmuludq         xmm7, [edx])    \
    01391         AS2(    movdqa          xmm6, xmm2)                             \
    01392         AS2(    pand            xmm2, xmm0)                             \
    01393         AS2(    psrld           xmm0, 16)                               \
    01394         AS2(    paddd           xmm4, xmm2)                             \
    01395         AS2(    paddd           xmm5, xmm0)                             \
    01396         AS2(    pand            xmm6, xmm7)                             \
    01397         AS2(    psrld           xmm7, 16)       \
    01398         SSE2_SaveShift(2*(n)-2)                 \
    01399         SSE2_FinalSave(2*(n)-1)                 \
    01400         AS1(    pop             esp)\
    01401         E
    01402
    01403 #define Squ_End(n)              SSE2_End(SquEpilogue, n)
    01404 #define Mul_End(n)              SSE2_End(MulEpilogue, n)
    01405 #define Top_End(n)              SSE2_End(TopEpilogue, n)
    01406
    01407 #define Squ_Column1(k, i)       \
    01408         Squ_SSE2_SaveShift(k)                                   \
    01409         AS2(    add                     esi, 16)        \
    01410         SSE2_FirstMultiply(1)\
    01411         Squ_Acc##i(i)   \
    01412         AS2(    paddd           xmm4, xmm4)             \
    01413         AS2(    paddd           xmm5, xmm5)             \
    01414         AS2(    movdqa          xmm3, [esi])                            \
    01415         AS2(    movq            xmm1, QWORD PTR [esi+8])        \
    01416         AS2(    pmuludq         xmm1, xmm3)             \
    01417         AS2(    pmuludq         xmm3, xmm3)             \
    01418         AS2(    movdqa          xmm0, [ebx])\
    01419         AS2(    movdqa          xmm2, xmm0)             \
    01420         AS2(    pand            xmm0, xmm1)             \
    01421         AS2(    psrld           xmm1, 16)               \
    01422         AS2(    paddd           xmm6, xmm0)             \
    01423         AS2(    paddd           xmm7, xmm1)             \
    01424         AS2(    pand            xmm2, xmm3)             \
    01425         AS2(    psrld           xmm3, 16)               \
    01426         AS2(    paddd           xmm6, xmm6)             \
    01427         AS2(    paddd           xmm7, xmm7)             \
    01428         AS2(    paddd           xmm4, xmm2)             \
    01429         AS2(    paddd           xmm5, xmm3)             \
    01430         AS2(    movq            xmm0, QWORD PTR [esp+4])\
    01431         AS2(    movq            xmm1, QWORD PTR [esp+12])\
    01432         AS2(    paddd           xmm4, xmm0)\
    01433         AS2(    paddd           xmm5, xmm1)\
    01434
    01435 #define Squ_Column0(k, i)       \
    01436         Squ_SSE2_SaveShift(k)                                   \
    01437         AS2(    add                     edi, 16)        \
    01438         AS2(    add                     edx, 16)        \
    01439         SSE2_FirstMultiply(1)\
    01440         Squ_Acc##i(i)   \
    01441         AS2(    paddd           xmm6, xmm6)             \
    01442         AS2(    paddd           xmm7, xmm7)             \
    01443         AS2(    paddd           xmm4, xmm4)             \
    01444         AS2(    paddd           xmm5, xmm5)             \
    01445         AS2(    movq            xmm0, QWORD PTR [esp+4])\
    01446         AS2(    movq            xmm1, QWORD PTR [esp+12])\
    01447         AS2(    paddd           xmm4, xmm0)\
    01448         AS2(    paddd           xmm5, xmm1)\
    01449
    01450 #define SSE2_MulAdd45                                           \
    01451         AS2(    movdqa          xmm7, [esi])    \
    01452         AS2(    movdqa          xmm0, [edi])    \
    01453         AS2(    pmuludq         xmm0, xmm7)                             \
    01454         AS2(    movdqa          xmm2, [ebx])            \
    01455         AS2(    pmuludq         xmm7, [edx])    \
    01456         AS2(    movdqa          xmm6, xmm2)                             \
    01457         AS2(    pand            xmm2, xmm0)                             \
    01458         AS2(    psrld           xmm0, 16)                               \
    01459         AS2(    paddd           xmm4, xmm2)                             \
    01460         AS2(    paddd           xmm5, xmm0)                             \
    01461         AS2(    pand            xmm6, xmm7)                             \
    01462         AS2(    psrld           xmm7, 16)
    01463
    01464 #define Mul_Begin(n)                                                    \
    01465         MulPrologue                                                                     \
    01466         AS2(    mov             esi, esp)\
    01467         AS2(    and             esp, 0xfffffff0)\
    01468         AS2(    sub             esp, 48*n+16)\
    01469         AS1(    push    esi)\
    01470         AS2(    xor             edx, edx)                                       \
    01471         ASL(1)                                                                          \
    01472         ASS(    pshufd  xmm0, [eax+edx], 3,1,2,0)       \
    01473         ASS(    pshufd  xmm1, [eax+edx], 2,0,3,1)       \
    01474         ASS(    pshufd  xmm2, [edi+edx], 3,1,2,0)       \
    01475         AS2(    movdqa  [esp+20+2*edx], xmm0)           \
    01476         AS2(    psrlq   xmm0, 32)                                       \
    01477         AS2(    movdqa  [esp+20+2*edx+16], xmm0)        \
    01478         AS2(    movdqa  [esp+20+16*n+2*edx], xmm1)              \
    01479         AS2(    psrlq   xmm1, 32)                                       \
    01480         AS2(    movdqa  [esp+20+16*n+2*edx+16], xmm1)   \
    01481         AS2(    movdqa  [esp+20+32*n+2*edx], xmm2)              \
    01482         AS2(    psrlq   xmm2, 32)                                       \
    01483         AS2(    movdqa  [esp+20+32*n+2*edx+16], xmm2)   \
    01484         AS2(    add             edx, 16)                                        \
    01485         AS2(    cmp             edx, 8*(n))                                     \
    01486         ASJ(    jne,    1, b)                                           \
    01487         AS2(    lea             edi, [esp+20])\
    01488         AS2(    lea             edx, [esp+20+16*n])\
    01489         AS2(    lea             esi, [esp+20+32*n])\
    01490         SSE2_FirstMultiply(0)                                                   \
    01491
    01492 #define Mul_Acc(i)                                                              \
    01493         ASL(LMul##i)                                                                            \
    01494         AS2(    movdqa          xmm1, [esi+i/2*(1-(i-2*(i/2))*2)*16])   \
    01495         AS2(    movdqa          xmm0, [edi-i/2*(1-(i-2*(i/2))*2)*16])   \
    01496         AS2(    movdqa          xmm2, [ebx])    \
    01497         AS2(    pmuludq         xmm0, xmm1)                             \
    01498         AS2(    pmuludq         xmm1, [edx-i/2*(1-(i-2*(i/2))*2)*16])   \
    01499         AS2(    movdqa          xmm3, xmm2)                     \
    01500         AS2(    pand            xmm2, xmm0)                     \
    01501         AS2(    psrld           xmm0, 16)                       \
    01502         AS2(    paddd           xmm4, xmm2)                     \
    01503         AS2(    paddd           xmm5, xmm0)                     \
    01504         AS2(    pand            xmm3, xmm1)                     \
    01505         AS2(    psrld           xmm1, 16)                       \
    01506         AS2(    paddd           xmm6, xmm3)                     \
    01507         AS2(    paddd           xmm7, xmm1)             \
    01508
    01509 #define Mul_Acc1(i)
    01510 #define Mul_Acc2(i)             ASC(call, LMul##i)
    01511 #define Mul_Acc3(i)             Mul_Acc2(i)
    01512 #define Mul_Acc4(i)             Mul_Acc2(i)
    01513 #define Mul_Acc5(i)             Mul_Acc2(i)
    01514 #define Mul_Acc6(i)             Mul_Acc2(i)
    01515 #define Mul_Acc7(i)             Mul_Acc2(i)
    01516 #define Mul_Acc8(i)             Mul_Acc2(i)
    01517 #define Mul_Acc9(i)             Mul_Acc2(i)
    01518 #define Mul_Acc10(i)    Mul_Acc2(i)
    01519 #define Mul_Acc11(i)    Mul_Acc2(i)
    01520 #define Mul_Acc12(i)    Mul_Acc2(i)
    01521 #define Mul_Acc13(i)    Mul_Acc2(i)
    01522 #define Mul_Acc14(i)    Mul_Acc2(i)
    01523 #define Mul_Acc15(i)    Mul_Acc2(i)
    01524 #define Mul_Acc16(i)    Mul_Acc2(i)
    01525
    01526 #define Mul_Column1(k, i)       \
    01527         SSE2_SaveShift(k)                                       \
    01528         AS2(    add                     esi, 16)        \
    01529         SSE2_MulAdd45\
    01530         Mul_Acc##i(i)   \
    01531
    01532 #define Mul_Column0(k, i)       \
    01533         SSE2_SaveShift(k)                                       \
    01534         AS2(    add                     edi, 16)        \
    01535         AS2(    add                     edx, 16)        \
    01536         SSE2_MulAdd45\
    01537         Mul_Acc##i(i)   \
    01538
    01539 #define Bot_Acc(i)                                                      \
    01540         AS2(    movdqa          xmm1, [esi+i/2*(1-(i-2*(i/2))*2)*16])   \
    01541         AS2(    movdqa          xmm0, [edi-i/2*(1-(i-2*(i/2))*2)*16])   \
    01542         AS2(    pmuludq         xmm0, xmm1)                             \
    01543         AS2(    pmuludq         xmm1, [edx-i/2*(1-(i-2*(i/2))*2)*16])           \
    01544         AS2(    paddq           xmm4, xmm0)                             \
    01545         AS2(    paddd           xmm6, xmm1)
    01546
    01547 #define Bot_SaveAcc(k)                                  \
    01548         SSE2_SaveShift(k)                                                       \
    01549         AS2(    add                     edi, 16)        \
    01550         AS2(    add                     edx, 16)        \
    01551         AS2(    movdqa          xmm6, [esi])    \
    01552         AS2(    movdqa          xmm0, [edi])    \
    01553         AS2(    pmuludq         xmm0, xmm6)                             \
    01554         AS2(    paddq           xmm4, xmm0)                             \
    01555         AS2(    psllq           xmm5, 16)                               \
    01556         AS2(    paddq           xmm4, xmm5)                             \
    01557         AS2(    pmuludq         xmm6, [edx])
    01558
    01559 #define Bot_End(n)                                                      \
    01560         AS2(    movhlps         xmm7, xmm6)                     \
    01561         AS2(    paddd           xmm6, xmm7)                     \
    01562         AS2(    psllq           xmm6, 32)                       \
    01563         AS2(    paddd           xmm4, xmm6)                     \
    01564         AS2(    movq            QWORD PTR [ecx+8*((n)-1)], xmm4)        \
    01565         AS1(    pop             esp)\
    01566         MulEpilogue
    01567
    01568 #define Top_Begin(n)                                                    \
    01569         TopPrologue                                                                     \
    01570         AS2(    mov             edx, esp)\
    01571         AS2(    and             esp, 0xfffffff0)\
    01572         AS2(    sub             esp, 48*n+16)\
    01573         AS1(    push    edx)\
    01574         AS2(    xor             edx, edx)                                       \
    01575         ASL(1)                                                                          \
    01576         ASS(    pshufd  xmm0, [eax+edx], 3,1,2,0)       \
    01577         ASS(    pshufd  xmm1, [eax+edx], 2,0,3,1)       \
    01578         ASS(    pshufd  xmm2, [edi+edx], 3,1,2,0)       \
    01579         AS2(    movdqa  [esp+20+2*edx], xmm0)           \
    01580         AS2(    psrlq   xmm0, 32)                                       \
    01581         AS2(    movdqa  [esp+20+2*edx+16], xmm0)        \
    01582         AS2(    movdqa  [esp+20+16*n+2*edx], xmm1)              \
    01583         AS2(    psrlq   xmm1, 32)                                       \
    01584         AS2(    movdqa  [esp+20+16*n+2*edx+16], xmm1)   \
    01585         AS2(    movdqa  [esp+20+32*n+2*edx], xmm2)              \
    01586         AS2(    psrlq   xmm2, 32)                                       \
    01587         AS2(    movdqa  [esp+20+32*n+2*edx+16], xmm2)   \
    01588         AS2(    add             edx, 16)                                        \
    01589         AS2(    cmp             edx, 8*(n))                                     \
    01590         ASJ(    jne,    1, b)                                           \
    01591         AS2(    mov             eax, esi)                                       \
    01592         AS2(    lea             edi, [esp+20+00*n+16*(n/2-1)])\
    01593         AS2(    lea             edx, [esp+20+16*n+16*(n/2-1)])\
    01594         AS2(    lea             esi, [esp+20+32*n+16*(n/2-1)])\
    01595         AS2(    pxor    xmm4, xmm4)\
    01596         AS2(    pxor    xmm5, xmm5)
    01597
    01598 #define Top_Acc(i)                                                      \
    01599         AS2(    movq            xmm0, QWORD PTR [esi+i/2*(1-(i-2*(i/2))*2)*16+8])       \
    01600         AS2(    pmuludq         xmm0, [edx-i/2*(1-(i-2*(i/2))*2)*16])   \
    01601         AS2(    psrlq           xmm0, 48)                               \
    01602         AS2(    paddd           xmm5, xmm0)\
    01603
    01604 #define Top_Column0(i)  \
    01605         AS2(    psllq           xmm5, 32)                               \
    01606         AS2(    add                     edi, 16)        \
    01607         AS2(    add                     edx, 16)        \
    01608         SSE2_MulAdd45\
    01609         Mul_Acc##i(i)   \
    01610
    01611 #define Top_Column1(i)  \
    01612         SSE2_SaveShift(0)                                       \
    01613         AS2(    add                     esi, 16)        \
    01614         SSE2_MulAdd45\
    01615         Mul_Acc##i(i)   \
    01616         AS2(    shr                     eax, 16)        \
    01617         AS2(    movd            xmm0, eax)\
    01618         AS2(    movd            xmm1, [ecx+4])\
    01619         AS2(    psrld           xmm1, 16)\
    01620         AS2(    pcmpgtd         xmm1, xmm0)\
    01621         AS2(    psrld           xmm1, 31)\
    01622         AS2(    paddd           xmm4, xmm1)\
    01623
    01624 void SSE2_Square4(word *C, const word *A)
    01625 {
    01626         Squ_Begin(2)
    01627         Squ_Column0(0, 1)
    01628         Squ_End(2)
    01629 }
    01630
    01631 void SSE2_Square8(word *C, const word *A)
    01632 {
    01633         Squ_Begin(4)
    01634 #ifndef __GNUC__
    01635         ASJ(    jmp,    0, f)
    01636         Squ_Acc(2)
    01637         AS1(    ret) ASL(0)
    01638 #endif
    01639         Squ_Column0(0, 1)
    01640         Squ_Column1(1, 1)
    01641         Squ_Column0(2, 2)
    01642         Squ_Column1(3, 1)
    01643         Squ_Column0(4, 1)
    01644         Squ_End(4)
    01645 }
    01646
    01647 void SSE2_Square16(word *C, const word *A)
    01648 {
    01649         Squ_Begin(8)
    01650 #ifndef __GNUC__
    01651         ASJ(    jmp,    0, f)
    01652         Squ_Acc(4) Squ_Acc(3) Squ_Acc(2)
    01653         AS1(    ret) ASL(0)
    01654 #endif
    01655         Squ_Column0(0, 1)
    01656         Squ_Column1(1, 1)
    01657         Squ_Column0(2, 2)
    01658         Squ_Column1(3, 2)
    01659         Squ_Column0(4, 3)
    01660         Squ_Column1(5, 3)
    01661         Squ_Column0(6, 4)
    01662         Squ_Column1(7, 3)
    01663         Squ_Column0(8, 3)
    01664         Squ_Column1(9, 2)
    01665         Squ_Column0(10, 2)
    01666         Squ_Column1(11, 1)
    01667         Squ_Column0(12, 1)
    01668         Squ_End(8)
    01669 }
    01670
    01671 void SSE2_Square32(word *C, const word *A)
    01672 {
    01673         Squ_Begin(16)
    01674         ASJ(    jmp,    0, f)
    01675         Squ_Acc(8) Squ_Acc(7) Squ_Acc(6) Squ_Acc(5) Squ_Acc(4) Squ_Acc(3) Squ_Acc(2)
    01676         AS1(    ret) ASL(0)
    01677         Squ_Column0(0, 1)
    01678         Squ_Column1(1, 1)
    01679         Squ_Column0(2, 2)
    01680         Squ_Column1(3, 2)
    01681         Squ_Column0(4, 3)
    01682         Squ_Column1(5, 3)
    01683         Squ_Column0(6, 4)
    01684         Squ_Column1(7, 4)
    01685         Squ_Column0(8, 5)
    01686         Squ_Column1(9, 5)
    01687         Squ_Column0(10, 6)
    01688         Squ_Column1(11, 6)
    01689         Squ_Column0(12, 7)
    01690         Squ_Column1(13, 7)
    01691         Squ_Column0(14, 8)
    01692         Squ_Column1(15, 7)
    01693         Squ_Column0(16, 7)
    01694         Squ_Column1(17, 6)
    01695         Squ_Column0(18, 6)
    01696         Squ_Column1(19, 5)
    01697         Squ_Column0(20, 5)
    01698         Squ_Column1(21, 4)
    01699         Squ_Column0(22, 4)
    01700         Squ_Column1(23, 3)
    01701         Squ_Column0(24, 3)
    01702         Squ_Column1(25, 2)
    01703         Squ_Column0(26, 2)
    01704         Squ_Column1(27, 1)
    01705         Squ_Column0(28, 1)
    01706         Squ_End(16)
    01707 }
    01708
    01709 void SSE2_Multiply4(word *C, const word *A, const word *B)
    01710 {
    01711         Mul_Begin(2)
    01712 #ifndef __GNUC__
    01713         ASJ(    jmp,    0, f)
    01714         Mul_Acc(2)
    01715         AS1(    ret) ASL(0)
    01716 #endif
    01717         Mul_Column0(0, 2)
    01718         Mul_End(2)
    01719 }
    01720
    01721 void SSE2_Multiply8(word *C, const word *A, const word *B)
    01722 {
    01723         Mul_Begin(4)
    01724 #ifndef __GNUC__
    01725         ASJ(    jmp,    0, f)
    01726         Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01727         AS1(    ret) ASL(0)
    01728 #endif
    01729         Mul_Column0(0, 2)
    01730         Mul_Column1(1, 3)
    01731         Mul_Column0(2, 4)
    01732         Mul_Column1(3, 3)
    01733         Mul_Column0(4, 2)
    01734         Mul_End(4)
    01735 }
    01736
    01737 void SSE2_Multiply16(word *C, const word *A, const word *B)
    01738 {
    01739         Mul_Begin(8)
    01740 #ifndef __GNUC__
    01741         ASJ(    jmp,    0, f)
    01742         Mul_Acc(8) Mul_Acc(7) Mul_Acc(6) Mul_Acc(5) Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01743         AS1(    ret) ASL(0)
    01744 #endif
    01745         Mul_Column0(0, 2)
    01746         Mul_Column1(1, 3)
    01747         Mul_Column0(2, 4)
    01748         Mul_Column1(3, 5)
    01749         Mul_Column0(4, 6)
    01750         Mul_Column1(5, 7)
    01751         Mul_Column0(6, 8)
    01752         Mul_Column1(7, 7)
    01753         Mul_Column0(8, 6)
    01754         Mul_Column1(9, 5)
    01755         Mul_Column0(10, 4)
    01756         Mul_Column1(11, 3)
    01757         Mul_Column0(12, 2)
    01758         Mul_End(8)
    01759 }
    01760
    01761 void SSE2_Multiply32(word *C, const word *A, const word *B)
    01762 {
    01763         Mul_Begin(16)
    01764         ASJ(    jmp,    0, f)
    01765         Mul_Acc(16) Mul_Acc(15) Mul_Acc(14) Mul_Acc(13) Mul_Acc(12) Mul_Acc(11) Mul_Acc(10) Mul_Acc(9) Mul_Acc(8) Mul_Acc(7) Mul_Acc(6) Mul_Acc(5) Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01766         AS1(    ret) ASL(0)
    01767         Mul_Column0(0, 2)
    01768         Mul_Column1(1, 3)
    01769         Mul_Column0(2, 4)
    01770         Mul_Column1(3, 5)
    01771         Mul_Column0(4, 6)
    01772         Mul_Column1(5, 7)
    01773         Mul_Column0(6, 8)
    01774         Mul_Column1(7, 9)
    01775         Mul_Column0(8, 10)
    01776         Mul_Column1(9, 11)
    01777         Mul_Column0(10, 12)
    01778         Mul_Column1(11, 13)
    01779         Mul_Column0(12, 14)
    01780         Mul_Column1(13, 15)
    01781         Mul_Column0(14, 16)
    01782         Mul_Column1(15, 15)
    01783         Mul_Column0(16, 14)
    01784         Mul_Column1(17, 13)
    01785         Mul_Column0(18, 12)
    01786         Mul_Column1(19, 11)
    01787         Mul_Column0(20, 10)
    01788         Mul_Column1(21, 9)
    01789         Mul_Column0(22, 8)
    01790         Mul_Column1(23, 7)
    01791         Mul_Column0(24, 6)
    01792         Mul_Column1(25, 5)
    01793         Mul_Column0(26, 4)
    01794         Mul_Column1(27, 3)
    01795         Mul_Column0(28, 2)
    01796         Mul_End(16)
    01797 }
    01798
    01799 void SSE2_MultiplyBottom4(word *C, const word *A, const word *B)
    01800 {
    01801         Mul_Begin(2)
    01802         Bot_SaveAcc(0) Bot_Acc(2)
    01803         Bot_End(2)
    01804 }
    01805
    01806 void SSE2_MultiplyBottom8(word *C, const word *A, const word *B)
    01807 {
    01808         Mul_Begin(4)
    01809 #ifndef __GNUC__
    01810         ASJ(    jmp,    0, f)
    01811         Mul_Acc(3) Mul_Acc(2)
    01812         AS1(    ret) ASL(0)
    01813 #endif
    01814         Mul_Column0(0, 2)
    01815         Mul_Column1(1, 3)
    01816         Bot_SaveAcc(2) Bot_Acc(4) Bot_Acc(3) Bot_Acc(2)
    01817         Bot_End(4)
    01818 }
    01819
    01820 void SSE2_MultiplyBottom16(word *C, const word *A, const word *B)
    01821 {
    01822         Mul_Begin(8)
    01823 #ifndef __GNUC__
    01824         ASJ(    jmp,    0, f)
    01825         Mul_Acc(7) Mul_Acc(6) Mul_Acc(5) Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01826         AS1(    ret) ASL(0)
    01827 #endif
    01828         Mul_Column0(0, 2)
    01829         Mul_Column1(1, 3)
    01830         Mul_Column0(2, 4)
    01831         Mul_Column1(3, 5)
    01832         Mul_Column0(4, 6)
    01833         Mul_Column1(5, 7)
    01834         Bot_SaveAcc(6) Bot_Acc(8) Bot_Acc(7) Bot_Acc(6) Bot_Acc(5) Bot_Acc(4) Bot_Acc(3) Bot_Acc(2)
    01835         Bot_End(8)
    01836 }
    01837
    01838 void SSE2_MultiplyBottom32(word *C, const word *A, const word *B)
    01839 {
    01840         Mul_Begin(16)
    01841 #ifndef __GNUC__
    01842         ASJ(    jmp,    0, f)
    01843         Mul_Acc(15) Mul_Acc(14) Mul_Acc(13) Mul_Acc(12) Mul_Acc(11) Mul_Acc(10) Mul_Acc(9) Mul_Acc(8) Mul_Acc(7) Mul_Acc(6) Mul_Acc(5) Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01844         AS1(    ret) ASL(0)
    01845 #endif
    01846         Mul_Column0(0, 2)
    01847         Mul_Column1(1, 3)
    01848         Mul_Column0(2, 4)
    01849         Mul_Column1(3, 5)
    01850         Mul_Column0(4, 6)
    01851         Mul_Column1(5, 7)
    01852         Mul_Column0(6, 8)
    01853         Mul_Column1(7, 9)
    01854         Mul_Column0(8, 10)
    01855         Mul_Column1(9, 11)
    01856         Mul_Column0(10, 12)
    01857         Mul_Column1(11, 13)
    01858         Mul_Column0(12, 14)
    01859         Mul_Column1(13, 15)
    01860         Bot_SaveAcc(14) Bot_Acc(16) Bot_Acc(15) Bot_Acc(14) Bot_Acc(13) Bot_Acc(12) Bot_Acc(11) Bot_Acc(10) Bot_Acc(9) Bot_Acc(8) Bot_Acc(7) Bot_Acc(6) Bot_Acc(5) Bot_Acc(4) Bot_Acc(3) Bot_Acc(2)
    01861         Bot_End(16)
    01862 }
    01863
    01864 void SSE2_MultiplyTop8(word *C, const word *A, const word *B, word L)
    01865 {
    01866         Top_Begin(4)
    01867         Top_Acc(3) Top_Acc(2) Top_Acc(1)
    01868 #ifndef __GNUC__
    01869         ASJ(    jmp,    0, f)
    01870         Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01871         AS1(    ret) ASL(0)
    01872 #endif
    01873         Top_Column0(4)
    01874         Top_Column1(3)
    01875         Mul_Column0(0, 2)
    01876         Top_End(2)
    01877 }
    01878
    01879 void SSE2_MultiplyTop16(word *C, const word *A, const word *B, word L)
    01880 {
    01881         Top_Begin(8)
    01882         Top_Acc(7) Top_Acc(6) Top_Acc(5) Top_Acc(4) Top_Acc(3) Top_Acc(2) Top_Acc(1)
    01883 #ifndef __GNUC__
    01884         ASJ(    jmp,    0, f)
    01885         Mul_Acc(8) Mul_Acc(7) Mul_Acc(6) Mul_Acc(5) Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01886         AS1(    ret) ASL(0)
    01887 #endif
    01888         Top_Column0(8)
    01889         Top_Column1(7)
    01890         Mul_Column0(0, 6)
    01891         Mul_Column1(1, 5)
    01892         Mul_Column0(2, 4)
    01893         Mul_Column1(3, 3)
    01894         Mul_Column0(4, 2)
    01895         Top_End(4)
    01896 }
    01897
    01898 void SSE2_MultiplyTop32(word *C, const word *A, const word *B, word L)
    01899 {
    01900         Top_Begin(16)
    01901         Top_Acc(15) Top_Acc(14) Top_Acc(13) Top_Acc(12) Top_Acc(11) Top_Acc(10) Top_Acc(9) Top_Acc(8) Top_Acc(7) Top_Acc(6) Top_Acc(5) Top_Acc(4) Top_Acc(3) Top_Acc(2) Top_Acc(1)
    01902 #ifndef __GNUC__
    01903         ASJ(    jmp,    0, f)
    01904         Mul_Acc(16) Mul_Acc(15) Mul_Acc(14) Mul_Acc(13) Mul_Acc(12) Mul_Acc(11) Mul_Acc(10) Mul_Acc(9) Mul_Acc(8) Mul_Acc(7) Mul_Acc(6) Mul_Acc(5) Mul_Acc(4) Mul_Acc(3) Mul_Acc(2)
    01905         AS1(    ret) ASL(0)
    01906 #endif
    01907         Top_Column0(16)
    01908         Top_Column1(15)
    01909         Mul_Column0(0, 14)
    01910         Mul_Column1(1, 13)
    01911         Mul_Column0(2, 12)
    01912         Mul_Column1(3, 11)
    01913         Mul_Column0(4, 10)
    01914         Mul_Column1(5, 9)
    01915         Mul_Column0(6, 8)
    01916         Mul_Column1(7, 7)
    01917         Mul_Column0(8, 6)
    01918         Mul_Column1(9, 5)
    01919         Mul_Column0(10, 4)
    01920         Mul_Column1(11, 3)
    01921         Mul_Column0(12, 2)
    01922         Top_End(8)
    01923 }
    01924
    01925 #endif  // #if CRYPTOPP_INTEGER_SSE2
    01926
    01927 // ********************************************************
    01928
    01929 typedef int (CRYPTOPP_FASTCALL * PAdd)(size_t N, word *C, const word *A, const word *B);
    01930 typedef void (* PMul)(word *C, const word *A, const word *B);
    01931 typedef void (* PSqu)(word *C, const word *A);
    01932 typedef void (* PMulTop)(word *C, const word *A, const word *B, word L);
    01933
    01934 #if CRYPTOPP_INTEGER_SSE2
    01935 static PAdd s_pAdd = &Baseline_Add, s_pSub = &Baseline_Sub;
    01936 static size_t s_recursionLimit = 8;
    01937 #else
    01938 static const size_t s_recursionLimit = 16;
    01939 #endif
    01940
    01941 static PMul s_pMul[9], s_pBot[9];
    01942 static PSqu s_pSqu[9];
    01943 static PMulTop s_pTop[9];
    01944
    01945 static void SetFunctionPointers()
    01946 {
    01947         s_pMul[0] = &Baseline_Multiply2;
    01948         s_pBot[0] = &Baseline_MultiplyBottom2;
    01949         s_pSqu[0] = &Baseline_Square2;
    01950         s_pTop[0] = &Baseline_MultiplyTop2;
    01951         s_pTop[1] = &Baseline_MultiplyTop4;
    01952
    01953 #if CRYPTOPP_INTEGER_SSE2
    01954         if (HasSSE2())
    01955         {
    01956                 if (IsP4())
    01957                 {
    01958                         s_pAdd = &SSE2_Add;
    01959                         s_pSub = &SSE2_Sub;
    01960                 }
    01961
    01962                 s_recursionLimit = 32;
    01963
    01964                 s_pMul[1] = &SSE2_Multiply4;
    01965                 s_pMul[2] = &SSE2_Multiply8;
    01966                 s_pMul[4] = &SSE2_Multiply16;
    01967                 s_pMul[8] = &SSE2_Multiply32;
    01968
    01969                 s_pBot[1] = &SSE2_MultiplyBottom4;
    01970                 s_pBot[2] = &SSE2_MultiplyBottom8;
    01971                 s_pBot[4] = &SSE2_MultiplyBottom16;
    01972                 s_pBot[8] = &SSE2_MultiplyBottom32;
    01973
    01974                 s_pSqu[1] = &SSE2_Square4;
    01975                 s_pSqu[2] = &SSE2_Square8;
    01976                 s_pSqu[4] = &SSE2_Square16;
    01977                 s_pSqu[8] = &SSE2_Square32;
    01978
    01979                 s_pTop[2] = &SSE2_MultiplyTop8;
    01980                 s_pTop[4] = &SSE2_MultiplyTop16;
    01981                 s_pTop[8] = &SSE2_MultiplyTop32;
    01982         }
    01983         else
    01984 #endif
    01985         {
    01986                 s_pMul[1] = &Baseline_Multiply4;
    01987                 s_pMul[2] = &Baseline_Multiply8;
    01988
    01989                 s_pBot[1] = &Baseline_MultiplyBottom4;
    01990                 s_pBot[2] = &Baseline_MultiplyBottom8;
    01991
    01992                 s_pSqu[1] = &Baseline_Square4;
    01993                 s_pSqu[2] = &Baseline_Square8;
    01994
    01995                 s_pTop[2] = &Baseline_MultiplyTop8;
    01996
    01997 #if     !CRYPTOPP_INTEGER_SSE2
    01998                 s_pMul[4] = &Baseline_Multiply16;
    01999                 s_pBot[4] = &Baseline_MultiplyBottom16;
    02000                 s_pSqu[4] = &Baseline_Square16;
    02001                 s_pTop[4] = &Baseline_MultiplyTop16;
    02002 #endif
    02003         }
    02004 }
    02005
    02006 inline int Add(word *C, const word *A, const word *B, size_t N)
    02007 {
    02008 #if CRYPTOPP_INTEGER_SSE2
    02009         return s_pAdd(N, C, A, B);
    02010 #else
    02011         return Baseline_Add(N, C, A, B);
    02012 #endif
    02013 }
    02014
    02015 inline int Subtract(word *C, const word *A, const word *B, size_t N)
    02016 {
    02017 #if CRYPTOPP_INTEGER_SSE2
    02018         return s_pSub(N, C, A, B);
    02019 #else
    02020         return Baseline_Sub(N, C, A, B);
    02021 #endif
    02022 }
    02023
    02024 // ********************************************************
    02025
    02026
    02027 #define A0              A
    02028 #define A1              (A+N2)
    02029 #define B0              B
    02030 #define B1              (B+N2)
    02031
    02032 #define T0              T
    02033 #define T1              (T+N2)
    02034 #define T2              (T+N)
    02035 #define T3              (T+N+N2)
    02036
    02037 #define R0              R
    02038 #define R1              (R+N2)
    02039 #define R2              (R+N)
    02040 #define R3              (R+N+N2)
    02041
    02042 // R[2*N] - result = A*B
    02043 // T[2*N] - temporary work space
    02044 // A[N] --- multiplier
    02045 // B[N] --- multiplicant
    02046
    02047 void RecursiveMultiply(word *R, word *T, const word *A, const word *B, size_t N)
    02048 {
    02049         assert(N>=2 && N%2==0);
    02050
    02051         if (N <= s_recursionLimit)
    02052                 s_pMul[N/4](R, A, B);
    02053         else
    02054         {
    02055                 const size_t N2 = N/2;
    02056
    02057                 size_t AN2 = Compare(A0, A1, N2) > 0 ?  0 : N2;
    02058                 Subtract(R0, A + AN2, A + (N2 ^ AN2), N2);
    02059
    02060                 size_t BN2 = Compare(B0, B1, N2) > 0 ?  0 : N2;
    02061                 Subtract(R1, B + BN2, B + (N2 ^ BN2), N2);
    02062
    02063                 RecursiveMultiply(R2, T2, A1, B1, N2);
    02064                 RecursiveMultiply(T0, T2, R0, R1, N2);
    02065                 RecursiveMultiply(R0, T2, A0, B0, N2);
    02066
    02067                 // now T[01] holds (A1-A0)*(B0-B1), R[01] holds A0*B0, R[23] holds A1*B1
    02068
    02069                 int c2 = Add(R2, R2, R1, N2);
    02070                 int c3 = c2;
    02071                 c2 += Add(R1, R2, R0, N2);
    02072                 c3 += Add(R2, R2, R3, N2);
    02073
    02074                 if (AN2 == BN2)
    02075                         c3 -= Subtract(R1, R1, T0, N);
    02076                 else
    02077                         c3 += Add(R1, R1, T0, N);
    02078
    02079                 c3 += Increment(R2, N2, c2);
    02080                 assert (c3 >= 0 && c3 <= 2);
    02081                 Increment(R3, N2, c3);
    02082         }
    02083 }
    02084
    02085 // R[2*N] - result = A*A
    02086 // T[2*N] - temporary work space
    02087 // A[N] --- number to be squared
    02088
    02089 void RecursiveSquare(word *R, word *T, const word *A, size_t N)
    02090 {
    02091         assert(N && N%2==0);
    02092
    02093         if (N <= s_recursionLimit)
    02094                 s_pSqu[N/4](R, A);
    02095         else
    02096         {
    02097                 const size_t N2 = N/2;
    02098
    02099                 RecursiveSquare(R0, T2, A0, N2);
    02100                 RecursiveSquare(R2, T2, A1, N2);
    02101                 RecursiveMultiply(T0, T2, A0, A1, N2);
    02102
    02103                 int carry = Add(R1, R1, T0, N);
    02104                 carry += Add(R1, R1, T0, N);
    02105                 Increment(R3, N2, carry);
    02106         }
    02107 }
    02108
    02109 // R[N] - bottom half of A*B
    02110 // T[3*N/2] - temporary work space
    02111 // A[N] - multiplier
    02112 // B[N] - multiplicant
    02113
    02114 void RecursiveMultiplyBottom(word *R, word *T, const word *A, const word *B, size_t N)
    02115 {
    02116         assert(N>=2 && N%2==0);
    02117
    02118         if (N <= s_recursionLimit)
    02119                 s_pBot[N/4](R, A, B);
    02120         else
    02121         {
    02122                 const size_t N2 = N/2;
    02123
    02124                 RecursiveMultiply(R, T, A0, B0, N2);
    02125                 RecursiveMultiplyBottom(T0, T1, A1, B0, N2);
    02126                 Add(R1, R1, T0, N2);
    02127                 RecursiveMultiplyBottom(T0, T1, A0, B1, N2);
    02128                 Add(R1, R1, T0, N2);
    02129         }
    02130 }
    02131
    02132 // R[N] --- upper half of A*B
    02133 // T[2*N] - temporary work space
    02134 // L[N] --- lower half of A*B
    02135 // A[N] --- multiplier
    02136 // B[N] --- multiplicant
    02137
    02138 void MultiplyTop(word *R, word *T, const word *L, const word *A, const word *B, size_t N)
    02139 {
    02140         assert(N>=2 && N%2==0);
    02141
    02142         if (N <= s_recursionLimit)
    02143                 s_pTop[N/4](R, A, B, L[N-1]);
    02144         else
    02145         {
    02146                 const size_t N2 = N/2;
    02147
    02148                 size_t AN2 = Compare(A0, A1, N2) > 0 ?  0 : N2;
    02149                 Subtract(R0, A + AN2, A + (N2 ^ AN2), N2);
    02150
    02151                 size_t BN2 = Compare(B0, B1, N2) > 0 ?  0 : N2;
    02152                 Subtract(R1, B + BN2, B + (N2 ^ BN2), N2);
    02153
    02154                 RecursiveMultiply(T0, T2, R0, R1, N2);
    02155                 RecursiveMultiply(R0, T2, A1, B1, N2);
    02156
    02157                 // now T[01] holds (A1-A0)*(B0-B1) = A1*B0+A0*B1-A1*B1-A0*B0, R[01] holds A1*B1
    02158
    02159                 int t, c3;
    02160                 int c2 = Subtract(T2, L+N2, L, N2);
    02161
    02162                 if (AN2 == BN2)
    02163                 {
    02164                         c2 -= Add(T2, T2, T0, N2);
    02165                         t = (Compare(T2, R0, N2) == -1);
    02166                         c3 = t - Subtract(T2, T2, T1, N2);
    02167                 }
    02168                 else
    02169                 {
    02170                         c2 += Subtract(T2, T2, T0, N2);
    02171                         t = (Compare(T2, R0, N2) == -1);
    02172                         c3 = t + Add(T2, T2, T1, N2);
    02173                 }
    02174
    02175                 c2 += t;
    02176                 if (c2 >= 0)
    02177                         c3 += Increment(T2, N2, c2);
    02178                 else
    02179                         c3 -= Decrement(T2, N2, -c2);
    02180                 c3 += Add(R0, T2, R1, N2);
    02181
    02182                 assert (c3 >= 0 && c3 <= 2);
    02183                 Increment(R1, N2, c3);
    02184         }
    02185 }
    02186
    02187 inline void Multiply(word *R, word *T, const word *A, const word *B, size_t N)
    02188 {
    02189         RecursiveMultiply(R, T, A, B, N);
    02190 }
    02191
    02192 inline void Square(word *R, word *T, const word *A, size_t N)
    02193 {
    02194         RecursiveSquare(R, T, A, N);
    02195 }
    02196
    02197 inline void MultiplyBottom(word *R, word *T, const word *A, const word *B, size_t N)
    02198 {
    02199         RecursiveMultiplyBottom(R, T, A, B, N);
    02200 }
    02201
    02202 // R[NA+NB] - result = A*B
    02203 // T[NA+NB] - temporary work space
    02204 // A[NA] ---- multiplier
    02205 // B[NB] ---- multiplicant
    02206
    02207 void AsymmetricMultiply(word *R, word *T, const word *A, size_t NA, const word *B, size_t NB)
    02208 {
    02209         if (NA == NB)
    02210         {
    02211                 if (A == B)
    02212                         Square(R, T, A, NA);
    02213                 else
    02214                         Multiply(R, T, A, B, NA);
    02215
    02216                 return;
    02217         }
    02218
    02219         if (NA > NB)
    02220         {
    02221                 std::swap(A, B);
    02222                 std::swap(NA, NB);
    02223         }
    02224
    02225         assert(NB % NA == 0);
    02226
    02227         if (NA==2 && !A[1])
    02228         {
    02229                 switch (A[0])
    02230                 {
    02231                 case 0:
    02232                         SetWords(R, 0, NB+2);
    02233                         return;
    02234                 case 1:
    02235                         CopyWords(R, B, NB);
    02236                         R[NB] = R[NB+1] = 0;
    02237                         return;
    02238                 default:
    02239                         R[NB] = LinearMultiply(R, B, A[0], NB);
    02240                         R[NB+1] = 0;
    02241                         return;
    02242                 }
    02243         }
    02244
    02245         size_t i;
    02246         if ((NB/NA)%2 == 0)
    02247         {
    02248                 Multiply(R, T, A, B, NA);
    02249                 CopyWords(T+2*NA, R+NA, NA);
    02250
    02251                 for (i=2*NA; i<NB; i+=2*NA)
    02252                         Multiply(T+NA+i, T, A, B+i, NA);
    02253                 for (i=NA; i<NB; i+=2*NA)
    02254                         Multiply(R+i, T, A, B+i, NA);
    02255         }
    02256         else
    02257         {
    02258                 for (i=0; i<NB; i+=2*NA)
    02259                         Multiply(R+i, T, A, B+i, NA);
    02260                 for (i=NA; i<NB; i+=2*NA)
    02261                         Multiply(T+NA+i, T, A, B+i, NA);
    02262         }
    02263
    02264         if (Add(R+NA, R+NA, T+2*NA, NB-NA))
    02265                 Increment(R+NB, NA);
    02266 }
    02267
    02268 // R[N] ----- result = A inverse mod 2**(WORD_BITS*N)
    02269 // T[3*N/2] - temporary work space
    02270 // A[N] ----- an odd number as input
    02271
    02272 void RecursiveInverseModPower2(word *R, word *T, const word *A, size_t N)
    02273 {
    02274         if (N==2)
    02275         {
    02276                 T[0] = AtomicInverseModPower2(A[0]);
    02277                 T[1] = 0;
    02278                 s_pBot[0](T+2, T, A);
    02279                 TwosComplement(T+2, 2);
    02280                 Increment(T+2, 2, 2);
    02281                 s_pBot[0](R, T, T+2);
    02282         }
    02283         else
    02284         {
    02285                 const size_t N2 = N/2;
    02286                 RecursiveInverseModPower2(R0, T0, A0, N2);
    02287                 T0[0] = 1;
    02288                 SetWords(T0+1, 0, N2-1);
    02289                 MultiplyTop(R1, T1, T0, R0, A0, N2);
    02290                 MultiplyBottom(T0, T1, R0, A1, N2);
    02291                 Add(T0, R1, T0, N2);
    02292                 TwosComplement(T0, N2);
    02293                 MultiplyBottom(R1, T1, R0, T0, N2);
    02294         }
    02295 }
    02296
    02297 // R[N] --- result = X/(2**(WORD_BITS*N)) mod M
    02298 // T[3*N] - temporary work space
    02299 // X[2*N] - number to be reduced
    02300 // M[N] --- modulus
    02301 // U[N] --- multiplicative inverse of M mod 2**(WORD_BITS*N)
    02302
    02303 void MontgomeryReduce(word *R, word *T, word *X, const word *M, const word *U, size_t N)
    02304 {
    02305 #if 1
    02306         MultiplyBottom(R, T, X, U, N);
    02307         MultiplyTop(T, T+N, X, R, M, N);
    02308         word borrow = Subtract(T, X+N, T, N);
    02309         // defend against timing attack by doing this Add even when not needed
    02310         word carry = Add(T+N, T, M, N);
    02311         assert(carry || !borrow);
    02312         CopyWords(R, T + (borrow ? N : 0), N);
    02313 #elif 0
    02314         const word u = 0-U[0];
    02315         Declare2Words(p)
    02316         for (size_t i=0; i<N; i++)
    02317         {
    02318                 const word t = u * X[i];
    02319                 word c = 0;
    02320                 for (size_t j=0; j<N; j+=2)
    02321                 {
    02322                         MultiplyWords(p, t, M[j]);
    02323                         Acc2WordsBy1(p, X[i+j]);
    02324                         Acc2WordsBy1(p, c);
    02325                         X[i+j] = LowWord(p);
    02326                         c = HighWord(p);
    02327                         MultiplyWords(p, t, M[j+1]);
    02328                         Acc2WordsBy1(p, X[i+j+1]);
    02329                         Acc2WordsBy1(p, c);
    02330                         X[i+j+1] = LowWord(p);
    02331                         c = HighWord(p);
    02332                 }
    02333
    02334                 if (Increment(X+N+i, N-i, c))
    02335                         while (!Subtract(X+N, X+N, M, N)) {}
    02336         }
    02337
    02338         memcpy(R, X+N, N*WORD_SIZE);
    02339 #else
    02340         __m64 u = _mm_cvtsi32_si64(0-U[0]), p;
    02341         for (size_t i=0; i<N; i++)
    02342         {
    02343                 __m64 t = _mm_cvtsi32_si64(X[i]);
    02344                 t = _mm_mul_su32(t, u);
    02345                 __m64 c = _mm_setzero_si64();
    02346                 for (size_t j=0; j<N; j+=2)
    02347                 {
    02348                         p = _mm_mul_su32(t, _mm_cvtsi32_si64(M[j]));
    02349                         p = _mm_add_si64(p, _mm_cvtsi32_si64(X[i+j]));
    02350                         c = _mm_add_si64(c, p);
    02351                         X[i+j] = _mm_cvtsi64_si32(c);
    02352                         c = _mm_srli_si64(c, 32);
    02353                         p = _mm_mul_su32(t, _mm_cvtsi32_si64(M[j+1]));
    02354                         p = _mm_add_si64(p, _mm_cvtsi32_si64(X[i+j+1]));
    02355                         c = _mm_add_si64(c, p);
    02356                         X[i+j+1] = _mm_cvtsi64_si32(c);
    02357                         c = _mm_srli_si64(c, 32);
    02358                 }
    02359
    02360                 if (Increment(X+N+i, N-i, _mm_cvtsi64_si32(c)))
    02361                         while (!Subtract(X+N, X+N, M, N)) {}
    02362         }
    02363
    02364         memcpy(R, X+N, N*WORD_SIZE);
    02365         _mm_empty();
    02366 #endif
    02367 }
    02368
    02369 // R[N] --- result = X/(2**(WORD_BITS*N/2)) mod M
    02370 // T[2*N] - temporary work space
    02371 // X[2*N] - number to be reduced
    02372 // M[N] --- modulus
    02373 // U[N/2] - multiplicative inverse of M mod 2**(WORD_BITS*N/2)
    02374 // V[N] --- 2**(WORD_BITS*3*N/2) mod M
    02375
    02376 void HalfMontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, const word *V, size_t N)
    02377 {
    02378         assert(N%2==0 && N>=4);
    02379
    02380 #define M0              M
    02381 #define M1              (M+N2)
    02382 #define V0              V
    02383 #define V1              (V+N2)
    02384
    02385 #define X0              X
    02386 #define X1              (X+N2)
    02387 #define X2              (X+N)
    02388 #define X3              (X+N+N2)
    02389
    02390         const size_t N2 = N/2;
    02391         Multiply(T0, T2, V0, X3, N2);
    02392         int c2 = Add(T0, T0, X0, N);
    02393         MultiplyBottom(T3, T2, T0, U, N2);
    02394         MultiplyTop(T2, R, T0, T3, M0, N2);
    02395         c2 -= Subtract(T2, T1, T2, N2);
    02396         Multiply(T0, R, T3, M1, N2);
    02397         c2 -= Subtract(T0, T2, T0, N2);
    02398         int c3 = -(int)Subtract(T1, X2, T1, N2);
    02399         Multiply(R0, T2, V1, X3, N2);
    02400         c3 += Add(R, R, T, N);
    02401
    02402         if (c2>0)
    02403                 c3 += Increment(R1, N2);
    02404         else if (c2<0)
    02405                 c3 -= Decrement(R1, N2, -c2);
    02406
    02407         assert(c3>=-1 && c3<=1);
    02408         if (c3>0)
    02409                 Subtract(R, R, M, N);
    02410         else if (c3<0)
    02411                 Add(R, R, M, N);
    02412
    02413 #undef M0
    02414 #undef M1
    02415 #undef V0
    02416 #undef V1
    02417
    02418 #undef X0
    02419 #undef X1
    02420 #undef X2
    02421 #undef X3
    02422 }
    02423
    02424 #undef A0
    02425 #undef A1
    02426 #undef B0
    02427 #undef B1
    02428
    02429 #undef T0
    02430 #undef T1
    02431 #undef T2
    02432 #undef T3
    02433
    02434 #undef R0
    02435 #undef R1
    02436 #undef R2
    02437 #undef R3
    02438
    02439 /*
    02440 // do a 3 word by 2 word divide, returns quotient and leaves remainder in A
    02441 static word SubatomicDivide(word *A, word B0, word B1)
    02442 {
    02443         // assert {A[2],A[1]} < {B1,B0}, so quotient can fit in a word
    02444         assert(A[2] < B1 || (A[2]==B1 && A[1] < B0));
    02445
    02446         // estimate the quotient: do a 2 word by 1 word divide
    02447         word Q;
    02448         if (B1+1 == 0)
    02449                 Q = A[2];
    02450         else
    02451                 Q = DWord(A[1], A[2]).DividedBy(B1+1);
    02452
    02453         // now subtract Q*B from A
    02454         DWord p = DWord::Multiply(B0, Q);
    02455         DWord u = (DWord) A[0] - p.GetLowHalf();
    02456         A[0] = u.GetLowHalf();
    02457         u = (DWord) A[1] - p.GetHighHalf() - u.GetHighHalfAsBorrow() - DWord::Multiply(B1, Q);
    02458         A[1] = u.GetLowHalf();
    02459         A[2] += u.GetHighHalf();
    02460
    02461         // Q <= actual quotient, so fix it
    02462         while (A[2] || A[1] > B1 || (A[1]==B1 && A[0]>=B0))
    02463         {
    02464                 u = (DWord) A[0] - B0;
    02465                 A[0] = u.GetLowHalf();
    02466                 u = (DWord) A[1] - B1 - u.GetHighHalfAsBorrow();
    02467                 A[1] = u.GetLowHalf();
    02468                 A[2] += u.GetHighHalf();
    02469                 Q++;
    02470                 assert(Q);      // shouldn't overflow
    02471         }
    02472
    02473         return Q;
    02474 }
    02475
    02476 // do a 4 word by 2 word divide, returns 2 word quotient in Q0 and Q1
    02477 static inline void AtomicDivide(word *Q, const word *A, const word *B)
    02478 {
    02479         if (!B[0] && !B[1]) // if divisor is 0, we assume divisor==2**(2*WORD_BITS)
    02480         {
    02481                 Q[0] = A[2];
    02482                 Q[1] = A[3];
    02483         }
    02484         else
    02485         {
    02486                 word T[4];
    02487                 T[0] = A[0]; T[1] = A[1]; T[2] = A[2]; T[3] = A[3];
    02488                 Q[1] = SubatomicDivide(T+1, B[0], B[1]);
    02489                 Q[0] = SubatomicDivide(T, B[0], B[1]);
    02490
    02491 #ifndef NDEBUG
    02492                 // multiply quotient and divisor and add remainder, make sure it equals dividend
    02493                 assert(!T[2] && !T[3] && (T[1] < B[1] || (T[1]==B[1] && T[0]<B[0])));
    02494                 word P[4];
    02495                 LowLevel::Multiply2(P, Q, B);
    02496                 Add(P, P, T, 4);
    02497                 assert(memcmp(P, A, 4*WORD_SIZE)==0);
    02498 #endif
    02499         }
    02500 }
    02501 */
    02502
    02503 static inline void AtomicDivide(word *Q, const word *A, const word *B)
    02504 {
    02505         word T[4];
    02506         DWord q = DivideFourWordsByTwo<word, DWord>(T, DWord(A[0], A[1]), DWord(A[2], A[3]), DWord(B[0], B[1]));
    02507         Q[0] = q.GetLowHalf();
    02508         Q[1] = q.GetHighHalf();
    02509
    02510 #ifndef NDEBUG
    02511         if (B[0] || B[1])
    02512         {
    02513                 // multiply quotient and divisor and add remainder, make sure it equals dividend
    02514                 assert(!T[2] && !T[3] && (T[1] < B[1] || (T[1]==B[1] && T[0]<B[0])));
    02515                 word P[4];
    02516                 s_pMul[0](P, Q, B);
    02517                 Add(P, P, T, 4);
    02518                 assert(memcmp(P, A, 4*WORD_SIZE)==0);
    02519         }
    02520 #endif
    02521 }
    02522
    02523 // for use by Divide(), corrects the underestimated quotient {Q1,Q0}
    02524 static void CorrectQuotientEstimate(word *R, word *T, word *Q, const word *B, size_t N)
    02525 {
    02526         assert(N && N%2==0);
    02527
    02528         AsymmetricMultiply(T, T+N+2, Q, 2, B, N);
    02529
    02530         word borrow = Subtract(R, R, T, N+2);
    02531         assert(!borrow && !R[N+1]);
    02532
    02533         while (R[N] || Compare(R, B, N) >= 0)
    02534         {
    02535                 R[N] -= Subtract(R, R, B, N);
    02536                 Q[1] += (++Q[0]==0);
    02537                 assert(Q[0] || Q[1]); // no overflow
    02538         }
    02539 }
    02540
    02541 // R[NB] -------- remainder = A%B
    02542 // Q[NA-NB+2] --- quotient      = A/B
    02543 // T[NA+3*(NB+2)] - temp work space
    02544 // A[NA] -------- dividend
    02545 // B[NB] -------- divisor
    02546
    02547 void Divide(word *R, word *Q, word *T, const word *A, size_t NA, const word *B, size_t NB)
    02548 {
    02549         assert(NA && NB && NA%2==0 && NB%2==0);
    02550         assert(B[NB-1] || B[NB-2]);
    02551         assert(NB <= NA);
    02552
    02553         // set up temporary work space
    02554         word *const TA=T;
    02555         word *const TB=T+NA+2;
    02556         word *const TP=T+NA+2+NB;
    02557
    02558         // copy B into TB and normalize it so that TB has highest bit set to 1
    02559         unsigned shiftWords = (B[NB-1]==0);
    02560         TB[0] = TB[NB-1] = 0;
    02561         CopyWords(TB+shiftWords, B, NB-shiftWords);
    02562         unsigned shiftBits = WORD_BITS - BitPrecision(TB[NB-1]);
    02563         assert(shiftBits < WORD_BITS);
    02564         ShiftWordsLeftByBits(TB, NB, shiftBits);
    02565
    02566         // copy A into TA and normalize it
    02567         TA[0] = TA[NA] = TA[NA+1] = 0;
    02568         CopyWords(TA+shiftWords, A, NA);
    02569         ShiftWordsLeftByBits(TA, NA+2, shiftBits);
    02570
    02571         if (TA[NA+1]==0 && TA[NA] <= 1)
    02572         {
    02573                 Q[NA-NB+1] = Q[NA-NB] = 0;
    02574                 while (TA[NA] || Compare(TA+NA-NB, TB, NB) >= 0)
    02575                 {
    02576                         TA[NA] -= Subtract(TA+NA-NB, TA+NA-NB, TB, NB);
    02577                         ++Q[NA-NB];
    02578                 }
    02579         }
    02580         else
    02581         {
    02582                 NA+=2;
    02583                 assert(Compare(TA+NA-NB, TB, NB) < 0);
    02584         }
    02585
    02586         word BT[2];
    02587         BT[0] = TB[NB-2] + 1;
    02588         BT[1] = TB[NB-1] + (BT[0]==0);
    02589
    02590         // start reducing TA mod TB, 2 words at a time
    02591         for (size_t i=NA-2; i>=NB; i-=2)
    02592         {
    02593                 AtomicDivide(Q+i-NB, TA+i-2, BT);
    02594                 CorrectQuotientEstimate(TA+i-NB, TP, Q+i-NB, TB, NB);
    02595         }
    02596
    02597         // copy TA into R, and denormalize it
    02598         CopyWords(R, TA+shiftWords, NB);
    02599         ShiftWordsRightByBits(R, NB, shiftBits);
    02600 }
    02601
    02602 static inline size_t EvenWordCount(const word *X, size_t N)
    02603 {
    02604         while (N && X[N-2]==0 && X[N-1]==0)
    02605                 N-=2;
    02606         return N;
    02607 }
    02608
    02609 // return k
    02610 // R[N] --- result = A^(-1) * 2^k mod M
    02611 // T[4*N] - temporary work space
    02612 // A[NA] -- number to take inverse of
    02613 // M[N] --- modulus
    02614
    02615 unsigned int AlmostInverse(word *R, word *T, const word *A, size_t NA, const word *M, size_t N)
    02616 {
    02617         assert(NA<=N && N && N%2==0);
    02618
    02619         word *b = T;
    02620         word *c = T+N;
    02621         word *f = T+2*N;
    02622         word *g = T+3*N;
    02623         size_t bcLen=2, fgLen=EvenWordCount(M, N);
    02624         unsigned int k=0, s=0;
    02625
    02626         SetWords(T, 0, 3*N);
    02627         b[0]=1;
    02628         CopyWords(f, A, NA);
    02629         CopyWords(g, M, N);
    02630
    02631         while (1)
    02632         {
    02633                 word t=f[0];
    02634                 while (!t)
    02635                 {
    02636                         if (EvenWordCount(f, fgLen)==0)
    02637                         {
    02638                                 SetWords(R, 0, N);
    02639                                 return 0;
    02640                         }
    02641
    02642                         ShiftWordsRightByWords(f, fgLen, 1);
    02643                         if (c[bcLen-1]) bcLen+=2;
    02644                         assert(bcLen <= N);
    02645                         ShiftWordsLeftByWords(c, bcLen, 1);
    02646                         k+=WORD_BITS;
    02647                         t=f[0];
    02648                 }
    02649
    02650                 unsigned int i=0;
    02651                 while (t%2 == 0)
    02652                 {
    02653                         t>>=1;
    02654                         i++;
    02655                 }
    02656                 k+=i;
    02657
    02658                 if (t==1 && f[1]==0 && EvenWordCount(f, fgLen)==2)
    02659                 {
    02660                         if (s%2==0)
    02661                                 CopyWords(R, b, N);
    02662                         else
    02663                                 Subtract(R, M, b, N);
    02664                         return k;
    02665                 }
    02666
    02667                 ShiftWordsRightByBits(f, fgLen, i);
    02668                 t=ShiftWordsLeftByBits(c, bcLen, i);
    02669                 if (t)
    02670                 {
    02671                         c[bcLen] = t;
    02672                         bcLen+=2;
    02673                         assert(bcLen <= N);
    02674                 }
    02675
    02676                 if (f[fgLen-2]==0 && g[fgLen-2]==0 && f[fgLen-1]==0 && g[fgLen-1]==0)
    02677                         fgLen-=2;
    02678
    02679                 if (Compare(f, g, fgLen)==-1)
    02680                 {
    02681                         std::swap(f, g);
    02682                         std::swap(b, c);
    02683                         s++;
    02684                 }
    02685
    02686                 Subtract(f, f, g, fgLen);
    02687
    02688                 if (Add(b, b, c, bcLen))
    02689                 {
    02690                         b[bcLen] = 1;
    02691                         bcLen+=2;
    02692                         assert(bcLen <= N);
    02693                 }
    02694         }
    02695 }
    02696
    02697 // R[N] - result = A/(2^k) mod M
    02698 // A[N] - input
    02699 // M[N] - modulus
    02700
    02701 void DivideByPower2Mod(word *R, const word *A, size_t k, const word *M, size_t N)
    02702 {
    02703         CopyWords(R, A, N);
    02704
    02705         while (k--)
    02706         {
    02707                 if (R[0]%2==0)
    02708                         ShiftWordsRightByBits(R, N, 1);
    02709                 else
    02710                 {
    02711                         word carry = Add(R, R, M, N);
    02712                         ShiftWordsRightByBits(R, N, 1);
    02713                         R[N-1] += carry<<(WORD_BITS-1);
    02714                 }
    02715         }
    02716 }
    02717
    02718 // R[N] - result = A*(2^k) mod M
    02719 // A[N] - input
    02720 // M[N] - modulus
    02721
    02722 void MultiplyByPower2Mod(word *R, const word *A, size_t k, const word *M, size_t N)
    02723 {
    02724         CopyWords(R, A, N);
    02725
    02726         while (k--)
    02727                 if (ShiftWordsLeftByBits(R, N, 1) || Compare(R, M, N)>=0)
    02728                         Subtract(R, R, M, N);
    02729 }
    02730
    02731 // ******************************************************************
    02732
    02733 InitializeInteger::InitializeInteger()
    02734 {
    02735         if (!g_pAssignIntToInteger)
    02736         {
    02737                 SetFunctionPointers();
    02738                 g_pAssignIntToInteger = AssignIntToInteger;
    02739         }
    02740 }
    02741
    02742 static const unsigned int RoundupSizeTable[] = {2, 2, 2, 4, 4, 8, 8, 8, 8};
    02743
    02744 static inline size_t RoundupSize(size_t n)
    02745 {
    02746         if (n<=8)
    02747                 return RoundupSizeTable[n];
    02748         else if (n<=16)
    02749                 return 16;
    02750         else if (n<=32)
    02751                 return 32;
    02752         else if (n<=64)
    02753                 return 64;
    02754         else return size_t(1) << BitPrecision(n-1);
    02755 }
    02756
    02757 Integer::Integer()
    02758         : reg(2), sign(POSITIVE)
    02759 {
    02760         reg[0] = reg[1] = 0;
    02761 }
    02762
    02763 Integer::Integer(const Integer& t)
    02764         : reg(RoundupSize(t.WordCount())), sign(t.sign)
    02765 {
    02766         CopyWords(reg, t.reg, reg.size());
    02767 }
    02768
    02769 Integer::Integer(Sign s, lword value)
    02770         : reg(2), sign(s)
    02771 {
    02772         reg[0] = word(value);
    02773         reg[1] = word(SafeRightShift<WORD_BITS>(value));
    02774 }
    02775
    02776 Integer::Integer(signed long value)
    02777         : reg(2)
    02778 {
    02779         if (value >= 0)
    02780                 sign = POSITIVE;
    02781         else
    02782         {
    02783                 sign = NEGATIVE;
    02784                 value = -value;
    02785         }
    02786         reg[0] = word(value);
    02787         reg[1] = word(SafeRightShift<WORD_BITS>((unsigned long)value));
    02788 }
    02789
    02790 Integer::Integer(Sign s, word high, word low)
    02791         : reg(2), sign(s)
    02792 {
    02793         reg[0] = low;
    02794         reg[1] = high;
    02795 }
    02796
    02797 bool Integer::IsConvertableToLong() const
    02798 {
    02799         if (ByteCount() > sizeof(long))
    02800                 return false;
    02801
    02802         unsigned long value = (unsigned long)reg[0];
    02803         value += SafeLeftShift<WORD_BITS, unsigned long>((unsigned long)reg[1]);
    02804
    02805         if (sign==POSITIVE)
    02806                 return (signed long)value >= 0;
    02807         else
    02808                 return -(signed long)value < 0;
    02809 }
    02810
    02811 signed long Integer::ConvertToLong() const
    02812 {
    02813         assert(IsConvertableToLong());
    02814
    02815         unsigned long value = (unsigned long)reg[0];
    02816         value += SafeLeftShift<WORD_BITS, unsigned long>((unsigned long)reg[1]);
    02817         return sign==POSITIVE ? value : -(signed long)value;
    02818 }
    02819
    02820 Integer::Integer(BufferedTransformation &encodedInteger, size_t byteCount, Signedness s)
    02821 {
    02822         Decode(encodedInteger, byteCount, s);
    02823 }
    02824
    02825 Integer::Integer(const byte *encodedInteger, size_t byteCount, Signedness s)
    02826 {
    02827         Decode(encodedInteger, byteCount, s);
    02828 }
    02829
    02830 Integer::Integer(BufferedTransformation &bt)
    02831 {
    02832         BERDecode(bt);
    02833 }
    02834
    02835 Integer::Integer(RandomNumberGenerator &rng, size_t bitcount)
    02836 {
    02837         Randomize(rng, bitcount);
    02838 }
    02839
    02840 Integer::Integer(RandomNumberGenerator &rng, const Integer &min, const Integer &max, RandomNumberType rnType, const Integer &equiv, const Integer &mod)
    02841 {
    02842         if (!Randomize(rng, min, max, rnType, equiv, mod))
    02843                 throw Integer::RandomNumberNotFound();
    02844 }
    02845
    02846 Integer Integer::Power2(size_t e)
    02847 {
    02848         Integer r((word)0, BitsToWords(e+1));
    02849         r.SetBit(e);
    02850         return r;
    02851 }
    02852
    02853 template <long i>
    02854 struct NewInteger
    02855 {
    02856         Integer * operator()() const
    02857         {
    02858                 return new Integer(i);
    02859         }
    02860 };
    02861
    02862 const Integer &Integer::Zero()
    02863 {
    02864         return Singleton<Integer>().Ref();
    02865 }
    02866
    02867 const Integer &Integer::One()
    02868 {
    02869         return Singleton<Integer, NewInteger<1> >().Ref();
    02870 }
    02871
    02872 const Integer &Integer::Two()
    02873 {
    02874         return Singleton<Integer, NewInteger<2> >().Ref();
    02875 }
    02876
    02877 bool Integer::operator!() const
    02878 {
    02879         return IsNegative() ? false : (reg[0]==0 && WordCount()==0);
    02880 }
    02881
    02882 Integer& Integer::operator=(const Integer& t)
    02883 {
    02884         if (this != &t)
    02885         {
    02886                 if (reg.size() != t.reg.size() || t.reg[t.reg.size()/2] == 0)
    02887                         reg.New(RoundupSize(t.WordCount()));
    02888                 CopyWords(reg, t.reg, reg.size());
    02889                 sign = t.sign;
    02890         }
    02891         return *this;
    02892 }
    02893
    02894 bool Integer::GetBit(size_t n) const
    02895 {
    02896         if (n/WORD_BITS >= reg.size())
    02897                 return 0;
    02898         else
    02899                 return bool((reg[n/WORD_BITS] >> (n % WORD_BITS)) & 1);
    02900 }
    02901
    02902 void Integer::SetBit(size_t n, bool value)
    02903 {
    02904         if (value)
    02905         {
    02906                 reg.CleanGrow(RoundupSize(BitsToWords(n+1)));
    02907                 reg[n/WORD_BITS] |= (word(1) << (n%WORD_BITS));
    02908         }
    02909         else
    02910         {
    02911                 if (n/WORD_BITS < reg.size())
    02912                         reg[n/WORD_BITS] &= ~(word(1) << (n%WORD_BITS));
    02913         }
    02914 }
    02915
    02916 byte Integer::GetByte(size_t n) const
    02917 {
    02918         if (n/WORD_SIZE >= reg.size())
    02919                 return 0;
    02920         else
    02921                 return byte(reg[n/WORD_SIZE] >> ((n%WORD_SIZE)*8));
    02922 }
    02923
    02924 void Integer::SetByte(size_t n, byte value)
    02925 {
    02926         reg.CleanGrow(RoundupSize(BytesToWords(n+1)));
    02927         reg[n/WORD_SIZE] &= ~(word(0xff) << 8*(n%WORD_SIZE));
    02928         reg[n/WORD_SIZE] |= (word(value) << 8*(n%WORD_SIZE));
    02929 }
    02930
    02931 lword Integer::GetBits(size_t i, size_t n) const
    02932 {
    02933         lword v = 0;
    02934         assert(n <= sizeof(v)*8);
    02935         for (unsigned int j=0; j<n; j++)
    02936                 v |= lword(GetBit(i+j)) << j;
    02937         return v;
    02938 }
    02939
    02940 Integer Integer::operator-() const
    02941 {
    02942         Integer result(*this);
    02943         result.Negate();
    02944         return result;
    02945 }
    02946
    02947 Integer Integer::AbsoluteValue() const
    02948 {
    02949         Integer result(*this);
    02950         result.sign = POSITIVE;
    02951         return result;
    02952 }
    02953
    02954 void Integer::swap(Integer &a)
    02955 {
    02956         reg.swap(a.reg);
    02957         std::swap(sign, a.sign);
    02958 }
    02959
    02960 Integer::Integer(word value, size_t length)
    02961         : reg(RoundupSize(length)), sign(POSITIVE)
    02962 {
    02963         reg[0] = value;
    02964         SetWords(reg+1, 0, reg.size()-1);
    02965 }
    02966
    02967 template <class T>
    02968 static Integer StringToInteger(const T *str)
    02969 {
    02970         int radix;
    02971         // GCC workaround
    02972         // std::char_traits<wchar_t>::length() not defined in GCC 3.2 and STLport 4.5.3
    02973         unsigned int length;
    02974         for (length = 0; str[length] != 0; length++) {}
    02975
    02976         Integer v;
    02977
    02978         if (length == 0)
    02979                 return v;
    02980
    02981         switch (str[length-1])
    02982         {
    02983         case 'h':
    02984         case 'H':
    02985                 radix=16;
    02986                 break;
    02987         case 'o':
    02988         case 'O':
    02989                 radix=8;
    02990                 break;
    02991         case 'b':
    02992         case 'B':
    02993                 radix=2;
    02994                 break;
    02995         default:
    02996                 radix=10;
    02997         }
    02998
    02999         if (length > 2 && str[0] == '0' && str[1] == 'x')
    03000                 radix = 16;
    03001
    03002         for (unsigned i=0; i<length; i++)
    03003         {
    03004                 int digit;
    03005
    03006                 if (str[i] >= '0' && str[i] <= '9')
    03007                         digit = str[i] - '0';
    03008                 else if (str[i] >= 'A' && str[i] <= 'F')
    03009                         digit = str[i] - 'A' + 10;
    03010                 else if (str[i] >= 'a' && str[i] <= 'f')
    03011                         digit = str[i] - 'a' + 10;
    03012                 else
    03013                         digit = radix;
    03014
    03015                 if (digit < radix)
    03016                 {
    03017                         v *= radix;
    03018                         v += digit;
    03019                 }
    03020         }
    03021
    03022         if (str[0] == '-')
    03023                 v.Negate();
    03024
    03025         return v;
    03026 }
    03027
    03028 Integer::Integer(const char *str)
    03029         : reg(2), sign(POSITIVE)
    03030 {
    03031         *this = StringToInteger(str);
    03032 }
    03033
    03034 Integer::Integer(const wchar_t *str)
    03035         : reg(2), sign(POSITIVE)
    03036 {
    03037         *this = StringToInteger(str);
    03038 }
    03039
    03040 unsigned int Integer::WordCount() const
    03041 {
    03042         return (unsigned int)CountWords(reg, reg.size());
    03043 }
    03044
    03045 unsigned int Integer::ByteCount() const
    03046 {
    03047         unsigned wordCount = WordCount();
    03048         if (wordCount)
    03049                 return (wordCount-1)*WORD_SIZE + BytePrecision(reg[wordCount-1]);
    03050         else
    03051                 return 0;
    03052 }
    03053
    03054 unsigned int Integer::BitCount() const
    03055 {
    03056         unsigned wordCount = WordCount();
    03057         if (wordCount)
    03058                 return (wordCount-1)*WORD_BITS + BitPrecision(reg[wordCount-1]);
    03059         else
    03060                 return 0;
    03061 }
    03062
    03063 void Integer::Decode(const byte *input, size_t inputLen, Signedness s)
    03064 {
    03065         StringStore store(input, inputLen);
    03066         Decode(store, inputLen, s);
    03067 }
    03068
    03069 void Integer::Decode(BufferedTransformation &bt, size_t inputLen, Signedness s)
    03070 {
    03071         assert(bt.MaxRetrievable() >= inputLen);
    03072
    03073         byte b;
    03074         bt.Peek(b);
    03075         sign = ((s==SIGNED) && (b & 0x80)) ? NEGATIVE : POSITIVE;
    03076
    03077         while (inputLen>0 && (sign==POSITIVE ? b==0 : b==0xff))
    03078         {
    03079                 bt.Skip(1);
    03080                 inputLen--;
    03081                 bt.Peek(b);
    03082         }
    03083
    03084         reg.CleanNew(RoundupSize(BytesToWords(inputLen)));
    03085
    03086         for (size_t i=inputLen; i > 0; i--)
    03087         {
    03088                 bt.Get(b);
    03089                 reg[(i-1)/WORD_SIZE] |= word(b) << ((i-1)%WORD_SIZE)*8;
    03090         }
    03091
    03092         if (sign == NEGATIVE)
    03093         {
    03094                 for (size_t i=inputLen; i<reg.size()*WORD_SIZE; i++)
    03095                         reg[i/WORD_SIZE] |= word(0xff) << (i%WORD_SIZE)*8;
    03096                 TwosComplement(reg, reg.size());
    03097         }
    03098 }
    03099
    03100 size_t Integer::MinEncodedSize(Signedness signedness) const
    03101 {
    03102         unsigned int outputLen = STDMAX(1U, ByteCount());
    03103         if (signedness == UNSIGNED)
    03104                 return outputLen;
    03105         if (NotNegative() && (GetByte(outputLen-1) & 0x80))
    03106                 outputLen++;
    03107         if (IsNegative() && *this < -Power2(outputLen*8-1))
    03108                 outputLen++;
    03109         return outputLen;
    03110 }
    03111
    03112 void Integer::Encode(byte *output, size_t outputLen, Signedness signedness) const
    03113 {
    03114         ArraySink sink(output, outputLen);
    03115         Encode(sink, outputLen, signedness);
    03116 }
    03117
    03118 void Integer::Encode(BufferedTransformation &bt, size_t outputLen, Signedness signedness) const
    03119 {
    03120         if (signedness == UNSIGNED || NotNegative())
    03121         {
    03122                 for (size_t i=outputLen; i > 0; i--)
    03123                         bt.Put(GetByte(i-1));
    03124         }
    03125         else
    03126         {
    03127                 // take two's complement of *this
    03128                 Integer temp = Integer::Power2(8*STDMAX((size_t)ByteCount(), outputLen)) + *this;
    03129                 temp.Encode(bt, outputLen, UNSIGNED);
    03130         }
    03131 }
    03132
    03133 void Integer::DEREncode(BufferedTransformation &bt) const
    03134 {
    03135         DERGeneralEncoder enc(bt, INTEGER);
    03136         Encode(enc, MinEncodedSize(SIGNED), SIGNED);
    03137         enc.MessageEnd();
    03138 }
    03139
    03140 void Integer::BERDecode(const byte *input, size_t len)
    03141 {
    03142         StringStore store(input, len);
    03143         BERDecode(store);
    03144 }
    03145
    03146 void Integer::BERDecode(BufferedTransformation &bt)
    03147 {
    03148         BERGeneralDecoder dec(bt, INTEGER);
    03149         if (!dec.IsDefiniteLength() || dec.MaxRetrievable() < dec.RemainingLength())
    03150                 BERDecodeError();
    03151         Decode(dec, (size_t)dec.RemainingLength(), SIGNED);
    03152         dec.MessageEnd();
    03153 }
    03154
    03155 void Integer::DEREncodeAsOctetString(BufferedTransformation &bt, size_t length) const
    03156 {
    03157         DERGeneralEncoder enc(bt, OCTET_STRING);
    03158         Encode(enc, length);
    03159         enc.MessageEnd();
    03160 }
    03161
    03162 void Integer::BERDecodeAsOctetString(BufferedTransformation &bt, size_t length)
    03163 {
    03164         BERGeneralDecoder dec(bt, OCTET_STRING);
    03165         if (!dec.IsDefiniteLength() || dec.RemainingLength() != length)
    03166                 BERDecodeError();
    03167         Decode(dec, length);
    03168         dec.MessageEnd();
    03169 }
    03170
    03171 size_t Integer::OpenPGPEncode(byte *output, size_t len) const
    03172 {
    03173         ArraySink sink(output, len);
    03174         return OpenPGPEncode(sink);
    03175 }
    03176
    03177 size_t Integer::OpenPGPEncode(BufferedTransformation &bt) const
    03178 {
    03179         word16 bitCount = BitCount();
    03180         bt.PutWord16(bitCount);
    03181         size_t byteCount = BitsToBytes(bitCount);
    03182         Encode(bt, byteCount);
    03183         return 2 + byteCount;
    03184 }
    03185
    03186 void Integer::OpenPGPDecode(const byte *input, size_t len)
    03187 {
    03188         StringStore store(input, len);
    03189         OpenPGPDecode(store);
    03190 }
    03191
    03192 void Integer::OpenPGPDecode(BufferedTransformation &bt)
    03193 {
    03194         word16 bitCount;
    03195         if (bt.GetWord16(bitCount) != 2 || bt.MaxRetrievable() < BitsToBytes(bitCount))
    03196                 throw OpenPGPDecodeErr();
    03197         Decode(bt, BitsToBytes(bitCount));
    03198 }
    03199
    03200 void Integer::Randomize(RandomNumberGenerator &rng, size_t nbits)
    03201 {
    03202         const size_t nbytes = nbits/8 + 1;
    03203         SecByteBlock buf(nbytes);
    03204         rng.GenerateBlock(buf, nbytes);
    03205         if (nbytes)
    03206                 buf[0] = (byte)Crop(buf[0], nbits % 8);
    03207         Decode(buf, nbytes, UNSIGNED);
    03208 }
    03209
    03210 void Integer::Randomize(RandomNumberGenerator &rng, const Integer &min, const Integer &max)
    03211 {
    03212         if (min > max)
    03213                 throw InvalidArgument("Integer: Min must be no greater than Max");
    03214
    03215         Integer range = max - min;
    03216         const unsigned int nbits = range.BitCount();
    03217
    03218         do
    03219         {
    03220                 Randomize(rng, nbits);
    03221         }
    03222         while (*this > range);
    03223
    03224         *this += min;
    03225 }
    03226
    03227 bool Integer::Randomize(RandomNumberGenerator &rng, const Integer &min, const Integer &max, RandomNumberType rnType, const Integer &equiv, const Integer &mod)
    03228 {
    03229         return GenerateRandomNoThrow(rng, MakeParameters("Min", min)("Max", max)("RandomNumberType", rnType)("EquivalentTo", equiv)("Mod", mod));
    03230 }
    03231
    03232 class KDF2_RNG : public RandomNumberGenerator
    03233 {
    03234 public:
    03235         KDF2_RNG(const byte *seed, size_t seedSize)
    03236                 : m_counter(0), m_counterAndSeed(seedSize + 4)
    03237         {
    03238                 memcpy(m_counterAndSeed + 4, seed, seedSize);
    03239         }
    03240
    03241         void GenerateBlock(byte *output, size_t size)
    03242         {
    03243                 PutWord(false, BIG_ENDIAN_ORDER, m_counterAndSeed, m_counter);
    03244                 ++m_counter;
    03245                 P1363_KDF2<SHA1>::DeriveKey(output, size, m_counterAndSeed, m_counterAndSeed.size(), NULL, 0);
    03246         }
    03247
    03248 private:
    03249         word32 m_counter;
    03250         SecByteBlock m_counterAndSeed;
    03251 };
    03252
    03253 bool Integer::GenerateRandomNoThrow(RandomNumberGenerator &i_rng, const NameValuePairs &params)
    03254 {
    03255         Integer min = params.GetValueWithDefault("Min", Integer::Zero());
    03256         Integer max;
    03257         if (!params.GetValue("Max", max))
    03258         {
    03259                 int bitLength;
    03260                 if (params.GetIntValue("BitLength", bitLength))
    03261                         max = Integer::Power2(bitLength);
    03262                 else
    03263                         throw InvalidArgument("Integer: missing Max argument");
    03264         }
    03265         if (min > max)
    03266                 throw InvalidArgument("Integer: Min must be no greater than Max");
    03267
    03268         Integer equiv = params.GetValueWithDefault("EquivalentTo", Integer::Zero());
    03269         Integer mod = params.GetValueWithDefault("Mod", Integer::One());
    03270
    03271         if (equiv.IsNegative() || equiv >= mod)
    03272                 throw InvalidArgument("Integer: invalid EquivalentTo and/or Mod argument");
    03273
    03274         Integer::RandomNumberType rnType = params.GetValueWithDefault("RandomNumberType", Integer::ANY);
    03275
    03276         member_ptr<KDF2_RNG> kdf2Rng;
    03277         ConstByteArrayParameter seed;
    03278         if (params.GetValue("Seed", seed))
    03279         {
    03280                 ByteQueue bq;
    03281                 DERSequenceEncoder seq(bq);
    03282                 min.DEREncode(seq);
    03283                 max.DEREncode(seq);
    03284                 equiv.DEREncode(seq);
    03285                 mod.DEREncode(seq);
    03286                 DEREncodeUnsigned(seq, rnType);
    03287                 DEREncodeOctetString(seq, seed.begin(), seed.size());
    03288                 seq.MessageEnd();
    03289
    03290                 SecByteBlock finalSeed((size_t)bq.MaxRetrievable());
    03291                 bq.Get(finalSeed, finalSeed.size());
    03292                 kdf2Rng.reset(new KDF2_RNG(finalSeed.begin(), finalSeed.size()));
    03293         }
    03294         RandomNumberGenerator &rng = kdf2Rng.get() ? (RandomNumberGenerator &)*kdf2Rng : i_rng;
    03295
    03296         switch (rnType)
    03297         {
    03298                 case ANY:
    03299                         if (mod == One())
    03300                                 Randomize(rng, min, max);
    03301                         else
    03302                         {
    03303                                 Integer min1 = min + (equiv-min)%mod;
    03304                                 if (max < min1)
    03305                                         return false;
    03306                                 Randomize(rng, Zero(), (max - min1) / mod);
    03307                                 *this *= mod;
    03308                                 *this += min1;
    03309                         }
    03310                         return true;
    03311
    03312                 case PRIME:
    03313                 {
    03314                         const PrimeSelector *pSelector = params.GetValueWithDefault(Name::PointerToPrimeSelector(), (const PrimeSelector *)NULL);
    03315
    03316                         int i;
    03317                         i = 0;
    03318                         while (1)
    03319                         {
    03320                                 if (++i==16)
    03321                                 {
    03322                                         // check if there are any suitable primes in [min, max]
    03323                                         Integer first = min;
    03324                                         if (FirstPrime(first, max, equiv, mod, pSelector))
    03325                                         {
    03326                                                 // if there is only one suitable prime, we're done
    03327                                                 *this = first;
    03328                                                 if (!FirstPrime(first, max, equiv, mod, pSelector))
    03329                                                         return true;
    03330                                         }
    03331                                         else
    03332                                                 return false;
    03333                                 }
    03334
    03335                                 Randomize(rng, min, max);
    03336                                 if (FirstPrime(*this, STDMIN(*this+mod*PrimeSearchInterval(max), max), equiv, mod, pSelector))
    03337                                         return true;
    03338                         }
    03339                 }
    03340
    03341                 default:
    03342                         throw InvalidArgument("Integer: invalid RandomNumberType argument");
    03343         }
    03344 }
    03345
    03346 std::istream& operator>>(std::istream& in, Integer &a)
    03347 {
    03348         char c;
    03349         unsigned int length = 0;
    03350         SecBlock<char> str(length + 16);
    03351
    03352         std::ws(in);
    03353
    03354         do
    03355         {
    03356                 in.read(&c, 1);
    03357                 str[length++] = c;
    03358                 if (length >= str.size())
    03359                         str.Grow(length + 16);
    03360         }
    03361         while (in && (c=='-' || c=='x' || (c>='0' && c<='9') || (c>='a' && c<='f') || (c>='A' && c<='F') || c=='h' || c=='H' || c=='o' || c=='O' || c==',' || c=='.'));
    03362
    03363         if (in.gcount())
    03364                 in.putback(c);
    03365         str[length-1] = '\0';
    03366         a = Integer(str);
    03367
    03368         return in;
    03369 }
    03370
    03371 std::ostream& operator<<(std::ostream& out, const Integer &a)
    03372 {
    03373         // Get relevant conversion specifications from ostream.
    03374         long f = out.flags() & std::ios::basefield; // Get base digits.
    03375         int base, block;
    03376         char suffix;
    03377         switch(f)
    03378         {
    03379         case std::ios::oct :
    03380                 base = 8;
    03381                 block = 8;
    03382                 suffix = 'o';
    03383                 break;
    03384         case std::ios::hex :
    03385                 base = 16;
    03386                 block = 4;
    03387                 suffix = 'h';
    03388                 break;
    03389         default :
    03390                 base = 10;
    03391                 block = 3;
    03392                 suffix = '.';
    03393         }
    03394
    03395         SecBlock<char> s(a.BitCount() / (BitPrecision(base)-1) + 1);
    03396         Integer temp1=a, temp2;
    03397         unsigned i=0;
    03398         const char vec[]="0123456789ABCDEF";
    03399
    03400         if (a.IsNegative())
    03401         {
    03402                 out << '-';
    03403                 temp1.Negate();
    03404         }
    03405
    03406         if (!a)
    03407                 out << '0';
    03408
    03409         while (!!temp1)
    03410         {
    03411                 word digit;
    03412                 Integer::Divide(digit, temp2, temp1, base);
    03413                 s[i++]=vec[digit];
    03414                 temp1=temp2;
    03415         }
    03416
    03417         while (i--)
    03418         {
    03419                 out << s[i];
    03420 //              if (i && !(i%block))
    03421 //                      out << ",";
    03422         }
    03423         return out << suffix;
    03424 }
    03425
    03426 Integer& Integer::operator++()
    03427 {
    03428         if (NotNegative())
    03429         {
    03430                 if (Increment(reg, reg.size()))
    03431                 {
    03432                         reg.CleanGrow(2*reg.size());
    03433                         reg[reg.size()/2]=1;
    03434                 }
    03435         }
    03436         else
    03437         {
    03438                 word borrow = Decrement(reg, reg.size());
    03439                 assert(!borrow);
    03440                 if (WordCount()==0)
    03441                         *this = Zero();
    03442         }
    03443         return *this;
    03444 }
    03445
    03446 Integer& Integer::operator--()
    03447 {
    03448         if (IsNegative())
    03449         {
    03450                 if (Increment(reg, reg.size()))
    03451                 {
    03452                         reg.CleanGrow(2*reg.size());
    03453                         reg[reg.size()/2]=1;
    03454                 }
    03455         }
    03456         else
    03457         {
    03458                 if (Decrement(reg, reg.size()))
    03459                         *this = -One();
    03460         }
    03461         return *this;
    03462 }
    03463
    03464 void PositiveAdd(Integer &sum, const Integer &a, const Integer& b)
    03465 {
    03466         int carry;
    03467         if (a.reg.size() == b.reg.size())
    03468                 carry = Add(sum.reg, a.reg, b.reg, a.reg.size());
    03469         else if (a.reg.size() > b.reg.size())
    03470         {
    03471                 carry = Add(sum.reg, a.reg, b.reg, b.reg.size());
    03472                 CopyWords(sum.reg+b.reg.size(), a.reg+b.reg.size(), a.reg.size()-b.reg.size());
    03473                 carry = Increment(sum.reg+b.reg.size(), a.reg.size()-b.reg.size(), carry);
    03474         }
    03475         else
    03476         {
    03477                 carry = Add(sum.reg, a.reg, b.reg, a.reg.size());
    03478                 CopyWords(sum.reg+a.reg.size(), b.reg+a.reg.size(), b.reg.size()-a.reg.size());
    03479                 carry = Increment(sum.reg+a.reg.size(), b.reg.size()-a.reg.size(), carry);
    03480         }
    03481
    03482         if (carry)
    03483         {
    03484                 sum.reg.CleanGrow(2*sum.reg.size());
    03485                 sum.reg[sum.reg.size()/2] = 1;
    03486         }
    03487         sum.sign = Integer::POSITIVE;
    03488 }
    03489
    03490 void PositiveSubtract(Integer &diff, const Integer &a, const Integer& b)
    03491 {
    03492         unsigned aSize = a.WordCount();
    03493         aSize += aSize%2;
    03494         unsigned bSize = b.WordCount();
    03495         bSize += bSize%2;
    03496
    03497         if (aSize == bSize)
    03498         {
    03499                 if (Compare(a.reg, b.reg, aSize) >= 0)
    03500                 {
    03501                         Subtract(diff.reg, a.reg, b.reg, aSize);
    03502                         diff.sign = Integer::POSITIVE;
    03503                 }
    03504                 else
    03505                 {
    03506                         Subtract(diff.reg, b.reg, a.reg, aSize);
    03507                         diff.sign = Integer::NEGATIVE;
    03508                 }
    03509         }
    03510         else if (aSize > bSize)
    03511         {
    03512                 word borrow = Subtract(diff.reg, a.reg, b.reg, bSize);
    03513                 CopyWords(diff.reg+bSize, a.reg+bSize, aSize-bSize);
    03514                 borrow = Decrement(diff.reg+bSize, aSize-bSize, borrow);
    03515                 assert(!borrow);
    03516                 diff.sign = Integer::POSITIVE;
    03517         }
    03518         else
    03519         {
    03520                 word borrow = Subtract(diff.reg, b.reg, a.reg, aSize);
    03521                 CopyWords(diff.reg+aSize, b.reg+aSize, bSize-aSize);
    03522                 borrow = Decrement(diff.reg+aSize, bSize-aSize, borrow);
    03523                 assert(!borrow);
    03524                 diff.sign = Integer::NEGATIVE;
    03525         }
    03526 }
    03527
    03528 // MSVC .NET 2003 workaround
    03529 template <class T> inline const T& STDMAX2(const T& a, const T& b)
    03530 {
    03531         return a < b ? b : a;
    03532 }
    03533
    03534 Integer Integer::Plus(const Integer& b) const
    03535 {
    03536         Integer sum((word)0, STDMAX2(reg.size(), b.reg.size()));
    03537         if (NotNegative())
    03538         {
    03539                 if (b.NotNegative())
    03540                         PositiveAdd(sum, *this, b);
    03541                 else
    03542                         PositiveSubtract(sum, *this, b);
    03543         }
    03544         else
    03545         {
    03546                 if (b.NotNegative())
    03547                         PositiveSubtract(sum, b, *this);
    03548                 else
    03549                 {
    03550                         PositiveAdd(sum, *this, b);
    03551                         sum.sign = Integer::NEGATIVE;
    03552                 }
    03553         }
    03554         return sum;
    03555 }
    03556
    03557 Integer& Integer::operator+=(const Integer& t)
    03558 {
    03559         reg.CleanGrow(t.reg.size());
    03560         if (NotNegative())
    03561         {
    03562                 if (t.NotNegative())
    03563                         PositiveAdd(*this, *this, t);
    03564                 else
    03565                         PositiveSubtract(*this, *this, t);
    03566         }
    03567         else
    03568         {
    03569                 if (t.NotNegative())
    03570                         PositiveSubtract(*this, t, *this);
    03571                 else
    03572                 {
    03573                         PositiveAdd(*this, *this, t);
    03574                         sign = Integer::NEGATIVE;
    03575                 }
    03576         }
    03577         return *this;
    03578 }
    03579
    03580 Integer Integer::Minus(const Integer& b) const
    03581 {
    03582         Integer diff((word)0, STDMAX2(reg.size(), b.reg.size()));
    03583         if (NotNegative())
    03584         {
    03585                 if (b.NotNegative())
    03586                         PositiveSubtract(diff, *this, b);
    03587                 else
    03588                         PositiveAdd(diff, *this, b);
    03589         }
    03590         else
    03591         {
    03592                 if (b.NotNegative())
    03593                 {
    03594                         PositiveAdd(diff, *this, b);
    03595                         diff.sign = Integer::NEGATIVE;
    03596                 }
    03597                 else
    03598                         PositiveSubtract(diff, b, *this);
    03599         }
    03600         return diff;
    03601 }
    03602
    03603 Integer& Integer::operator-=(const Integer& t)
    03604 {
    03605         reg.CleanGrow(t.reg.size());
    03606         if (NotNegative())
    03607         {
    03608                 if (t.NotNegative())
    03609                         PositiveSubtract(*this, *this, t);
    03610                 else
    03611                         PositiveAdd(*this, *this, t);
    03612         }
    03613         else
    03614         {
    03615                 if (t.NotNegative())
    03616                 {
    03617                         PositiveAdd(*this, *this, t);
    03618                         sign = Integer::NEGATIVE;
    03619                 }
    03620                 else
    03621                         PositiveSubtract(*this, t, *this);
    03622         }
    03623         return *this;
    03624 }
    03625
    03626 Integer& Integer::operator<<=(size_t n)
    03627 {
    03628         const size_t wordCount = WordCount();
    03629         const size_t shiftWords = n / WORD_BITS;
    03630         const unsigned int shiftBits = (unsigned int)(n % WORD_BITS);
    03631
    03632         reg.CleanGrow(RoundupSize(wordCount+BitsToWords(n)));
    03633         ShiftWordsLeftByWords(reg, wordCount + shiftWords, shiftWords);
    03634         ShiftWordsLeftByBits(reg+shiftWords, wordCount+BitsToWords(shiftBits), shiftBits);
    03635         return *this;
    03636 }
    03637
    03638 Integer& Integer::operator>>=(size_t n)
    03639 {
    03640         const size_t wordCount = WordCount();
    03641         const size_t shiftWords = n / WORD_BITS;
    03642         const unsigned int shiftBits = (unsigned int)(n % WORD_BITS);
    03643
    03644         ShiftWordsRightByWords(reg, wordCount, shiftWords);
    03645         if (wordCount > shiftWords)
    03646                 ShiftWordsRightByBits(reg, wordCount-shiftWords, shiftBits);
    03647         if (IsNegative() && WordCount()==0)   // avoid -0
    03648                 *this = Zero();
    03649         return *this;
    03650 }
    03651
    03652 void PositiveMultiply(Integer &product, const Integer &a, const Integer &b)
    03653 {
    03654         size_t aSize = RoundupSize(a.WordCount());
    03655         size_t bSize = RoundupSize(b.WordCount());
    03656
    03657         product.reg.CleanNew(RoundupSize(aSize+bSize));
    03658         product.sign = Integer::POSITIVE;
    03659
    03660         IntegerSecBlock workspace(aSize + bSize);
    03661         AsymmetricMultiply(product.reg, workspace, a.reg, aSize, b.reg, bSize);
    03662 }
    03663
    03664 void Multiply(Integer &product, const Integer &a, const Integer &b)
    03665 {
    03666         PositiveMultiply(product, a, b);
    03667
    03668         if (a.NotNegative() != b.NotNegative())
    03669                 product.Negate();
    03670 }
    03671
    03672 Integer Integer::Times(const Integer &b) const
    03673 {
    03674         Integer product;
    03675         Multiply(product, *this, b);
    03676         return product;
    03677 }
    03678
    03679 /*
    03680 void PositiveDivide(Integer &remainder, Integer &quotient,
    03681                                    const Integer &dividend, const Integer &divisor)
    03682 {
    03683         remainder.reg.CleanNew(divisor.reg.size());
    03684         remainder.sign = Integer::POSITIVE;
    03685         quotient.reg.New(0);
    03686         quotient.sign = Integer::POSITIVE;
    03687         unsigned i=dividend.BitCount();
    03688         while (i--)
    03689         {
    03690                 word overflow = ShiftWordsLeftByBits(remainder.reg, remainder.reg.size(), 1);
    03691                 remainder.reg[0] |= dividend[i];
    03692                 if (overflow || remainder >= divisor)
    03693                 {
    03694                         Subtract(remainder.reg, remainder.reg, divisor.reg, remainder.reg.size());
    03695                         quotient.SetBit(i);
    03696                 }
    03697         }
    03698 }
    03699 */
    03700
    03701 void PositiveDivide(Integer &remainder, Integer &quotient,
    03702                                    const Integer &a, const Integer &b)
    03703 {
    03704         unsigned aSize = a.WordCount();
    03705         unsigned bSize = b.WordCount();
    03706
    03707         if (!bSize)
    03708                 throw Integer::DivideByZero();
    03709
    03710         if (a.PositiveCompare(b) == -1)
    03711         {
    03712                 remainder = a;
    03713                 remainder.sign = Integer::POSITIVE;
    03714                 quotient = Integer::Zero();
    03715                 return;
    03716         }
    03717
    03718         aSize += aSize%2;       // round up to next even number
    03719         bSize += bSize%2;
    03720
    03721         remainder.reg.CleanNew(RoundupSize(bSize));
    03722         remainder.sign = Integer::POSITIVE;
    03723         quotient.reg.CleanNew(RoundupSize(aSize-bSize+2));
    03724         quotient.sign = Integer::POSITIVE;
    03725
    03726         IntegerSecBlock T(aSize+3*(bSize+2));
    03727         Divide(remainder.reg, quotient.reg, T, a.reg, aSize, b.reg, bSize);
    03728 }
    03729
    03730 void Integer::Divide(Integer &remainder, Integer &quotient, const Integer &dividend, const Integer &divisor)
    03731 {
    03732         PositiveDivide(remainder, quotient, dividend, divisor);
    03733
    03734         if (dividend.IsNegative())
    03735         {
    03736                 quotient.Negate();
    03737                 if (remainder.NotZero())
    03738                 {
    03739                         --quotient;
    03740                         remainder = divisor.AbsoluteValue() - remainder;
    03741                 }
    03742         }
    03743
    03744         if (divisor.IsNegative())
    03745                 quotient.Negate();
    03746 }
    03747
    03748 void Integer::DivideByPowerOf2(Integer &r, Integer &q, const Integer &a, unsigned int n)
    03749 {
    03750         q = a;
    03751         q >>= n;
    03752
    03753         const size_t wordCount = BitsToWords(n);
    03754         if (wordCount <= a.WordCount())
    03755         {
    03756                 r.reg.resize(RoundupSize(wordCount));
    03757                 CopyWords(r.reg, a.reg, wordCount);
    03758                 SetWords(r.reg+wordCount, 0, r.reg.size()-wordCount);
    03759                 if (n % WORD_BITS != 0)
    03760                         r.reg[wordCount-1] %= (word(1) << (n % WORD_BITS));
    03761         }
    03762         else
    03763         {
    03764                 r.reg.resize(RoundupSize(a.WordCount()));
    03765                 CopyWords(r.reg, a.reg, r.reg.size());
    03766         }
    03767         r.sign = POSITIVE;
    03768
    03769         if (a.IsNegative() && r.NotZero())
    03770         {
    03771                 --q;
    03772                 r = Power2(n) - r;
    03773         }
    03774 }
    03775
    03776 Integer Integer::DividedBy(const Integer &b) const
    03777 {
    03778         Integer remainder, quotient;
    03779         Integer::Divide(remainder, quotient, *this, b);
    03780         return quotient;
    03781 }
    03782
    03783 Integer Integer::Modulo(const Integer &b) const
    03784 {
    03785         Integer remainder, quotient;
    03786         Integer::Divide(remainder, quotient, *this, b);
    03787         return remainder;
    03788 }
    03789
    03790 void Integer::Divide(word &remainder, Integer &quotient, const Integer &dividend, word divisor)
    03791 {
    03792         if (!divisor)
    03793                 throw Integer::DivideByZero();
    03794
    03795         assert(divisor);
    03796
    03797         if ((divisor & (divisor-1)) == 0)       // divisor is a power of 2
    03798         {
    03799                 quotient = dividend >> (BitPrecision(divisor)-1);
    03800                 remainder = dividend.reg[0] & (divisor-1);
    03801                 return;
    03802         }
    03803
    03804         unsigned int i = dividend.WordCount();
    03805         quotient.reg.CleanNew(RoundupSize(i));
    03806         remainder = 0;
    03807         while (i--)
    03808         {
    03809                 quotient.reg[i] = DWord(dividend.reg[i], remainder) / divisor;
    03810                 remainder = DWord(dividend.reg[i], remainder) % divisor;
    03811         }
    03812
    03813         if (dividend.NotNegative())
    03814                 quotient.sign = POSITIVE;
    03815         else
    03816         {
    03817                 quotient.sign = NEGATIVE;
    03818                 if (remainder)
    03819                 {
    03820                         --quotient;
    03821                         remainder = divisor - remainder;
    03822                 }
    03823         }
    03824 }
    03825
    03826 Integer Integer::DividedBy(word b) const
    03827 {
    03828         word remainder;
    03829         Integer quotient;
    03830         Integer::Divide(remainder, quotient, *this, b);
    03831         return quotient;
    03832 }
    03833
    03834 word Integer::Modulo(word divisor) const
    03835 {
    03836         if (!divisor)
    03837                 throw Integer::DivideByZero();
    03838
    03839         assert(divisor);
    03840
    03841         word remainder;
    03842
    03843         if ((divisor & (divisor-1)) == 0)       // divisor is a power of 2
    03844                 remainder = reg[0] & (divisor-1);
    03845         else
    03846         {
    03847                 unsigned int i = WordCount();
    03848
    03849                 if (divisor <= 5)
    03850                 {
    03851                         DWord sum(0, 0);
    03852                         while (i--)
    03853                                 sum += reg[i];
    03854                         remainder = sum % divisor;
    03855                 }
    03856                 else
    03857                 {
    03858                         remainder = 0;
    03859                         while (i--)
    03860                                 remainder = DWord(reg[i], remainder) % divisor;
    03861                 }
    03862         }
    03863
    03864         if (IsNegative() && remainder)
    03865                 remainder = divisor - remainder;
    03866
    03867         return remainder;
    03868 }
    03869
    03870 void Integer::Negate()
    03871 {
    03872         if (!!(*this))  // don't flip sign if *this==0
    03873                 sign = Sign(1-sign);
    03874 }
    03875
    03876 int Integer::PositiveCompare(const Integer& t) const
    03877 {
    03878         unsigned size = WordCount(), tSize = t.WordCount();
    03879
    03880         if (size == tSize)
    03881                 return CryptoPP::Compare(reg, t.reg, size);
    03882         else
    03883                 return size > tSize ? 1 : -1;
    03884 }
    03885
    03886 int Integer::Compare(const Integer& t) const
    03887 {
    03888         if (NotNegative())
    03889         {
    03890                 if (t.NotNegative())
    03891                         return PositiveCompare(t);
    03892                 else
    03893                         return 1;
    03894         }
    03895         else
    03896         {
    03897                 if (t.NotNegative())
    03898                         return -1;
    03899                 else
    03900                         return -PositiveCompare(t);
    03901         }
    03902 }
    03903
    03904 Integer Integer::SquareRoot() const
    03905 {
    03906         if (!IsPositive())
    03907                 return Zero();
    03908
    03909         // overestimate square root
    03910         Integer x, y = Power2((BitCount()+1)/2);
    03911         assert(y*y >= *this);
    03912
    03913         do
    03914         {
    03915                 x = y;
    03916                 y = (x + *this/x) >> 1;
    03917         } while (y<x);
    03918
    03919         return x;
    03920 }
    03921
    03922 bool Integer::IsSquare() const
    03923 {
    03924         Integer r = SquareRoot();
    03925         return *this == r.Squared();
    03926 }
    03927
    03928 bool Integer::IsUnit() const
    03929 {
    03930         return (WordCount() == 1) && (reg[0] == 1);
    03931 }
    03932
    03933 Integer Integer::MultiplicativeInverse() const
    03934 {
    03935         return IsUnit() ? *this : Zero();
    03936 }
    03937
    03938 Integer a_times_b_mod_c(const Integer &x, const Integer& y, const Integer& m)
    03939 {
    03940         return x*y%m;
    03941 }
    03942
    03943 Integer a_exp_b_mod_c(const Integer &x, const Integer& e, const Integer& m)
    03944 {
    03945         ModularArithmetic mr(m);
    03946         return mr.Exponentiate(x, e);
    03947 }
    03948
    03949 Integer Integer::Gcd(const Integer &a, const Integer &b)
    03950 {
    03951         return EuclideanDomainOf<Integer>().Gcd(a, b);
    03952 }
    03953
    03954 Integer Integer::InverseMod(const Integer &m) const
    03955 {
    03956         assert(m.NotNegative());
    03957
    03958         if (IsNegative() || *this>=m)
    03959                 return (*this%m).InverseMod(m);
    03960
    03961         if (m.IsEven())
    03962         {
    03963                 if (!m || IsEven())
    03964                         return Zero();  // no inverse
    03965                 if (*this == One())
    03966                         return One();
    03967
    03968                 Integer u = m.InverseMod(*this);
    03969                 return !u ? Zero() : (m*(*this-u)+1)/(*this);
    03970         }
    03971
    03972         SecBlock<word> T(m.reg.size() * 4);
    03973         Integer r((word)0, m.reg.size());
    03974         unsigned k = AlmostInverse(r.reg, T, reg, reg.size(), m.reg, m.reg.size());
    03975         DivideByPower2Mod(r.reg, r.reg, k, m.reg, m.reg.size());
    03976         return r;
    03977 }
    03978
    03979 word Integer::InverseMod(word mod) const
    03980 {
    03981         word g0 = mod, g1 = *this % mod;
    03982         word v0 = 0, v1 = 1;
    03983         word y;
    03984
    03985         while (g1)
    03986         {
    03987                 if (g1 == 1)
    03988                         return v1;
    03989                 y = g0 / g1;
    03990                 g0 = g0 % g1;
    03991                 v0 += y * v1;
    03992
    03993                 if (!g0)
    03994                         break;
    03995                 if (g0 == 1)
    03996                         return mod-v0;
    03997                 y = g1 / g0;
    03998                 g1 = g1 % g0;
    03999                 v1 += y * v0;
    04000         }
    04001         return 0;
    04002 }
    04003
    04004 // ********************************************************
    04005
    04006 ModularArithmetic::ModularArithmetic(BufferedTransformation &bt)
    04007 {
    04008         BERSequenceDecoder seq(bt);
    04009         OID oid(seq);
    04010         if (oid != ASN1::prime_field())
    04011                 BERDecodeError();
    04012         m_modulus.BERDecode(seq);
    04013         seq.MessageEnd();
    04014         m_result.reg.resize(m_modulus.reg.size());
    04015 }
    04016
    04017 void ModularArithmetic::DEREncode(BufferedTransformation &bt) const
    04018 {
    04019         DERSequenceEncoder seq(bt);
    04020         ASN1::prime_field().DEREncode(seq);
    04021         m_modulus.DEREncode(seq);
    04022         seq.MessageEnd();
    04023 }
    04024
    04025 void ModularArithmetic::DEREncodeElement(BufferedTransformation &out, const Element &a) const
    04026 {
    04027         a.DEREncodeAsOctetString(out, MaxElementByteLength());
    04028 }
    04029
    04030 void ModularArithmetic::BERDecodeElement(BufferedTransformation &in, Element &a) const
    04031 {
    04032         a.BERDecodeAsOctetString(in, MaxElementByteLength());
    04033 }
    04034
    04035 const Integer& ModularArithmetic::Half(const Integer &a) const
    04036 {
    04037         if (a.reg.size()==m_modulus.reg.size())
    04038         {
    04039                 CryptoPP::DivideByPower2Mod(m_result.reg.begin(), a.reg, 1, m_modulus.reg, a.reg.size());
    04040                 return m_result;
    04041         }
    04042         else
    04043                 return m_result1 = (a.IsEven() ? (a >> 1) : ((a+m_modulus) >> 1));
    04044 }
    04045
    04046 const Integer& ModularArithmetic::Add(const Integer &a, const Integer &b) const
    04047 {
    04048         if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
    04049         {
    04050                 if (CryptoPP::Add(m_result.reg.begin(), a.reg, b.reg, a.reg.size())
    04051                         || Compare(m_result.reg, m_modulus.reg, a.reg.size()) >= 0)
    04052                 {
    04053                         CryptoPP::Subtract(m_result.reg.begin(), m_result.reg, m_modulus.reg, a.reg.size());
    04054                 }
    04055                 return m_result;
    04056         }
    04057         else
    04058         {
    04059                 m_result1 = a+b;
    04060                 if (m_result1 >= m_modulus)
    04061                         m_result1 -= m_modulus;
    04062                 return m_result1;
    04063         }
    04064 }
    04065
    04066 Integer& ModularArithmetic::Accumulate(Integer &a, const Integer &b) const
    04067 {
    04068         if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
    04069         {
    04070                 if (CryptoPP::Add(a.reg, a.reg, b.reg, a.reg.size())
    04071                         || Compare(a.reg, m_modulus.reg, a.reg.size()) >= 0)
    04072                 {
    04073                         CryptoPP::Subtract(a.reg, a.reg, m_modulus.reg, a.reg.size());
    04074                 }
    04075         }
    04076         else
    04077         {
    04078                 a+=b;
    04079                 if (a>=m_modulus)
    04080                         a-=m_modulus;
    04081         }
    04082
    04083         return a;
    04084 }
    04085
    04086 const Integer& ModularArithmetic::Subtract(const Integer &a, const Integer &b) const
    04087 {
    04088         if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
    04089         {
    04090                 if (CryptoPP::Subtract(m_result.reg.begin(), a.reg, b.reg, a.reg.size()))
    04091                         CryptoPP::Add(m_result.reg.begin(), m_result.reg, m_modulus.reg, a.reg.size());
    04092                 return m_result;
    04093         }
    04094         else
    04095         {
    04096                 m_result1 = a-b;
    04097                 if (m_result1.IsNegative())
    04098                         m_result1 += m_modulus;
    04099                 return m_result1;
    04100         }
    04101 }
    04102
    04103 Integer& ModularArithmetic::Reduce(Integer &a, const Integer &b) const
    04104 {
    04105         if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
    04106         {
    04107                 if (CryptoPP::Subtract(a.reg, a.reg, b.reg, a.reg.size()))
    04108                         CryptoPP::Add(a.reg, a.reg, m_modulus.reg, a.reg.size());
    04109         }
    04110         else
    04111         {
    04112                 a-=b;
    04113                 if (a.IsNegative())
    04114                         a+=m_modulus;
    04115         }
    04116
    04117         return a;
    04118 }
    04119
    04120 const Integer& ModularArithmetic::Inverse(const Integer &a) const
    04121 {
    04122         if (!a)
    04123                 return a;
    04124
    04125         CopyWords(m_result.reg.begin(), m_modulus.reg, m_modulus.reg.size());
    04126         if (CryptoPP::Subtract(m_result.reg.begin(), m_result.reg, a.reg, a.reg.size()))
    04127                 Decrement(m_result.reg.begin()+a.reg.size(), m_modulus.reg.size()-a.reg.size());
    04128
    04129         return m_result;
    04130 }
    04131
    04132 Integer ModularArithmetic::CascadeExponentiate(const Integer &x, const Integer &e1, const Integer &y, const Integer &e2) const
    04133 {
    04134         if (m_modulus.IsOdd())
    04135         {
    04136                 MontgomeryRepresentation dr(m_modulus);
    04137                 return dr.ConvertOut(dr.CascadeExponentiate(dr.ConvertIn(x), e1, dr.ConvertIn(y), e2));
    04138         }
    04139         else
    04140                 return AbstractRing<Integer>::CascadeExponentiate(x, e1, y, e2);
    04141 }
    04142
    04143 void ModularArithmetic::SimultaneousExponentiate(Integer *results, const Integer &base, const Integer *exponents, unsigned int exponentsCount) const
    04144 {
    04145         if (m_modulus.IsOdd())
    04146         {
    04147                 MontgomeryRepresentation dr(m_modulus);
    04148                 dr.SimultaneousExponentiate(results, dr.ConvertIn(base), exponents, exponentsCount);
    04149                 for (unsigned int i=0; i<exponentsCount; i++)
    04150                         results[i] = dr.ConvertOut(results[i]);
    04151         }
    04152         else
    04153                 AbstractRing<Integer>::SimultaneousExponentiate(results, base, exponents, exponentsCount);
    04154 }
    04155
    04156 MontgomeryRepresentation::MontgomeryRepresentation(const Integer &m)    // modulus must be odd
    04157         : ModularArithmetic(m),
    04158           m_u((word)0, m_modulus.reg.size()),
    04159           m_workspace(5*m_modulus.reg.size())
    04160 {
    04161         if (!m_modulus.IsOdd())
    04162                 throw InvalidArgument("MontgomeryRepresentation: Montgomery representation requires an odd modulus");
    04163
    04164         RecursiveInverseModPower2(m_u.reg, m_workspace, m_modulus.reg, m_modulus.reg.size());
    04165 }
    04166
    04167 const Integer& MontgomeryRepresentation::Multiply(const Integer &a, const Integer &b) const
    04168 {
    04169         word *const T = m_workspace.begin();
    04170         word *const R = m_result.reg.begin();
    04171         const size_t N = m_modulus.reg.size();
    04172         assert(a.reg.size()<=N && b.reg.size()<=N);
    04173
    04174         AsymmetricMultiply(T, T+2*N, a.reg, a.reg.size(), b.reg, b.reg.size());
    04175         SetWords(T+a.reg.size()+b.reg.size(), 0, 2*N-a.reg.size()-b.reg.size());
    04176         MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
    04177         return m_result;
    04178 }
    04179
    04180 const Integer& MontgomeryRepresentation::Square(const Integer &a) const
    04181 {
    04182         word *const T = m_workspace.begin();
    04183         word *const R = m_result.reg.begin();
    04184         const size_t N = m_modulus.reg.size();
    04185         assert(a.reg.size()<=N);
    04186
    04187         CryptoPP::Square(T, T+2*N, a.reg, a.reg.size());
    04188         SetWords(T+2*a.reg.size(), 0, 2*N-2*a.reg.size());
    04189         MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
    04190         return m_result;
    04191 }
    04192
    04193 Integer MontgomeryRepresentation::ConvertOut(const Integer &a) const
    04194 {
    04195         word *const T = m_workspace.begin();
    04196         word *const R = m_result.reg.begin();
    04197         const size_t N = m_modulus.reg.size();
    04198         assert(a.reg.size()<=N);
    04199
    04200         CopyWords(T, a.reg, a.reg.size());
    04201         SetWords(T+a.reg.size(), 0, 2*N-a.reg.size());
    04202         MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
    04203         return m_result;
    04204 }
    04205
    04206 const Integer& MontgomeryRepresentation::MultiplicativeInverse(const Integer &a) const
    04207 {
    04208 //        return (EuclideanMultiplicativeInverse(a, modulus)<<(2*WORD_BITS*modulus.reg.size()))%modulus;
    04209         word *const T = m_workspace.begin();
    04210         word *const R = m_result.reg.begin();
    04211         const size_t N = m_modulus.reg.size();
    04212         assert(a.reg.size()<=N);
    04213
    04214         CopyWords(T, a.reg, a.reg.size());
    04215         SetWords(T+a.reg.size(), 0, 2*N-a.reg.size());
    04216         MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
    04217         unsigned k = AlmostInverse(R, T, R, N, m_modulus.reg, N);
    04218
    04219 //      cout << "k=" << k << " N*32=" << 32*N << endl;
    04220
    04221         if (k>N*WORD_BITS)
    04222                 DivideByPower2Mod(R, R, k-N*WORD_BITS, m_modulus.reg, N);
    04223         else
    04224                 MultiplyByPower2Mod(R, R, N*WORD_BITS-k, m_modulus.reg, N);
    04225
    04226         return m_result;
    04227 }
    04228
    04229 NAMESPACE_END
    04230
    04231 #endif