~zooko/cryptopp/trunk

« back to all changes in this revision

Viewing changes to integer.cpp

  • Committer: weidai
  • Date: 2003-07-29 01:18:33 UTC
  • Revision ID: svn-v4:57ff6487-cd31-0410-9ec3-f628ee90f5f0:trunk/c5:118
fix potential threading problem with initialization of static objects

Show diffs side-by-side

added added

removed removed

Lines of Context:
17
17
 
18
18
#include <iostream>
19
19
 
20
 
#ifdef _M_X64
21
 
#include <Intrin.h>
22
 
#endif
23
 
 
24
20
#ifdef SSE2_INTRINSICS_AVAILABLE
25
 
        #ifdef __GNUC__
26
 
                #include <xmmintrin.h>
27
 
                #include <signal.h>
28
 
                #include <setjmp.h>
29
 
                #ifdef CRYPTOPP_MEMALIGN_AVAILABLE
30
 
                        #include <malloc.h>
31
 
                #else
32
 
                        #include <stdlib.h>
33
 
                #endif
34
 
        #else
35
 
                #include <emmintrin.h>
36
 
        #endif
 
21
#include <emmintrin.h>
37
22
#elif defined(_MSC_VER) && defined(_M_IX86)
38
 
        #pragma message("You do not seem to have the Visual C++ Processor Pack installed, so use of SSE2 intrinsics will be disabled.")
39
 
#elif defined(__GNUC__) && defined(__i386__)
40
 
        #warning "You do not have GCC 3.3 or later, or did not specify -msse2 compiler option, so use of SSE2 intrinsics will be disabled."
 
23
#pragma message("You do no seem to have the Visual C++ Processor Pack installed, so use of SSE2 intrinsics will be disabled.")
41
24
#endif
42
25
 
43
26
NAMESPACE_BEGIN(CryptoPP)
44
27
 
45
 
bool AssignIntToInteger(const std::type_info &valueType, void *pInteger, const void *pInt)
 
28
bool FunctionAssignIntToInteger(const std::type_info &valueType, void *pInteger, const void *pInt)
46
29
{
47
30
        if (valueType != typeid(Integer))
48
31
                return false;
50
33
        return true;
51
34
}
52
35
 
53
 
#ifdef SSE2_INTRINSICS_AVAILABLE
 
36
static const char s_RunAtStartup = (AssignIntToInteger = FunctionAssignIntToInteger, 0);
 
37
 
 
38
#if defined(SSE2_INTRINSICS_AVAILABLE) || defined(_MSC_VER)
54
39
template <class T>
55
40
CPP_TYPENAME AllocatorBase<T>::pointer AlignedAllocator<T>::allocate(size_type n, const void *)
56
41
{
57
 
        CheckSize(n);
58
 
        if (n == 0)
59
 
                return NULL;
 
42
#ifdef SSE2_INTRINSICS_AVAILABLE
60
43
        if (n >= 4)
61
 
        {
62
 
                void *p;
63
 
        #ifdef CRYPTOPP_MM_MALLOC_AVAILABLE
64
 
                while (!(p = _mm_malloc(sizeof(T)*n, 16)))
65
 
        #elif defined(CRYPTOPP_MEMALIGN_AVAILABLE)
66
 
                while (!(p = memalign(16, sizeof(T)*n)))
67
 
        #elif defined(CRYPTOPP_MALLOC_ALIGNMENT_IS_16)
68
 
                while (!(p = malloc(sizeof(T)*n)))
69
 
        #else
70
 
                while (!(p = (byte *)malloc(sizeof(T)*n + 8)))  // assume malloc alignment is at least 8
71
 
        #endif
72
 
                        CallNewHandler();
73
 
 
74
 
        #ifdef CRYPTOPP_NO_ALIGNED_ALLOC
75
 
                assert(m_pBlock == NULL);
76
 
                m_pBlock = p;
77
 
                if (!IsAlignedOn(p, 16))
78
 
                {
79
 
                        assert(IsAlignedOn(p, 8));
80
 
                        p = (byte *)p + 8;
81
 
                }
82
 
        #endif
83
 
 
84
 
                assert(IsAlignedOn(p, 16));
85
 
                return (T*)p;
86
 
        }
87
 
        return new T[n];
 
44
                return (T *)_mm_malloc(sizeof(T)*n, 16);
 
45
        else
 
46
#endif
 
47
                return new T[n];
88
48
}
89
49
 
90
50
template <class T>
91
51
void AlignedAllocator<T>::deallocate(void *p, size_type n)
92
52
{
93
53
        memset(p, 0, n*sizeof(T));
 
54
#ifdef SSE2_INTRINSICS_AVAILABLE
94
55
        if (n >= 4)
95
 
        {
96
 
                #ifdef CRYPTOPP_MM_MALLOC_AVAILABLE
97
 
                        _mm_free(p);
98
 
                #elif defined(CRYPTOPP_NO_ALIGNED_ALLOC)
99
 
                        assert(m_pBlock == p || (byte *)m_pBlock+8 == p);
100
 
                        free(m_pBlock);
101
 
                        m_pBlock = NULL;
102
 
                #else
103
 
                        free(p);
104
 
                #endif
105
 
        }
 
56
                _mm_free(p);
106
57
        else
107
 
                delete [] (T *)p;
 
58
#endif
 
59
                delete [] p;
108
60
}
109
 
 
110
 
template class CRYPTOPP_DLL AlignedAllocator<word>;
111
61
#endif
112
62
 
113
 
static int Compare(const word *A, const word *B, size_t N)
 
63
static int Compare(const word *A, const word *B, unsigned int N)
114
64
{
115
65
        while (N--)
116
66
                if (A[N] > B[N])
121
71
        return 0;
122
72
}
123
73
 
124
 
static int Increment(word *A, size_t N, word B=1)
 
74
static word Increment(word *A, unsigned int N, word B=1)
125
75
{
126
76
        assert(N);
127
77
        word t = A[0];
134
84
        return 1;
135
85
}
136
86
 
137
 
static int Decrement(word *A, size_t N, word B=1)
 
87
static word Decrement(word *A, unsigned int N, word B=1)
138
88
{
139
89
        assert(N);
140
90
        word t = A[0];
147
97
        return 1;
148
98
}
149
99
 
150
 
static void TwosComplement(word *A, size_t N)
 
100
static void TwosComplement(word *A, unsigned int N)
151
101
{
152
102
        Decrement(A, N);
153
103
        for (unsigned i=0; i<N; i++)
208
158
                        __asm__("mulq %3" : "=d" (r.m_halfs.high), "=a" (r.m_halfs.low) : "a" (a), "rm" (b) : "cc");
209
159
                #elif defined(__mips64)
210
160
                        __asm__("dmultu %2,%3" : "=h" (r.m_halfs.high), "=l" (r.m_halfs.low) : "r" (a), "r" (b));
211
 
                #elif defined(_M_X64)
212
 
                        r.m_halfs.low = _umul128(a, b, &r.m_halfs.high);
213
161
                #elif defined(_M_IX86)
214
162
                        // for testing
215
163
                        word64 t = (word64)a * b;
460
408
class Portable
461
409
{
462
410
public:
463
 
        static int Add(word *C, const word *A, const word *B, size_t N);
464
 
        static int Subtract(word *C, const word *A, const word *B, size_t N);
 
411
        static word Add(word *C, const word *A, const word *B, unsigned int N);
 
412
        static word Subtract(word *C, const word *A, const word *B, unsigned int N);
465
413
 
466
414
        static inline void Multiply2(word *C, const word *A, const word *B);
467
415
        static inline word Multiply2Add(word *C, const word *A, const word *B);
480
428
        static inline unsigned int SquareRecursionLimit() {return 4;}
481
429
};
482
430
 
483
 
int Portable::Add(word *C, const word *A, const word *B, size_t N)
 
431
word Portable::Add(word *C, const word *A, const word *B, unsigned int N)
484
432
{
485
433
        assert (N%2 == 0);
486
434
 
492
440
                u = DWord(A[i+1]) + B[i+1] + u.GetHighHalf();
493
441
                C[i+1] = u.GetLowHalf();
494
442
        }
495
 
        return int(u.GetHighHalf());
 
443
        return u.GetHighHalf();
496
444
}
497
445
 
498
 
int Portable::Subtract(word *C, const word *A, const word *B, size_t N)
 
446
word Portable::Subtract(word *C, const word *A, const word *B, unsigned int N)
499
447
{
500
448
        assert (N%2 == 0);
501
449
 
507
455
                u = (DWord) A[i+1] - B[i+1] - u.GetHighHalfAsBorrow();
508
456
                C[i+1] = u.GetLowHalf();
509
457
        }
510
 
        return int(0-u.GetHighHalf());
 
458
        return 0-u.GetHighHalf();
511
459
}
512
460
 
513
461
void Portable::Multiply2(word *C, const word *A, const word *B)
692
640
 
693
641
void Portable::Square4(word *R, const word *A)
694
642
{
695
 
#ifdef _MSC_VER
696
 
        // VC60 workaround: MSVC 6.0 has an optimization bug that makes
697
 
        // (dword)A*B where either A or B has been cast to a dword before
698
 
        // very expensive. Revisit this function when this
699
 
        // bug is fixed.
700
 
        Multiply4(R, A, A);
701
 
#else
702
643
        const word *B = A;
703
644
        DWord p, q;
704
645
        word c, d, e;
725
666
        p = DWord::MultiplyAndAdd(A[3], A[3], d);
726
667
        R[6] = p.GetLowHalf();
727
668
        R[7] = e + p.GetHighHalf();
728
 
#endif
729
669
}
730
670
 
731
671
void Portable::Multiply8(word *R, const word *A, const word *B)
893
833
#undef SquAcc
894
834
#undef SaveSquAcc
895
835
 
896
 
#ifdef CRYPTOPP_X86ASM_AVAILABLE
897
 
 
898
 
// ************** x86 feature detection ***************
899
 
 
900
 
static bool s_sse2Enabled = true;
901
 
 
902
 
static void CpuId(word32 input, word32 *output)
903
 
{
904
 
#ifdef __GNUC__
905
 
        __asm__
906
 
        (
907
 
                // save ebx in case -fPIC is being used
908
 
                "push %%ebx; cpuid; mov %%ebx, %%edi; pop %%ebx"
909
 
                : "=a" (output[0]), "=D" (output[1]), "=c" (output[2]), "=d" (output[3])
910
 
                : "a" (input)
911
 
        );
912
 
#else
913
 
        __asm
914
 
        {
915
 
                mov eax, input
 
836
// CodeWarrior defines _MSC_VER
 
837
#if defined(_MSC_VER) && !defined(__MWERKS__) && defined(_M_IX86) && (_M_IX86<=700)
 
838
 
 
839
class PentiumOptimized : public Portable
 
840
{
 
841
public:
 
842
        static word __fastcall Add(word *C, const word *A, const word *B, unsigned int N);
 
843
        static word __fastcall Subtract(word *C, const word *A, const word *B, unsigned int N);
 
844
// TODO test this with .NET #if _MSC_VER < 1300
 
845
        static inline void Square4(word *R, const word *A)
 
846
        {
 
847
                // VC60 workaround: MSVC 6.0 has an optimization bug that makes
 
848
                // (dword)A*B where either A or B has been cast to a dword before
 
849
                // very expensive. Revisit this function when this
 
850
                // bug is fixed.
 
851
                Multiply4(R, A, A);
 
852
        }
 
853
//#endif
 
854
};
 
855
 
 
856
typedef PentiumOptimized LowLevel;
 
857
 
 
858
__declspec(naked) word __fastcall PentiumOptimized::Add(word *C, const word *A, const word *B, unsigned int N)
 
859
{
 
860
        __asm
 
861
        {
 
862
                push ebp
 
863
                push ebx
 
864
                push esi
 
865
                push edi
 
866
 
 
867
                mov esi, [esp+24]       ; N
 
868
                mov ebx, [esp+20]       ; B
 
869
 
 
870
                // now: ebx = B, ecx = C, edx = A, esi = N
 
871
 
 
872
                sub ecx, edx    // hold the distance between C & A so we can add this to A to get C
 
873
                xor eax, eax    // clear eax
 
874
 
 
875
                sub eax, esi    // eax is a negative index from end of B
 
876
                lea ebx, [ebx+4*esi]    // ebx is end of B
 
877
 
 
878
                sar eax, 1              // unit of eax is now dwords; this also clears the carry flag
 
879
                jz      loopend         // if no dwords then nothing to do
 
880
 
 
881
loopstart:
 
882
                mov    esi,[edx]                        // load lower word of A
 
883
                mov    ebp,[edx+4]                      // load higher word of A
 
884
 
 
885
                mov    edi,[ebx+8*eax]          // load lower word of B
 
886
                lea    edx,[edx+8]                      // advance A and C
 
887
 
 
888
                adc    esi,edi                          // add lower words
 
889
                mov    edi,[ebx+8*eax+4]        // load higher word of B
 
890
 
 
891
                adc    ebp,edi                          // add higher words
 
892
                inc    eax                                      // advance B
 
893
 
 
894
                mov    [edx+ecx-8],esi          // store lower word result
 
895
                mov    [edx+ecx-4],ebp          // store higher word result
 
896
 
 
897
                jnz    loopstart                        // loop until eax overflows and becomes zero
 
898
 
 
899
loopend:
 
900
                adc eax, 0              // store carry into eax (return result register)
 
901
                pop edi
 
902
                pop esi
 
903
                pop ebx
 
904
                pop ebp
 
905
                ret 8
 
906
        }
 
907
}
 
908
 
 
909
__declspec(naked) word __fastcall PentiumOptimized::Subtract(word *C, const word *A, const word *B, unsigned int N)
 
910
{
 
911
        __asm
 
912
        {
 
913
                push ebp
 
914
                push ebx
 
915
                push esi
 
916
                push edi
 
917
 
 
918
                mov esi, [esp+24]       ; N
 
919
                mov ebx, [esp+20]       ; B
 
920
 
 
921
                sub ecx, edx
 
922
                xor eax, eax
 
923
 
 
924
                sub eax, esi
 
925
                lea ebx, [ebx+4*esi]
 
926
 
 
927
                sar eax, 1
 
928
                jz      loopend
 
929
 
 
930
loopstart:
 
931
                mov    esi,[edx]
 
932
                mov    ebp,[edx+4]
 
933
 
 
934
                mov    edi,[ebx+8*eax]
 
935
                lea    edx,[edx+8]
 
936
 
 
937
                sbb    esi,edi
 
938
                mov    edi,[ebx+8*eax+4]
 
939
 
 
940
                sbb    ebp,edi
 
941
                inc    eax
 
942
 
 
943
                mov    [edx+ecx-8],esi
 
944
                mov    [edx+ecx-4],ebp
 
945
 
 
946
                jnz    loopstart
 
947
 
 
948
loopend:
 
949
                adc eax, 0
 
950
                pop edi
 
951
                pop esi
 
952
                pop ebx
 
953
                pop ebp
 
954
                ret 8
 
955
        }
 
956
}
 
957
 
 
958
#ifdef SSE2_INTRINSICS_AVAILABLE
 
959
 
 
960
static bool GetSSE2Capability()
 
961
{
 
962
        word32 b;
 
963
 
 
964
        __asm
 
965
        {
 
966
                mov             eax, 1
916
967
                cpuid
917
 
                mov edi, output
918
 
                mov [edi], eax
919
 
                mov [edi+4], ebx
920
 
                mov [edi+8], ecx
921
 
                mov [edi+12], edx
922
 
        }
923
 
#endif
924
 
}
925
 
 
926
 
#ifdef SSE2_INTRINSICS_AVAILABLE
927
 
#ifndef _MSC_VER
928
 
static jmp_buf s_env;
929
 
static void SigIllHandler(int)
930
 
{
931
 
        longjmp(s_env, 1);
932
 
}
933
 
#endif
934
 
 
935
 
static bool HasSSE2()
936
 
{
937
 
        if (!s_sse2Enabled)
938
 
                return false;
939
 
 
940
 
        word32 cpuid[4];
941
 
        CpuId(1, cpuid);
942
 
        if ((cpuid[3] & (1 << 26)) == 0)
943
 
                return false;
944
 
 
945
 
#ifdef _MSC_VER
946
 
    __try
947
 
        {
948
 
        __asm xorpd xmm0, xmm0        // executing SSE2 instruction
949
 
        }
950
 
    __except (1)
951
 
        {
952
 
                return false;
953
 
    }
954
 
        return true;
955
 
#else
956
 
        typedef void (*SigHandler)(int);
957
 
 
958
 
        SigHandler oldHandler = signal(SIGILL, SigIllHandler);
959
 
        if (oldHandler == SIG_ERR)
960
 
                return false;
961
 
 
962
 
        bool result = true;
963
 
        if (setjmp(s_env))
964
 
                result = false;
965
 
        else
966
 
                __asm __volatile ("xorps %xmm0, %xmm0");
967
 
 
968
 
        signal(SIGILL, oldHandler);
969
 
        return result;
970
 
#endif
971
 
}
972
 
#endif
973
 
 
974
 
static bool IsP4()
975
 
{
976
 
        word32 cpuid[4];
977
 
 
978
 
        CpuId(0, cpuid);
979
 
        std::swap(cpuid[2], cpuid[3]);
980
 
        if (memcmp(cpuid+1, "GenuineIntel", 12) != 0)
981
 
                return false;
982
 
 
983
 
        CpuId(1, cpuid);
984
 
        return ((cpuid[0] >> 8) & 0xf) == 0xf;
985
 
}
986
 
 
987
 
// ************** Pentium/P4 optimizations ***************
988
 
 
989
 
class PentiumOptimized : public Portable
990
 
{
991
 
public:
992
 
        static int Add(word *C, const word *A, const word *B, size_t N);
993
 
        static int Subtract(word *C, const word *A, const word *B, size_t N);
994
 
        static void Multiply4(word *C, const word *A, const word *B);
995
 
        static void Multiply8(word *C, const word *A, const word *B);
996
 
        static void Multiply8Bottom(word *C, const word *A, const word *B);
997
 
};
998
 
 
999
 
class P4Optimized
1000
 
{
1001
 
public:
1002
 
        static int Add(word *C, const word *A, const word *B, size_t N);
1003
 
        static int Subtract(word *C, const word *A, const word *B, size_t N);
1004
 
#ifdef SSE2_INTRINSICS_AVAILABLE
1005
 
        static void Multiply4(word *C, const word *A, const word *B);
1006
 
        static void Multiply8(word *C, const word *A, const word *B);
1007
 
        static void Multiply8Bottom(word *C, const word *A, const word *B);
1008
 
#endif
1009
 
};
1010
 
 
1011
 
typedef int (* PAddSub)(word *C, const word *A, const word *B, size_t N);
1012
 
typedef void (* PMul)(word *C, const word *A, const word *B);
1013
 
 
1014
 
static PAddSub s_pAdd, s_pSub;
1015
 
#ifdef SSE2_INTRINSICS_AVAILABLE
1016
 
static PMul s_pMul4, s_pMul8, s_pMul8B;
1017
 
#endif
1018
 
 
1019
 
static void SetPentiumFunctionPointers()
1020
 
{
1021
 
        if (IsP4())
1022
 
        {
1023
 
                s_pAdd = &P4Optimized::Add;
1024
 
                s_pSub = &P4Optimized::Subtract;
1025
 
        }
1026
 
        else
1027
 
        {
1028
 
                s_pAdd = &PentiumOptimized::Add;
1029
 
                s_pSub = &PentiumOptimized::Subtract;
1030
 
        }
1031
 
 
1032
 
#ifdef SSE2_INTRINSICS_AVAILABLE
1033
 
        if (HasSSE2())
1034
 
        {
1035
 
                s_pMul4 = &P4Optimized::Multiply4;
1036
 
                s_pMul8 = &P4Optimized::Multiply8;
1037
 
                s_pMul8B = &P4Optimized::Multiply8Bottom;
1038
 
        }
1039
 
        else
1040
 
        {
1041
 
                s_pMul4 = &PentiumOptimized::Multiply4;
1042
 
                s_pMul8 = &PentiumOptimized::Multiply8;
1043
 
                s_pMul8B = &PentiumOptimized::Multiply8Bottom;
1044
 
        }
1045
 
#endif
1046
 
}
 
968
                mov             b, edx
 
969
        }
 
970
 
 
971
        return (b & (1 << 26)) != 0;
 
972
}
 
973
 
 
974
bool g_sse2DetectionDone = false, g_sse2Detected, g_sse2Enabled = true;
1047
975
 
1048
976
void DisableSSE2()
1049
977
{
1050
 
        s_sse2Enabled = false;
1051
 
        SetPentiumFunctionPointers();
1052
 
}
1053
 
 
1054
 
class LowLevel : public PentiumOptimized
 
978
        g_sse2Enabled = false;
 
979
}
 
980
 
 
981
static inline bool HasSSE2()
 
982
{
 
983
        if (g_sse2Enabled && !g_sse2DetectionDone)
 
984
        {
 
985
                g_sse2Detected = GetSSE2Capability();
 
986
                g_sse2DetectionDone = true;
 
987
        }
 
988
        return g_sse2Enabled && g_sse2Detected;
 
989
}
 
990
 
 
991
class P4Optimized : public PentiumOptimized
1055
992
{
1056
993
public:
1057
 
        inline static int Add(word *C, const word *A, const word *B, size_t N)
1058
 
                {return s_pAdd(C, A, B, N);}
1059
 
        inline static int Subtract(word *C, const word *A, const word *B, size_t N)
1060
 
                {return s_pSub(C, A, B, N);}
1061
 
        inline static void Square4(word *R, const word *A)
1062
 
                {Multiply4(R, A, A);}
1063
 
#ifdef SSE2_INTRINSICS_AVAILABLE
1064
 
        inline static void Multiply4(word *C, const word *A, const word *B)
1065
 
                {s_pMul4(C, A, B);}
1066
 
        inline static void Multiply8(word *C, const word *A, const word *B)
1067
 
                {s_pMul8(C, A, B);}
1068
 
        inline static void Multiply8Bottom(word *C, const word *A, const word *B)
1069
 
                {s_pMul8B(C, A, B);}
1070
 
#endif
 
994
        static word __fastcall Add(word *C, const word *A, const word *B, unsigned int N);
 
995
        static word __fastcall Subtract(word *C, const word *A, const word *B, unsigned int N);
 
996
        static void Multiply4(word *C, const word *A, const word *B);
 
997
        static void Multiply8(word *C, const word *A, const word *B);
 
998
        static inline void Square4(word *R, const word *A)
 
999
        {
 
1000
                Multiply4(R, A, A);
 
1001
        }
 
1002
        static void Multiply8Bottom(word *C, const word *A, const word *B);
1071
1003
};
1072
1004
 
1073
 
// use some tricks to share assembly code between MSVC and GCC
1074
 
#ifdef _MSC_VER
1075
 
        #define CRYPTOPP_NAKED __declspec(naked)
1076
 
        #define AS1(x) __asm x
1077
 
        #define AS2(x, y) __asm x, y
1078
 
        #define AddPrologue \
1079
 
                __asm   push ebp \
1080
 
                __asm   push ebx \
1081
 
                __asm   push esi \
1082
 
                __asm   push edi \
1083
 
                __asm   mov             ecx, [esp+20] \
1084
 
                __asm   mov             edx, [esp+24] \
1085
 
                __asm   mov             ebx, [esp+28] \
1086
 
                __asm   mov             esi, [esp+32]
1087
 
        #define AddEpilogue \
1088
 
                __asm   pop edi \
1089
 
                __asm   pop esi \
1090
 
                __asm   pop ebx \
1091
 
                __asm   pop ebp \
1092
 
                __asm   ret
1093
 
        #define MulPrologue \
1094
 
                __asm   push ebp \
1095
 
                __asm   push ebx \
1096
 
                __asm   push esi \
1097
 
                __asm   push edi \
1098
 
                __asm   mov ecx, [esp+28] \
1099
 
                __asm   mov esi, [esp+24] \
1100
 
                __asm   push [esp+20]
1101
 
        #define MulEpilogue \
1102
 
                __asm   add esp, 4 \
1103
 
                __asm   pop edi \
1104
 
                __asm   pop esi \
1105
 
                __asm   pop ebx \
1106
 
                __asm   pop ebp \
1107
 
                __asm   ret
1108
 
#else
1109
 
        #define CRYPTOPP_NAKED
1110
 
        #define AS1(x) #x ";"
1111
 
        #define AS2(x, y) #x ", " #y ";"
1112
 
        #define AddPrologue \
1113
 
                __asm__ __volatile__ \
1114
 
                ( \
1115
 
                        "push %%ebx;"   /* save this manually, in case of -fPIC */ \
1116
 
                        "mov %2, %%ebx;" \
1117
 
                        ".intel_syntax noprefix;" \
1118
 
                        "push ebp;"
1119
 
        #define AddEpilogue \
1120
 
                        "pop ebp;" \
1121
 
                        ".att_syntax prefix;" \
1122
 
                        "pop %%ebx;" \
1123
 
                                        : \
1124
 
                                        : "c" (C), "d" (A), "m" (B), "S" (N) \
1125
 
                                        : "%edi", "memory", "cc" \
1126
 
                );
1127
 
        #define MulPrologue \
1128
 
                __asm__ __volatile__ \
1129
 
                ( \
1130
 
                        "push %%ebx;"   /* save this manually, in case of -fPIC */ \
1131
 
                        "push %%ebp;" \
1132
 
                        "push %0;" \
1133
 
                        ".intel_syntax noprefix;"
1134
 
        #define MulEpilogue \
1135
 
                        "add esp, 4;" \
1136
 
                        "pop ebp;" \
1137
 
                        "pop ebx;" \
1138
 
                        ".att_syntax prefix;" \
1139
 
                        : \
1140
 
                        : "rm" (Z), "S" (X), "c" (Y) \
1141
 
                        : "%eax", "%edx", "%edi", "memory", "cc" \
1142
 
                );
1143
 
#endif
1144
 
 
1145
 
CRYPTOPP_NAKED int PentiumOptimized::Add(word *C, const word *A, const word *B, size_t N)
1146
 
{
1147
 
        AddPrologue
1148
 
 
1149
 
        // now: ebx = B, ecx = C, edx = A, esi = N
1150
 
        AS2(    sub ecx, edx)   // hold the distance between C & A so we can add this to A to get C
1151
 
        AS2(    xor eax, eax)   // clear eax
1152
 
 
1153
 
        AS2(    sub eax, esi)   // eax is a negative index from end of B
1154
 
        AS2(    lea ebx, [ebx+4*esi])   // ebx is end of B
1155
 
 
1156
 
        AS2(    sar eax, 1)             // unit of eax is now dwords; this also clears the carry flag
1157
 
        AS1(    jz      loopendAdd)             // if no dwords then nothing to do
1158
 
 
1159
 
        AS1(loopstartAdd:)
1160
 
        AS2(    mov    esi,[edx])                       // load lower word of A
1161
 
        AS2(    mov    ebp,[edx+4])                     // load higher word of A
1162
 
 
1163
 
        AS2(    mov    edi,[ebx+8*eax])         // load lower word of B
1164
 
        AS2(    lea    edx,[edx+8])                     // advance A and C
1165
 
 
1166
 
        AS2(    adc    esi,edi)                         // add lower words
1167
 
        AS2(    mov    edi,[ebx+8*eax+4])       // load higher word of B
1168
 
 
1169
 
        AS2(    adc    ebp,edi)                         // add higher words
1170
 
        AS1(    inc    eax)                                     // advance B
1171
 
 
1172
 
        AS2(    mov    [edx+ecx-8],esi)         // store lower word result
1173
 
        AS2(    mov    [edx+ecx-4],ebp)         // store higher word result
1174
 
 
1175
 
        AS1(    jnz    loopstartAdd)                    // loop until eax overflows and becomes zero
1176
 
 
1177
 
        AS1(loopendAdd:)
1178
 
        AS2(    adc eax, 0)             // store carry into eax (return result register)
1179
 
 
1180
 
        AddEpilogue
1181
 
}
1182
 
 
1183
 
CRYPTOPP_NAKED int PentiumOptimized::Subtract(word *C, const word *A, const word *B, size_t N)
1184
 
{
1185
 
        AddPrologue
1186
 
 
1187
 
        // now: ebx = B, ecx = C, edx = A, esi = N
1188
 
        AS2(    sub ecx, edx)   // hold the distance between C & A so we can add this to A to get C
1189
 
        AS2(    xor eax, eax)   // clear eax
1190
 
 
1191
 
        AS2(    sub eax, esi)   // eax is a negative index from end of B
1192
 
        AS2(    lea ebx, [ebx+4*esi])   // ebx is end of B
1193
 
 
1194
 
        AS2(    sar eax, 1)             // unit of eax is now dwords; this also clears the carry flag
1195
 
        AS1(    jz      loopendSub)             // if no dwords then nothing to do
1196
 
 
1197
 
        AS1(loopstartSub:)
1198
 
        AS2(    mov    esi,[edx])                       // load lower word of A
1199
 
        AS2(    mov    ebp,[edx+4])                     // load higher word of A
1200
 
 
1201
 
        AS2(    mov    edi,[ebx+8*eax])         // load lower word of B
1202
 
        AS2(    lea    edx,[edx+8])                     // advance A and C
1203
 
 
1204
 
        AS2(    sbb    esi,edi)                         // subtract lower words
1205
 
        AS2(    mov    edi,[ebx+8*eax+4])       // load higher word of B
1206
 
 
1207
 
        AS2(    sbb    ebp,edi)                         // subtract higher words
1208
 
        AS1(    inc    eax)                                     // advance B
1209
 
 
1210
 
        AS2(    mov    [edx+ecx-8],esi)         // store lower word result
1211
 
        AS2(    mov    [edx+ecx-4],ebp)         // store higher word result
1212
 
 
1213
 
        AS1(    jnz    loopstartSub)                    // loop until eax overflows and becomes zero
1214
 
 
1215
 
        AS1(loopendSub:)
1216
 
        AS2(    adc eax, 0)             // store carry into eax (return result register)
1217
 
 
1218
 
        AddEpilogue
1219
 
}
1220
 
 
1221
 
// On Pentium 4, the adc and sbb instructions are very expensive, so avoid them.
1222
 
 
1223
 
CRYPTOPP_NAKED int P4Optimized::Add(word *C, const word *A, const word *B, size_t N)
1224
 
{
1225
 
        AddPrologue
1226
 
 
1227
 
        // now: ebx = B, ecx = C, edx = A, esi = N
1228
 
        AS2(    xor             eax, eax)
1229
 
        AS1(    neg             esi)
1230
 
        AS1(    jz              loopendAddP4)           // if no dwords then nothing to do
1231
 
 
1232
 
        AS2(    mov             edi, [edx])
1233
 
        AS2(    mov             ebp, [ebx])
1234
 
        AS1(    jmp             carry1AddP4)
1235
 
 
1236
 
        AS1(loopstartAddP4:)
1237
 
        AS2(    mov             edi, [edx+8])
1238
 
        AS2(    add             ecx, 8)
1239
 
        AS2(    add             edx, 8)
1240
 
        AS2(    mov             ebp, [ebx])
1241
 
        AS2(    add             edi, eax)
1242
 
        AS1(    jc              carry1AddP4)
1243
 
        AS2(    xor             eax, eax)
1244
 
 
1245
 
        AS1(carry1AddP4:)
1246
 
        AS2(    add             edi, ebp)
1247
 
        AS2(    mov             ebp, 1)
1248
 
        AS2(    mov             [ecx], edi)
1249
 
        AS2(    mov             edi, [edx+4])
1250
 
        AS2(    cmovc   eax, ebp)
1251
 
        AS2(    mov             ebp, [ebx+4])
1252
 
        AS2(    add             ebx, 8)
1253
 
        AS2(    add             edi, eax)
1254
 
        AS1(    jc              carry2AddP4)
1255
 
        AS2(    xor             eax, eax)
1256
 
 
1257
 
        AS1(carry2AddP4:)
1258
 
        AS2(    add             edi, ebp)
1259
 
        AS2(    mov             ebp, 1)
1260
 
        AS2(    cmovc   eax, ebp)
1261
 
        AS2(    mov             [ecx+4], edi)
1262
 
        AS2(    add             esi, 2)
1263
 
        AS1(    jnz             loopstartAddP4)
1264
 
 
1265
 
        AS1(loopendAddP4:)
1266
 
 
1267
 
        AddEpilogue
1268
 
}
1269
 
 
1270
 
CRYPTOPP_NAKED int P4Optimized::Subtract(word *C, const word *A, const word *B, size_t N)
1271
 
{
1272
 
        AddPrologue
1273
 
 
1274
 
        // now: ebx = B, ecx = C, edx = A, esi = N
1275
 
        AS2(    xor             eax, eax)
1276
 
        AS1(    neg             esi)
1277
 
        AS1(    jz              loopendSubP4)           // if no dwords then nothing to do
1278
 
 
1279
 
        AS2(    mov             edi, [edx])
1280
 
        AS2(    mov             ebp, [ebx])
1281
 
        AS1(    jmp             carry1SubP4)
1282
 
 
1283
 
        AS1(loopstartSubP4:)
1284
 
        AS2(    mov             edi, [edx+8])
1285
 
        AS2(    add             edx, 8)
1286
 
        AS2(    add             ecx, 8)
1287
 
        AS2(    mov             ebp, [ebx])
1288
 
        AS2(    sub             edi, eax)
1289
 
        AS1(    jc              carry1SubP4)
1290
 
        AS2(    xor             eax, eax)
1291
 
 
1292
 
        AS1(carry1SubP4:)
1293
 
        AS2(    sub             edi, ebp)
1294
 
        AS2(    mov             ebp, 1)
1295
 
        AS2(    mov             [ecx], edi)
1296
 
        AS2(    mov             edi, [edx+4])
1297
 
        AS2(    cmovc   eax, ebp)
1298
 
        AS2(    mov             ebp, [ebx+4])
1299
 
        AS2(    add             ebx, 8)
1300
 
        AS2(    sub             edi, eax)
1301
 
        AS1(    jc              carry2SubP4)
1302
 
        AS2(    xor             eax, eax)
1303
 
 
1304
 
        AS1(carry2SubP4:)
1305
 
        AS2(    sub             edi, ebp)
1306
 
        AS2(    mov             ebp, 1)
1307
 
        AS2(    cmovc   eax, ebp)
1308
 
        AS2(    mov             [ecx+4], edi)
1309
 
        AS2(    add             esi, 2)
1310
 
        AS1(    jnz             loopstartSubP4)
1311
 
 
1312
 
        AS1(loopendSubP4:)
1313
 
 
1314
 
        AddEpilogue
1315
 
}
1316
 
 
1317
 
// multiply assembly code originally contributed by Leonard Janke
1318
 
 
1319
 
#define MulStartup \
1320
 
        AS2(xor ebp, ebp) \
1321
 
        AS2(xor edi, edi) \
1322
 
        AS2(xor ebx, ebx) 
1323
 
 
1324
 
#define MulShiftCarry \
1325
 
        AS2(mov ebp, edx) \
1326
 
        AS2(mov edi, ebx) \
1327
 
        AS2(xor ebx, ebx)
1328
 
 
1329
 
#define MulAccumulateBottom(i,j) \
1330
 
        AS2(mov eax, [ecx+4*j]) \
1331
 
        AS2(imul eax, dword ptr [esi+4*i]) \
1332
 
        AS2(add ebp, eax)
1333
 
 
1334
 
#define MulAccumulate(i,j) \
1335
 
        AS2(mov eax, [ecx+4*j]) \
1336
 
        AS1(mul dword ptr [esi+4*i]) \
1337
 
        AS2(add ebp, eax) \
1338
 
        AS2(adc edi, edx) \
1339
 
        AS2(adc bl, bh)
1340
 
 
1341
 
#define MulStoreDigit(i)  \
1342
 
        AS2(mov edx, edi) \
1343
 
        AS2(mov edi, [esp]) \
1344
 
        AS2(mov [edi+4*i], ebp)
1345
 
 
1346
 
#define MulLastDiagonal(digits) \
1347
 
        AS2(mov eax, [ecx+4*(digits-1)]) \
1348
 
        AS1(mul dword ptr [esi+4*(digits-1)]) \
1349
 
        AS2(add ebp, eax) \
1350
 
        AS2(adc edx, edi) \
1351
 
        AS2(mov edi, [esp]) \
1352
 
        AS2(mov [edi+4*(2*digits-2)], ebp) \
1353
 
        AS2(mov [edi+4*(2*digits-1)], edx)
1354
 
 
1355
 
CRYPTOPP_NAKED void PentiumOptimized::Multiply4(word* Z, const word* X, const word* Y)
1356
 
{
1357
 
        MulPrologue
1358
 
        // now: [esp] = Z, esi = X, ecx = Y
1359
 
        MulStartup
1360
 
        MulAccumulate(0,0)
1361
 
        MulStoreDigit(0)
1362
 
        MulShiftCarry
1363
 
 
1364
 
        MulAccumulate(1,0)
1365
 
        MulAccumulate(0,1)
1366
 
        MulStoreDigit(1)
1367
 
        MulShiftCarry
1368
 
 
1369
 
        MulAccumulate(2,0)
1370
 
        MulAccumulate(1,1)
1371
 
        MulAccumulate(0,2)
1372
 
        MulStoreDigit(2)
1373
 
        MulShiftCarry
1374
 
 
1375
 
        MulAccumulate(3,0)
1376
 
        MulAccumulate(2,1)
1377
 
        MulAccumulate(1,2)
1378
 
        MulAccumulate(0,3)
1379
 
        MulStoreDigit(3)
1380
 
        MulShiftCarry
1381
 
 
1382
 
        MulAccumulate(3,1)
1383
 
        MulAccumulate(2,2)
1384
 
        MulAccumulate(1,3)
1385
 
        MulStoreDigit(4)
1386
 
        MulShiftCarry
1387
 
 
1388
 
        MulAccumulate(3,2)
1389
 
        MulAccumulate(2,3)
1390
 
        MulStoreDigit(5)
1391
 
        MulShiftCarry
1392
 
 
1393
 
        MulLastDiagonal(4)
1394
 
        MulEpilogue
1395
 
}
1396
 
 
1397
 
CRYPTOPP_NAKED void PentiumOptimized::Multiply8(word* Z, const word* X, const word* Y)
1398
 
{
1399
 
        MulPrologue
1400
 
        // now: [esp] = Z, esi = X, ecx = Y
1401
 
        MulStartup
1402
 
        MulAccumulate(0,0)
1403
 
        MulStoreDigit(0)
1404
 
        MulShiftCarry
1405
 
 
1406
 
        MulAccumulate(1,0)
1407
 
        MulAccumulate(0,1)
1408
 
        MulStoreDigit(1)
1409
 
        MulShiftCarry
1410
 
 
1411
 
        MulAccumulate(2,0)
1412
 
        MulAccumulate(1,1)
1413
 
        MulAccumulate(0,2)
1414
 
        MulStoreDigit(2)
1415
 
        MulShiftCarry
1416
 
 
1417
 
        MulAccumulate(3,0)
1418
 
        MulAccumulate(2,1)
1419
 
        MulAccumulate(1,2)
1420
 
        MulAccumulate(0,3)
1421
 
        MulStoreDigit(3)
1422
 
        MulShiftCarry
1423
 
 
1424
 
        MulAccumulate(4,0)
1425
 
        MulAccumulate(3,1)
1426
 
        MulAccumulate(2,2)
1427
 
        MulAccumulate(1,3)
1428
 
        MulAccumulate(0,4)
1429
 
        MulStoreDigit(4)
1430
 
        MulShiftCarry
1431
 
 
1432
 
        MulAccumulate(5,0)
1433
 
        MulAccumulate(4,1)
1434
 
        MulAccumulate(3,2)
1435
 
        MulAccumulate(2,3)
1436
 
        MulAccumulate(1,4)
1437
 
        MulAccumulate(0,5)
1438
 
        MulStoreDigit(5)
1439
 
        MulShiftCarry
1440
 
 
1441
 
        MulAccumulate(6,0)
1442
 
        MulAccumulate(5,1)
1443
 
        MulAccumulate(4,2)
1444
 
        MulAccumulate(3,3)
1445
 
        MulAccumulate(2,4)
1446
 
        MulAccumulate(1,5)
1447
 
        MulAccumulate(0,6)
1448
 
        MulStoreDigit(6)
1449
 
        MulShiftCarry
1450
 
 
1451
 
        MulAccumulate(7,0)
1452
 
        MulAccumulate(6,1)
1453
 
        MulAccumulate(5,2)
1454
 
        MulAccumulate(4,3)
1455
 
        MulAccumulate(3,4)
1456
 
        MulAccumulate(2,5)
1457
 
        MulAccumulate(1,6)
1458
 
        MulAccumulate(0,7)
1459
 
        MulStoreDigit(7)
1460
 
        MulShiftCarry
1461
 
 
1462
 
        MulAccumulate(7,1)
1463
 
        MulAccumulate(6,2)
1464
 
        MulAccumulate(5,3)
1465
 
        MulAccumulate(4,4)
1466
 
        MulAccumulate(3,5)
1467
 
        MulAccumulate(2,6)
1468
 
        MulAccumulate(1,7)
1469
 
        MulStoreDigit(8)
1470
 
        MulShiftCarry
1471
 
 
1472
 
        MulAccumulate(7,2)
1473
 
        MulAccumulate(6,3)
1474
 
        MulAccumulate(5,4)
1475
 
        MulAccumulate(4,5)
1476
 
        MulAccumulate(3,6)
1477
 
        MulAccumulate(2,7)
1478
 
        MulStoreDigit(9)
1479
 
        MulShiftCarry
1480
 
 
1481
 
        MulAccumulate(7,3)
1482
 
        MulAccumulate(6,4)
1483
 
        MulAccumulate(5,5)
1484
 
        MulAccumulate(4,6)
1485
 
        MulAccumulate(3,7)
1486
 
        MulStoreDigit(10)
1487
 
        MulShiftCarry
1488
 
 
1489
 
        MulAccumulate(7,4)
1490
 
        MulAccumulate(6,5)
1491
 
        MulAccumulate(5,6)
1492
 
        MulAccumulate(4,7)
1493
 
        MulStoreDigit(11)
1494
 
        MulShiftCarry
1495
 
 
1496
 
        MulAccumulate(7,5)
1497
 
        MulAccumulate(6,6)
1498
 
        MulAccumulate(5,7)
1499
 
        MulStoreDigit(12)
1500
 
        MulShiftCarry
1501
 
 
1502
 
        MulAccumulate(7,6)
1503
 
        MulAccumulate(6,7)
1504
 
        MulStoreDigit(13)
1505
 
        MulShiftCarry
1506
 
 
1507
 
        MulLastDiagonal(8)
1508
 
        MulEpilogue
1509
 
}
1510
 
 
1511
 
CRYPTOPP_NAKED void PentiumOptimized::Multiply8Bottom(word* Z, const word* X, const word* Y)
1512
 
{
1513
 
        MulPrologue
1514
 
        // now: [esp] = Z, esi = X, ecx = Y
1515
 
        MulStartup
1516
 
        MulAccumulate(0,0)
1517
 
        MulStoreDigit(0)
1518
 
        MulShiftCarry
1519
 
 
1520
 
        MulAccumulate(1,0)
1521
 
        MulAccumulate(0,1)
1522
 
        MulStoreDigit(1)
1523
 
        MulShiftCarry
1524
 
 
1525
 
        MulAccumulate(2,0)
1526
 
        MulAccumulate(1,1)
1527
 
        MulAccumulate(0,2)
1528
 
        MulStoreDigit(2)
1529
 
        MulShiftCarry
1530
 
 
1531
 
        MulAccumulate(3,0)
1532
 
        MulAccumulate(2,1)
1533
 
        MulAccumulate(1,2)
1534
 
        MulAccumulate(0,3)
1535
 
        MulStoreDigit(3)
1536
 
        MulShiftCarry
1537
 
 
1538
 
        MulAccumulate(4,0)
1539
 
        MulAccumulate(3,1)
1540
 
        MulAccumulate(2,2)
1541
 
        MulAccumulate(1,3)
1542
 
        MulAccumulate(0,4)
1543
 
        MulStoreDigit(4)
1544
 
        MulShiftCarry
1545
 
 
1546
 
        MulAccumulate(5,0)
1547
 
        MulAccumulate(4,1)
1548
 
        MulAccumulate(3,2)
1549
 
        MulAccumulate(2,3)
1550
 
        MulAccumulate(1,4)
1551
 
        MulAccumulate(0,5)
1552
 
        MulStoreDigit(5)
1553
 
        MulShiftCarry
1554
 
 
1555
 
        MulAccumulate(6,0)
1556
 
        MulAccumulate(5,1)
1557
 
        MulAccumulate(4,2)
1558
 
        MulAccumulate(3,3)
1559
 
        MulAccumulate(2,4)
1560
 
        MulAccumulate(1,5)
1561
 
        MulAccumulate(0,6)
1562
 
        MulStoreDigit(6)
1563
 
        MulShiftCarry
1564
 
 
1565
 
        MulAccumulateBottom(7,0)
1566
 
        MulAccumulateBottom(6,1)
1567
 
        MulAccumulateBottom(5,2)
1568
 
        MulAccumulateBottom(4,3)
1569
 
        MulAccumulateBottom(3,4)
1570
 
        MulAccumulateBottom(2,5)
1571
 
        MulAccumulateBottom(1,6)
1572
 
        MulAccumulateBottom(0,7)
1573
 
        MulStoreDigit(7)
1574
 
        MulEpilogue
1575
 
}
1576
 
 
1577
 
#undef AS1
1578
 
#undef AS2
1579
 
 
1580
 
#else   // not x86 - no processor specific code at this layer
1581
 
 
1582
 
typedef Portable LowLevel;
1583
 
 
1584
 
#endif
1585
 
 
1586
 
#ifdef SSE2_INTRINSICS_AVAILABLE
1587
 
 
1588
 
#ifdef __GNUC__
1589
 
#define CRYPTOPP_FASTCALL
1590
 
#else
1591
 
#define CRYPTOPP_FASTCALL __fastcall
1592
 
#endif
1593
 
 
1594
 
static void CRYPTOPP_FASTCALL P4_Mul(__m128i *C, const __m128i *A, const __m128i *B)
 
1005
static void __fastcall P4_Mul(__m128i *C, const __m128i *A, const __m128i *B)
1595
1006
{
1596
1007
        __m128i a3210 = _mm_load_si128(A);
1597
1008
        __m128i b3210 = _mm_load_si128(B);
1661
1072
 
1662
1073
        __m64 s1, s2;
1663
1074
 
1664
 
        __m64 w1 = _mm_cvtsi32_si64(w[1]);
 
1075
        __m64 w1 = _m_from_int(w[1]);
1665
1076
        __m64 w4 = mw[2];
1666
1077
        __m64 w6 = mw[3];
1667
1078
        __m64 w8 = mw[4];
1672
1083
        __m64 w18 = mw[9];
1673
1084
        __m64 w20 = mw[10];
1674
1085
        __m64 w22 = mw[11];
1675
 
        __m64 w26 = _mm_cvtsi32_si64(w[26]);
 
1086
        __m64 w26 = _m_from_int(w[26]);
1676
1087
 
1677
1088
        s1 = _mm_add_si64(w1, w4);
1678
 
        C[1] = _mm_cvtsi64_si32(s1);
1679
 
        s1 = _mm_srli_si64(s1, 32);
 
1089
        C[1] = _m_to_int(s1);
 
1090
        s1 = _m_psrlqi(s1, 32);
1680
1091
 
1681
1092
        s2 = _mm_add_si64(w6, w8);
1682
1093
        s1 = _mm_add_si64(s1, s2);
1683
 
        C[2] = _mm_cvtsi64_si32(s1);
1684
 
        s1 = _mm_srli_si64(s1, 32);
 
1094
        C[2] = _m_to_int(s1);
 
1095
        s1 = _m_psrlqi(s1, 32);
1685
1096
 
1686
1097
        s2 = _mm_add_si64(w10, w12);
1687
1098
        s1 = _mm_add_si64(s1, s2);
1688
 
        C[3] = _mm_cvtsi64_si32(s1);
1689
 
        s1 = _mm_srli_si64(s1, 32);
 
1099
        C[3] = _m_to_int(s1);
 
1100
        s1 = _m_psrlqi(s1, 32);
1690
1101
 
1691
1102
        s2 = _mm_add_si64(w14, w16);
1692
1103
        s1 = _mm_add_si64(s1, s2);
1693
 
        C[4] = _mm_cvtsi64_si32(s1);
1694
 
        s1 = _mm_srli_si64(s1, 32);
 
1104
        C[4] = _m_to_int(s1);
 
1105
        s1 = _m_psrlqi(s1, 32);
1695
1106
 
1696
1107
        s2 = _mm_add_si64(w18, w20);
1697
1108
        s1 = _mm_add_si64(s1, s2);
1698
 
        C[5] = _mm_cvtsi64_si32(s1);
1699
 
        s1 = _mm_srli_si64(s1, 32);
 
1109
        C[5] = _m_to_int(s1);
 
1110
        s1 = _m_psrlqi(s1, 32);
1700
1111
 
1701
1112
        s2 = _mm_add_si64(w22, w26);
1702
1113
        s1 = _mm_add_si64(s1, s2);
1703
 
        C[6] = _mm_cvtsi64_si32(s1);
1704
 
        s1 = _mm_srli_si64(s1, 32);
 
1114
        C[6] = _m_to_int(s1);
 
1115
        s1 = _m_psrlqi(s1, 32);
1705
1116
 
1706
 
        C[7] = _mm_cvtsi64_si32(s1) + w[27];
 
1117
        C[7] = _m_to_int(s1) + w[27];
1707
1118
        _mm_empty();
1708
1119
}
1709
1120
 
1731
1142
 
1732
1143
        __m64 s1, s2, s3, s4;
1733
1144
 
1734
 
        __m64 w1 = _mm_cvtsi32_si64(w[1]);
 
1145
        __m64 w1 = _m_from_int(w[1]);
1735
1146
        __m64 w4 = mw[2];
1736
1147
        __m64 w6 = mw[3];
1737
1148
        __m64 w8 = mw[4];
1742
1153
        __m64 w18 = mw[9];
1743
1154
        __m64 w20 = mw[10];
1744
1155
        __m64 w22 = mw[11];
1745
 
        __m64 w26 = _mm_cvtsi32_si64(w[26]);
1746
 
        __m64 w27 = _mm_cvtsi32_si64(w[27]);
 
1156
        __m64 w26 = _m_from_int(w[26]);
 
1157
        __m64 w27 = _m_from_int(w[27]);
1747
1158
 
1748
 
        __m64 x0 = _mm_cvtsi32_si64(x[0]);
1749
 
        __m64 x1 = _mm_cvtsi32_si64(x[1]);
 
1159
        __m64 x0 = _m_from_int(x[0]);
 
1160
        __m64 x1 = _m_from_int(x[1]);
1750
1161
        __m64 x4 = mx[2];
1751
1162
        __m64 x6 = mx[3];
1752
1163
        __m64 x8 = mx[4];
1757
1168
        __m64 x18 = mx[9];
1758
1169
        __m64 x20 = mx[10];
1759
1170
        __m64 x22 = mx[11];
1760
 
        __m64 x26 = _mm_cvtsi32_si64(x[26]);
1761
 
        __m64 x27 = _mm_cvtsi32_si64(x[27]);
 
1171
        __m64 x26 = _m_from_int(x[26]);
 
1172
        __m64 x27 = _m_from_int(x[27]);
1762
1173
 
1763
 
        __m64 y0 = _mm_cvtsi32_si64(y[0]);
1764
 
        __m64 y1 = _mm_cvtsi32_si64(y[1]);
 
1174
        __m64 y0 = _m_from_int(y[0]);
 
1175
        __m64 y1 = _m_from_int(y[1]);
1765
1176
        __m64 y4 = my[2];
1766
1177
        __m64 y6 = my[3];
1767
1178
        __m64 y8 = my[4];
1772
1183
        __m64 y18 = my[9];
1773
1184
        __m64 y20 = my[10];
1774
1185
        __m64 y22 = my[11];
1775
 
        __m64 y26 = _mm_cvtsi32_si64(y[26]);
1776
 
        __m64 y27 = _mm_cvtsi32_si64(y[27]);
 
1186
        __m64 y26 = _m_from_int(y[26]);
 
1187
        __m64 y27 = _m_from_int(y[27]);
1777
1188
 
1778
 
        __m64 z0 = _mm_cvtsi32_si64(z[0]);
1779
 
        __m64 z1 = _mm_cvtsi32_si64(z[1]);
 
1189
        __m64 z0 = _m_from_int(z[0]);
 
1190
        __m64 z1 = _m_from_int(z[1]);
1780
1191
        __m64 z4 = mz[2];
1781
1192
        __m64 z6 = mz[3];
1782
1193
        __m64 z8 = mz[4];
1787
1198
        __m64 z18 = mz[9];
1788
1199
        __m64 z20 = mz[10];
1789
1200
        __m64 z22 = mz[11];
1790
 
        __m64 z26 = _mm_cvtsi32_si64(z[26]);
 
1201
        __m64 z26 = _m_from_int(z[26]);
1791
1202
 
1792
1203
        s1 = _mm_add_si64(w1, w4);
1793
 
        C[1] = _mm_cvtsi64_si32(s1);
1794
 
        s1 = _mm_srli_si64(s1, 32);
 
1204
        C[1] = _m_to_int(s1);
 
1205
        s1 = _m_psrlqi(s1, 32);
1795
1206
 
1796
1207
        s2 = _mm_add_si64(w6, w8);
1797
1208
        s1 = _mm_add_si64(s1, s2);
1798
 
        C[2] = _mm_cvtsi64_si32(s1);
1799
 
        s1 = _mm_srli_si64(s1, 32);
 
1209
        C[2] = _m_to_int(s1);
 
1210
        s1 = _m_psrlqi(s1, 32);
1800
1211
 
1801
1212
        s2 = _mm_add_si64(w10, w12);
1802
1213
        s1 = _mm_add_si64(s1, s2);
1803
 
        C[3] = _mm_cvtsi64_si32(s1);
1804
 
        s1 = _mm_srli_si64(s1, 32);
 
1214
        C[3] = _m_to_int(s1);
 
1215
        s1 = _m_psrlqi(s1, 32);
1805
1216
 
1806
1217
        s3 = _mm_add_si64(x0, y0);
1807
1218
        s2 = _mm_add_si64(w14, w16);
1808
1219
        s1 = _mm_add_si64(s1, s3);
1809
1220
        s1 = _mm_add_si64(s1, s2);
1810
 
        C[4] = _mm_cvtsi64_si32(s1);
1811
 
        s1 = _mm_srli_si64(s1, 32);
 
1221
        C[4] = _m_to_int(s1);
 
1222
        s1 = _m_psrlqi(s1, 32);
1812
1223
 
1813
1224
        s3 = _mm_add_si64(x1, y1);
1814
1225
        s4 = _mm_add_si64(x4, y4);
1816
1227
        s3 = _mm_add_si64(s3, s4);
1817
1228
        s1 = _mm_add_si64(s1, w20);
1818
1229
        s1 = _mm_add_si64(s1, s3);
1819
 
        C[5] = _mm_cvtsi64_si32(s1);
1820
 
        s1 = _mm_srli_si64(s1, 32);
 
1230
        C[5] = _m_to_int(s1);
 
1231
        s1 = _m_psrlqi(s1, 32);
1821
1232
 
1822
1233
        s3 = _mm_add_si64(x6, y6);
1823
1234
        s4 = _mm_add_si64(x8, y8);
1825
1236
        s3 = _mm_add_si64(s3, s4);
1826
1237
        s1 = _mm_add_si64(s1, w26);
1827
1238
        s1 = _mm_add_si64(s1, s3);
1828
 
        C[6] = _mm_cvtsi64_si32(s1);
1829
 
        s1 = _mm_srli_si64(s1, 32);
 
1239
        C[6] = _m_to_int(s1);
 
1240
        s1 = _m_psrlqi(s1, 32);
1830
1241
 
1831
1242
        s3 = _mm_add_si64(x10, y10);
1832
1243
        s4 = _mm_add_si64(x12, y12);
1833
1244
        s1 = _mm_add_si64(s1, w27);
1834
1245
        s3 = _mm_add_si64(s3, s4);
1835
1246
        s1 = _mm_add_si64(s1, s3);
1836
 
        C[7] = _mm_cvtsi64_si32(s1);
1837
 
        s1 = _mm_srli_si64(s1, 32);
 
1247
        C[7] = _m_to_int(s1);
 
1248
        s1 = _m_psrlqi(s1, 32);
1838
1249
 
1839
1250
        s3 = _mm_add_si64(x14, y14);
1840
1251
        s4 = _mm_add_si64(x16, y16);
1841
1252
        s1 = _mm_add_si64(s1, z0);
1842
1253
        s3 = _mm_add_si64(s3, s4);
1843
1254
        s1 = _mm_add_si64(s1, s3);
1844
 
        C[8] = _mm_cvtsi64_si32(s1);
1845
 
        s1 = _mm_srli_si64(s1, 32);
 
1255
        C[8] = _m_to_int(s1);
 
1256
        s1 = _m_psrlqi(s1, 32);
1846
1257
 
1847
1258
        s3 = _mm_add_si64(x18, y18);
1848
1259
        s4 = _mm_add_si64(x20, y20);
1850
1261
        s3 = _mm_add_si64(s3, s4);
1851
1262
        s1 = _mm_add_si64(s1, z4);
1852
1263
        s1 = _mm_add_si64(s1, s3);
1853
 
        C[9] = _mm_cvtsi64_si32(s1);
1854
 
        s1 = _mm_srli_si64(s1, 32);
 
1264
        C[9] = _m_to_int(s1);
 
1265
        s1 = _m_psrlqi(s1, 32);
1855
1266
 
1856
1267
        s3 = _mm_add_si64(x22, y22);
1857
1268
        s4 = _mm_add_si64(x26, y26);
1859
1270
        s3 = _mm_add_si64(s3, s4);
1860
1271
        s1 = _mm_add_si64(s1, z8);
1861
1272
        s1 = _mm_add_si64(s1, s3);
1862
 
        C[10] = _mm_cvtsi64_si32(s1);
1863
 
        s1 = _mm_srli_si64(s1, 32);
 
1273
        C[10] = _m_to_int(s1);
 
1274
        s1 = _m_psrlqi(s1, 32);
1864
1275
 
1865
1276
        s3 = _mm_add_si64(x27, y27);
1866
1277
        s1 = _mm_add_si64(s1, z10);
1867
1278
        s1 = _mm_add_si64(s1, z12);
1868
1279
        s1 = _mm_add_si64(s1, s3);
1869
 
        C[11] = _mm_cvtsi64_si32(s1);
1870
 
        s1 = _mm_srli_si64(s1, 32);
 
1280
        C[11] = _m_to_int(s1);
 
1281
        s1 = _m_psrlqi(s1, 32);
1871
1282
 
1872
1283
        s3 = _mm_add_si64(z14, z16);
1873
1284
        s1 = _mm_add_si64(s1, s3);
1874
 
        C[12] = _mm_cvtsi64_si32(s1);
1875
 
        s1 = _mm_srli_si64(s1, 32);
 
1285
        C[12] = _m_to_int(s1);
 
1286
        s1 = _m_psrlqi(s1, 32);
1876
1287
 
1877
1288
        s3 = _mm_add_si64(z18, z20);
1878
1289
        s1 = _mm_add_si64(s1, s3);
1879
 
        C[13] = _mm_cvtsi64_si32(s1);
1880
 
        s1 = _mm_srli_si64(s1, 32);
 
1290
        C[13] = _m_to_int(s1);
 
1291
        s1 = _m_psrlqi(s1, 32);
1881
1292
 
1882
1293
        s3 = _mm_add_si64(z22, z26);
1883
1294
        s1 = _mm_add_si64(s1, s3);
1884
 
        C[14] = _mm_cvtsi64_si32(s1);
1885
 
        s1 = _mm_srli_si64(s1, 32);
 
1295
        C[14] = _m_to_int(s1);
 
1296
        s1 = _m_psrlqi(s1, 32);
1886
1297
 
1887
 
        C[15] = z[27] + _mm_cvtsi64_si32(s1);
 
1298
        C[15] = z[27] + _m_to_int(s1);
1888
1299
        _mm_empty();
1889
1300
}
1890
1301
 
1908
1319
 
1909
1320
        __m64 s1, s2, s3, s4;
1910
1321
 
1911
 
        __m64 w1 = _mm_cvtsi32_si64(w[1]);
 
1322
        __m64 w1 = _m_from_int(w[1]);
1912
1323
        __m64 w4 = mw[2];
1913
1324
        __m64 w6 = mw[3];
1914
1325
        __m64 w8 = mw[4];
1919
1330
        __m64 w18 = mw[9];
1920
1331
        __m64 w20 = mw[10];
1921
1332
        __m64 w22 = mw[11];
1922
 
        __m64 w26 = _mm_cvtsi32_si64(w[26]);
 
1333
        __m64 w26 = _m_from_int(w[26]);
1923
1334
 
1924
 
        __m64 x0 = _mm_cvtsi32_si64(x[0]);
1925
 
        __m64 x1 = _mm_cvtsi32_si64(x[1]);
 
1335
        __m64 x0 = _m_from_int(x[0]);
 
1336
        __m64 x1 = _m_from_int(x[1]);
1926
1337
        __m64 x4 = mx[2];
1927
1338
        __m64 x6 = mx[3];
1928
1339
        __m64 x8 = mx[4];
1929
1340
 
1930
 
        __m64 y0 = _mm_cvtsi32_si64(y[0]);
1931
 
        __m64 y1 = _mm_cvtsi32_si64(y[1]);
 
1341
        __m64 y0 = _m_from_int(y[0]);
 
1342
        __m64 y1 = _m_from_int(y[1]);
1932
1343
        __m64 y4 = my[2];
1933
1344
        __m64 y6 = my[3];
1934
1345
        __m64 y8 = my[4];
1935
1346
 
1936
1347
        s1 = _mm_add_si64(w1, w4);
1937
 
        C[1] = _mm_cvtsi64_si32(s1);
1938
 
        s1 = _mm_srli_si64(s1, 32);
 
1348
        C[1] = _m_to_int(s1);
 
1349
        s1 = _m_psrlqi(s1, 32);
1939
1350
 
1940
1351
        s2 = _mm_add_si64(w6, w8);
1941
1352
        s1 = _mm_add_si64(s1, s2);
1942
 
        C[2] = _mm_cvtsi64_si32(s1);
1943
 
        s1 = _mm_srli_si64(s1, 32);
 
1353
        C[2] = _m_to_int(s1);
 
1354
        s1 = _m_psrlqi(s1, 32);
1944
1355
 
1945
1356
        s2 = _mm_add_si64(w10, w12);
1946
1357
        s1 = _mm_add_si64(s1, s2);
1947
 
        C[3] = _mm_cvtsi64_si32(s1);
1948
 
        s1 = _mm_srli_si64(s1, 32);
 
1358
        C[3] = _m_to_int(s1);
 
1359
        s1 = _m_psrlqi(s1, 32);
1949
1360
 
1950
1361
        s3 = _mm_add_si64(x0, y0);
1951
1362
        s2 = _mm_add_si64(w14, w16);
1952
1363
        s1 = _mm_add_si64(s1, s3);
1953
1364
        s1 = _mm_add_si64(s1, s2);
1954
 
        C[4] = _mm_cvtsi64_si32(s1);
1955
 
        s1 = _mm_srli_si64(s1, 32);
 
1365
        C[4] = _m_to_int(s1);
 
1366
        s1 = _m_psrlqi(s1, 32);
1956
1367
 
1957
1368
        s3 = _mm_add_si64(x1, y1);
1958
1369
        s4 = _mm_add_si64(x4, y4);
1960
1371
        s3 = _mm_add_si64(s3, s4);
1961
1372
        s1 = _mm_add_si64(s1, w20);
1962
1373
        s1 = _mm_add_si64(s1, s3);
1963
 
        C[5] = _mm_cvtsi64_si32(s1);
1964
 
        s1 = _mm_srli_si64(s1, 32);
 
1374
        C[5] = _m_to_int(s1);
 
1375
        s1 = _m_psrlqi(s1, 32);
1965
1376
 
1966
1377
        s3 = _mm_add_si64(x6, y6);
1967
1378
        s4 = _mm_add_si64(x8, y8);
1969
1380
        s3 = _mm_add_si64(s3, s4);
1970
1381
        s1 = _mm_add_si64(s1, w26);
1971
1382
        s1 = _mm_add_si64(s1, s3);
1972
 
        C[6] = _mm_cvtsi64_si32(s1);
1973
 
        s1 = _mm_srli_si64(s1, 32);
 
1383
        C[6] = _m_to_int(s1);
 
1384
        s1 = _m_psrlqi(s1, 32);
1974
1385
 
1975
 
        C[7] = _mm_cvtsi64_si32(s1) + w[27] + x[10] + y[10] + x[12] + y[12];
 
1386
        C[7] = _m_to_int(s1) + w[27] + x[10] + y[10] + x[12] + y[12];
1976
1387
        _mm_empty();
1977
1388
}
1978
1389
 
 
1390
__declspec(naked) word __fastcall P4Optimized::Add(word *C, const word *A, const word *B, unsigned int N)
 
1391
{
 
1392
        __asm
 
1393
        {
 
1394
                sub             esp, 16
 
1395
                xor             eax, eax
 
1396
                mov             [esp], edi
 
1397
                mov             [esp+4], esi
 
1398
                mov             [esp+8], ebx
 
1399
                mov             [esp+12], ebp
 
1400
 
 
1401
                mov             ebx, [esp+20]   // B
 
1402
                mov             esi, [esp+24]   // N
 
1403
 
 
1404
                // now: ebx = B, ecx = C, edx = A, esi = N
 
1405
 
 
1406
                neg             esi
 
1407
                jz              loopend         // if no dwords then nothing to do
 
1408
 
 
1409
                mov             edi, [edx]
 
1410
                mov             ebp, [ebx]
 
1411
 
 
1412
loopstart:
 
1413
                add             edi, eax
 
1414
                jc              carry1
 
1415
 
 
1416
                xor             eax, eax
 
1417
 
 
1418
carry1continue:
 
1419
                add             edi, ebp
 
1420
                mov             ebp, 1
 
1421
                mov             [ecx], edi
 
1422
                mov             edi, [edx+4]
 
1423
                cmovc   eax, ebp
 
1424
                mov             ebp, [ebx+4]
 
1425
                lea             ebx, [ebx+8]
 
1426
                add             edi, eax
 
1427
                jc              carry2
 
1428
 
 
1429
                xor             eax, eax
 
1430
 
 
1431
carry2continue:
 
1432
                add             edi, ebp
 
1433
                mov             ebp, 1
 
1434
                cmovc   eax, ebp
 
1435
                mov             [ecx+4], edi
 
1436
                add             ecx, 8
 
1437
                mov             edi, [edx+8]
 
1438
                add             edx, 8
 
1439
                add             esi, 2
 
1440
                mov             ebp, [ebx]
 
1441
                jnz             loopstart
 
1442
 
 
1443
loopend:
 
1444
                mov             edi, [esp]
 
1445
                mov             esi, [esp+4]
 
1446
                mov             ebx, [esp+8]
 
1447
                mov             ebp, [esp+12]
 
1448
                add             esp, 16
 
1449
                ret             8
 
1450
 
 
1451
carry1:
 
1452
                mov             eax, 1
 
1453
                jmp             carry1continue
 
1454
 
 
1455
carry2:
 
1456
                mov             eax, 1
 
1457
                jmp             carry2continue
 
1458
        }
 
1459
}
 
1460
 
 
1461
__declspec(naked) word __fastcall P4Optimized::Subtract(word *C, const word *A, const word *B, unsigned int N)
 
1462
{
 
1463
        __asm
 
1464
        {
 
1465
                sub             esp, 16
 
1466
                xor             eax, eax
 
1467
                mov             [esp], edi
 
1468
                mov             [esp+4], esi
 
1469
                mov             [esp+8], ebx
 
1470
                mov             [esp+12], ebp
 
1471
 
 
1472
                mov             ebx, [esp+20]   // B
 
1473
                mov             esi, [esp+24]   // N
 
1474
 
 
1475
                // now: ebx = B, ecx = C, edx = A, esi = N
 
1476
 
 
1477
                neg             esi
 
1478
                jz              loopend         // if no dwords then nothing to do
 
1479
 
 
1480
                mov             edi, [edx]
 
1481
                mov             ebp, [ebx]
 
1482
 
 
1483
loopstart:
 
1484
                sub             edi, eax
 
1485
                jc              carry1
 
1486
 
 
1487
                xor             eax, eax
 
1488
 
 
1489
carry1continue:
 
1490
                sub             edi, ebp
 
1491
                mov             ebp, 1
 
1492
                mov             [ecx], edi
 
1493
                mov             edi, [edx+4]
 
1494
                cmovc   eax, ebp
 
1495
                mov             ebp, [ebx+4]
 
1496
                lea             ebx, [ebx+8]
 
1497
                sub             edi, eax
 
1498
                jc              carry2
 
1499
 
 
1500
                xor             eax, eax
 
1501
 
 
1502
carry2continue:
 
1503
                sub             edi, ebp
 
1504
                mov             ebp, 1
 
1505
                cmovc   eax, ebp
 
1506
                mov             [ecx+4], edi
 
1507
                add             ecx, 8
 
1508
                mov             edi, [edx+8]
 
1509
                add             edx, 8
 
1510
                add             esi, 2
 
1511
                mov             ebp, [ebx]
 
1512
                jnz             loopstart
 
1513
 
 
1514
loopend:
 
1515
                mov             edi, [esp]
 
1516
                mov             esi, [esp+4]
 
1517
                mov             ebx, [esp+8]
 
1518
                mov             ebp, [esp+12]
 
1519
                add             esp, 16
 
1520
                ret             8
 
1521
 
 
1522
carry1:
 
1523
                mov             eax, 1
 
1524
                jmp             carry1continue
 
1525
 
 
1526
carry2:
 
1527
                mov             eax, 1
 
1528
                jmp             carry2continue
 
1529
        }
 
1530
}
 
1531
 
1979
1532
#endif  // #ifdef SSE2_INTRINSICS_AVAILABLE
1980
1533
 
 
1534
#elif defined(__GNUC__) && defined(__i386__)
 
1535
 
 
1536
class PentiumOptimized : public Portable
 
1537
{
 
1538
public:
 
1539
#ifndef __pic__         // -fpic uses up a register, leaving too few for the asm code
 
1540
        static word Add(word *C, const word *A, const word *B, unsigned int N);
 
1541
        static word Subtract(word *C, const word *A, const word *B, unsigned int N);
 
1542
#endif
 
1543
        static void Square4(word *R, const word *A);
 
1544
        static void Multiply4(word *C, const word *A, const word *B);
 
1545
        static void Multiply8(word *C, const word *A, const word *B);
 
1546
};
 
1547
 
 
1548
typedef PentiumOptimized LowLevel;
 
1549
 
 
1550
// Add and Subtract assembly code originally contributed by Alister Lee
 
1551
 
 
1552
#ifndef __pic__
 
1553
__attribute__((regparm(3))) word PentiumOptimized::Add(word *C, const word *A, const word *B, unsigned int N)
 
1554
{
 
1555
        assert (N%2 == 0);
 
1556
 
 
1557
        register word carry, temp;
 
1558
 
 
1559
        __asm__ __volatile__(
 
1560
                        "push %%ebp;"
 
1561
                        "sub %3, %2;"
 
1562
                        "xor %0, %0;"
 
1563
                        "sub %4, %0;"
 
1564
                        "lea (%1,%4,4), %1;"
 
1565
                        "sar $1, %0;"
 
1566
                        "jz 1f;"
 
1567
 
 
1568
                "0:;"
 
1569
                        "mov 0(%3), %4;"
 
1570
                        "mov 4(%3), %%ebp;"
 
1571
                        "mov (%1,%0,8), %5;"
 
1572
                        "lea 8(%3), %3;"
 
1573
                        "adc %5, %4;"
 
1574
                        "mov 4(%1,%0,8), %5;"
 
1575
                        "adc %5, %%ebp;"
 
1576
                        "inc %0;"
 
1577
                        "mov %4, -8(%3, %2);"
 
1578
                        "mov %%ebp, -4(%3, %2);"
 
1579
                        "jnz 0b;"
 
1580
 
 
1581
                "1:;"
 
1582
                        "adc $0, %0;"
 
1583
                        "pop %%ebp;"
 
1584
 
 
1585
                : "=aSD" (carry), "+r" (B), "+r" (C), "+r" (A), "+r" (N), "=r" (temp)
 
1586
                : : "cc", "memory");
 
1587
 
 
1588
        return carry;
 
1589
}
 
1590
 
 
1591
__attribute__((regparm(3))) word PentiumOptimized::Subtract(word *C, const word *A, const word *B, unsigned int N)
 
1592
{
 
1593
        assert (N%2 == 0);
 
1594
 
 
1595
        register word carry, temp;
 
1596
 
 
1597
        __asm__ __volatile__(
 
1598
                        "push %%ebp;"
 
1599
                        "sub %3, %2;"
 
1600
                        "xor %0, %0;"
 
1601
                        "sub %4, %0;"
 
1602
                        "lea (%1,%4,4), %1;"
 
1603
                        "sar $1, %0;"
 
1604
                        "jz 1f;"
 
1605
 
 
1606
                "0:;"
 
1607
                        "mov 0(%3), %4;"
 
1608
                        "mov 4(%3), %%ebp;"
 
1609
                        "mov (%1,%0,8), %5;"
 
1610
                        "lea 8(%3), %3;"
 
1611
                        "sbb %5, %4;"
 
1612
                        "mov 4(%1,%0,8), %5;"
 
1613
                        "sbb %5, %%ebp;"
 
1614
                        "inc %0;"
 
1615
                        "mov %4, -8(%3, %2);"
 
1616
                        "mov %%ebp, -4(%3, %2);"
 
1617
                        "jnz 0b;"
 
1618
 
 
1619
                "1:;"
 
1620
                        "adc $0, %0;"
 
1621
                        "pop %%ebp;"
 
1622
 
 
1623
                : "=aSD" (carry), "+r" (B), "+r" (C), "+r" (A), "+r" (N), "=r" (temp)
 
1624
                : : "cc", "memory");
 
1625
 
 
1626
        return carry;
 
1627
}
 
1628
#endif  // __pic__
 
1629
 
 
1630
// Comba square and multiply assembly code originally contributed by Leonard Janke
 
1631
 
 
1632
#define SqrStartup \
 
1633
  "push %%ebp\n\t" \
 
1634
  "push %%esi\n\t" \
 
1635
  "push %%ebx\n\t" \
 
1636
  "xor %%ebp, %%ebp\n\t" \
 
1637
  "xor %%ebx, %%ebx\n\t" \
 
1638
  "xor %%ecx, %%ecx\n\t" 
 
1639
 
 
1640
#define SqrShiftCarry \
 
1641
  "mov %%ebx, %%ebp\n\t" \
 
1642
  "mov %%ecx, %%ebx\n\t" \
 
1643
  "xor %%ecx, %%ecx\n\t"
 
1644
 
 
1645
#define SqrAccumulate(i,j) \
 
1646
  "mov 4*"#j"(%%esi), %%eax\n\t" \
 
1647
  "mull 4*"#i"(%%esi)\n\t" \
 
1648
  "add %%eax, %%ebp\n\t" \
 
1649
  "adc %%edx, %%ebx\n\t" \
 
1650
  "adc %%ch, %%cl\n\t" \
 
1651
  "add %%eax, %%ebp\n\t" \
 
1652
  "adc %%edx, %%ebx\n\t" \
 
1653
  "adc %%ch, %%cl\n\t"
 
1654
 
 
1655
#define SqrAccumulateCentre(i) \
 
1656
  "mov 4*"#i"(%%esi), %%eax\n\t" \
 
1657
  "mull 4*"#i"(%%esi)\n\t" \
 
1658
  "add %%eax, %%ebp\n\t" \
 
1659
  "adc %%edx, %%ebx\n\t" \
 
1660
  "adc %%ch, %%cl\n\t" 
 
1661
 
 
1662
#define SqrStoreDigit(X)  \
 
1663
  "mov %%ebp, 4*"#X"(%%edi)\n\t" \
 
1664
 
 
1665
#define SqrLastDiagonal(digits) \
 
1666
  "mov 4*("#digits"-1)(%%esi), %%eax\n\t" \
 
1667
  "mull 4*("#digits"-1)(%%esi)\n\t" \
 
1668
  "add %%eax, %%ebp\n\t" \
 
1669
  "adc %%edx, %%ebx\n\t" \
 
1670
  "mov %%ebp, 4*(2*"#digits"-2)(%%edi)\n\t" \
 
1671
  "mov %%ebx, 4*(2*"#digits"-1)(%%edi)\n\t" 
 
1672
 
 
1673
#define SqrCleanup \
 
1674
  "pop %%ebx\n\t" \
 
1675
  "pop %%esi\n\t" \
 
1676
  "pop %%ebp\n\t" 
 
1677
 
 
1678
void PentiumOptimized::Square4(word* Y, const word* X)
 
1679
{
 
1680
        __asm__ __volatile__(
 
1681
                SqrStartup
 
1682
 
 
1683
                SqrAccumulateCentre(0)
 
1684
                SqrStoreDigit(0)
 
1685
                SqrShiftCarry
 
1686
 
 
1687
                SqrAccumulate(1,0)
 
1688
                SqrStoreDigit(1)
 
1689
                SqrShiftCarry
 
1690
 
 
1691
                SqrAccumulate(2,0)
 
1692
                SqrAccumulateCentre(1)
 
1693
                SqrStoreDigit(2)
 
1694
                SqrShiftCarry
 
1695
 
 
1696
                SqrAccumulate(3,0)
 
1697
                SqrAccumulate(2,1)
 
1698
                SqrStoreDigit(3)
 
1699
                SqrShiftCarry
 
1700
 
 
1701
                SqrAccumulate(3,1)
 
1702
                SqrAccumulateCentre(2)
 
1703
                SqrStoreDigit(4)
 
1704
                SqrShiftCarry
 
1705
 
 
1706
                SqrAccumulate(3,2)
 
1707
                SqrStoreDigit(5)
 
1708
                SqrShiftCarry
 
1709
 
 
1710
                SqrLastDiagonal(4)
 
1711
 
 
1712
                SqrCleanup
 
1713
 
 
1714
                :
 
1715
                : "D" (Y), "S" (X)
 
1716
                : "eax",  "ecx", "edx", "ebp",   "memory"
 
1717
        );
 
1718
}
 
1719
 
 
1720
#define MulStartup \
 
1721
  "push %%ebp\n\t" \
 
1722
  "push %%esi\n\t" \
 
1723
  "push %%ebx\n\t" \
 
1724
  "push %%edi\n\t" \
 
1725
  "mov %%eax, %%ebx \n\t" \
 
1726
  "xor %%ebp, %%ebp\n\t" \
 
1727
  "xor %%edi, %%edi\n\t" \
 
1728
  "xor %%ecx, %%ecx\n\t" 
 
1729
 
 
1730
#define MulShiftCarry \
 
1731
  "mov %%edx, %%ebp\n\t" \
 
1732
  "mov %%ecx, %%edi\n\t" \
 
1733
  "xor %%ecx, %%ecx\n\t"
 
1734
 
 
1735
#define MulAccumulate(i,j) \
 
1736
  "mov 4*"#j"(%%ebx), %%eax\n\t" \
 
1737
  "mull 4*"#i"(%%esi)\n\t" \
 
1738
  "add %%eax, %%ebp\n\t" \
 
1739
  "adc %%edx, %%edi\n\t" \
 
1740
  "adc %%ch, %%cl\n\t"
 
1741
 
 
1742
#define MulStoreDigit(X)  \
 
1743
  "mov %%edi, %%edx \n\t" \
 
1744
  "mov (%%esp), %%edi \n\t" \
 
1745
  "mov %%ebp, 4*"#X"(%%edi)\n\t" \
 
1746
  "mov %%edi, (%%esp)\n\t" 
 
1747
 
 
1748
#define MulLastDiagonal(digits) \
 
1749
  "mov 4*("#digits"-1)(%%ebx), %%eax\n\t" \
 
1750
  "mull 4*("#digits"-1)(%%esi)\n\t" \
 
1751
  "add %%eax, %%ebp\n\t" \
 
1752
  "adc %%edi, %%edx\n\t" \
 
1753
  "mov (%%esp), %%edi\n\t" \
 
1754
  "mov %%ebp, 4*(2*"#digits"-2)(%%edi)\n\t" \
 
1755
  "mov %%edx, 4*(2*"#digits"-1)(%%edi)\n\t" 
 
1756
 
 
1757
#define MulCleanup \
 
1758
  "pop %%edi\n\t" \
 
1759
  "pop %%ebx\n\t" \
 
1760
  "pop %%esi\n\t" \
 
1761
  "pop %%ebp\n\t" 
 
1762
 
 
1763
void PentiumOptimized::Multiply4(word* Z, const word* X, const word* Y)
 
1764
{
 
1765
        __asm__ __volatile__(
 
1766
                MulStartup
 
1767
                MulAccumulate(0,0)
 
1768
                MulStoreDigit(0)
 
1769
                MulShiftCarry
 
1770
 
 
1771
                MulAccumulate(1,0)
 
1772
                MulAccumulate(0,1)
 
1773
                MulStoreDigit(1)
 
1774
                MulShiftCarry
 
1775
 
 
1776
                MulAccumulate(2,0)
 
1777
                MulAccumulate(1,1)
 
1778
                MulAccumulate(0,2)
 
1779
                MulStoreDigit(2)
 
1780
                MulShiftCarry
 
1781
 
 
1782
                MulAccumulate(3,0)
 
1783
                MulAccumulate(2,1)
 
1784
                MulAccumulate(1,2)
 
1785
                MulAccumulate(0,3)
 
1786
                MulStoreDigit(3)
 
1787
                MulShiftCarry
 
1788
 
 
1789
                MulAccumulate(3,1)
 
1790
                MulAccumulate(2,2)
 
1791
                MulAccumulate(1,3)
 
1792
                MulStoreDigit(4)
 
1793
                MulShiftCarry
 
1794
 
 
1795
                MulAccumulate(3,2)
 
1796
                MulAccumulate(2,3)
 
1797
                MulStoreDigit(5)
 
1798
                MulShiftCarry
 
1799
 
 
1800
                MulLastDiagonal(4)
 
1801
 
 
1802
                MulCleanup
 
1803
 
 
1804
                : 
 
1805
                : "D" (Z), "S" (X), "a" (Y)
 
1806
                : "%ecx", "%edx",  "memory"
 
1807
        );
 
1808
}
 
1809
 
 
1810
void PentiumOptimized::Multiply8(word* Z, const word* X, const word* Y)
 
1811
{
 
1812
        __asm__ __volatile__(
 
1813
                MulStartup
 
1814
                MulAccumulate(0,0)
 
1815
                MulStoreDigit(0)
 
1816
                MulShiftCarry
 
1817
 
 
1818
                MulAccumulate(1,0)
 
1819
                MulAccumulate(0,1)
 
1820
                MulStoreDigit(1)
 
1821
                MulShiftCarry
 
1822
 
 
1823
                MulAccumulate(2,0)
 
1824
                MulAccumulate(1,1)
 
1825
                MulAccumulate(0,2)
 
1826
                MulStoreDigit(2)
 
1827
                MulShiftCarry
 
1828
 
 
1829
                MulAccumulate(3,0)
 
1830
                MulAccumulate(2,1)
 
1831
                MulAccumulate(1,2)
 
1832
                MulAccumulate(0,3)
 
1833
                MulStoreDigit(3)
 
1834
                MulShiftCarry
 
1835
 
 
1836
                MulAccumulate(4,0)
 
1837
                MulAccumulate(3,1)
 
1838
                MulAccumulate(2,2)
 
1839
                MulAccumulate(1,3)
 
1840
                MulAccumulate(0,4)
 
1841
                MulStoreDigit(4)
 
1842
                MulShiftCarry
 
1843
 
 
1844
                MulAccumulate(5,0)
 
1845
                MulAccumulate(4,1)
 
1846
                MulAccumulate(3,2)
 
1847
                MulAccumulate(2,3)
 
1848
                MulAccumulate(1,4)
 
1849
                MulAccumulate(0,5)
 
1850
                MulStoreDigit(5)
 
1851
                MulShiftCarry
 
1852
 
 
1853
                MulAccumulate(6,0)
 
1854
                MulAccumulate(5,1)
 
1855
                MulAccumulate(4,2)
 
1856
                MulAccumulate(3,3)
 
1857
                MulAccumulate(2,4)
 
1858
                MulAccumulate(1,5)
 
1859
                MulAccumulate(0,6)
 
1860
                MulStoreDigit(6)
 
1861
                MulShiftCarry
 
1862
 
 
1863
                MulAccumulate(7,0)
 
1864
                MulAccumulate(6,1)
 
1865
                MulAccumulate(5,2)
 
1866
                MulAccumulate(4,3)
 
1867
                MulAccumulate(3,4)
 
1868
                MulAccumulate(2,5)
 
1869
                MulAccumulate(1,6)
 
1870
                MulAccumulate(0,7)
 
1871
                MulStoreDigit(7)
 
1872
                MulShiftCarry
 
1873
 
 
1874
                MulAccumulate(7,1)
 
1875
                MulAccumulate(6,2)
 
1876
                MulAccumulate(5,3)
 
1877
                MulAccumulate(4,4)
 
1878
                MulAccumulate(3,5)
 
1879
                MulAccumulate(2,6)
 
1880
                MulAccumulate(1,7)
 
1881
                MulStoreDigit(8)
 
1882
                MulShiftCarry
 
1883
 
 
1884
                MulAccumulate(7,2)
 
1885
                MulAccumulate(6,3)
 
1886
                MulAccumulate(5,4)
 
1887
                MulAccumulate(4,5)
 
1888
                MulAccumulate(3,6)
 
1889
                MulAccumulate(2,7)
 
1890
                MulStoreDigit(9)
 
1891
                MulShiftCarry
 
1892
 
 
1893
                MulAccumulate(7,3)
 
1894
                MulAccumulate(6,4)
 
1895
                MulAccumulate(5,5)
 
1896
                MulAccumulate(4,6)
 
1897
                MulAccumulate(3,7)
 
1898
                MulStoreDigit(10)
 
1899
                MulShiftCarry
 
1900
 
 
1901
                MulAccumulate(7,4)
 
1902
                MulAccumulate(6,5)
 
1903
                MulAccumulate(5,6)
 
1904
                MulAccumulate(4,7)
 
1905
                MulStoreDigit(11)
 
1906
                MulShiftCarry
 
1907
 
 
1908
                MulAccumulate(7,5)
 
1909
                MulAccumulate(6,6)
 
1910
                MulAccumulate(5,7)
 
1911
                MulStoreDigit(12)
 
1912
                MulShiftCarry
 
1913
 
 
1914
                MulAccumulate(7,6)
 
1915
                MulAccumulate(6,7)
 
1916
                MulStoreDigit(13)
 
1917
                MulShiftCarry
 
1918
 
 
1919
                MulLastDiagonal(8)
 
1920
 
 
1921
                MulCleanup
 
1922
 
 
1923
                : 
 
1924
                : "D" (Z), "S" (X), "a" (Y)
 
1925
                : "%ecx", "%edx",  "memory"
 
1926
        );
 
1927
}
 
1928
 
 
1929
#else   // no processor specific code at this layer
 
1930
 
 
1931
typedef Portable LowLevel;
 
1932
 
 
1933
#endif
 
1934
 
1981
1935
// ********************************************************
1982
1936
 
1983
1937
#define A0              A
1995
1949
#define R2              (R+N)
1996
1950
#define R3              (R+N+N2)
1997
1951
 
 
1952
//VC60 workaround: compiler bug triggered without the extra dummy parameters
 
1953
 
1998
1954
// R[2*N] - result = A*B
1999
1955
// T[2*N] - temporary work space
2000
1956
// A[N] --- multiplier
2001
1957
// B[N] --- multiplicant
2002
1958
 
2003
 
void RecursiveMultiply(word *R, word *T, const word *A, const word *B, size_t N)
 
1959
template <class P>
 
1960
void DoRecursiveMultiply(word *R, word *T, const word *A, const word *B, unsigned int N, const P *dummy=NULL);
 
1961
 
 
1962
template <class P>
 
1963
inline void RecursiveMultiply(word *R, word *T, const word *A, const word *B, unsigned int N, const P *dummy=NULL)
2004
1964
{
2005
1965
        assert(N>=2 && N%2==0);
2006
1966
 
2007
 
        if (LowLevel::MultiplyRecursionLimit() >= 8 && N==8)
2008
 
                LowLevel::Multiply8(R, A, B);
2009
 
        else if (LowLevel::MultiplyRecursionLimit() >= 4 && N==4)
2010
 
                LowLevel::Multiply4(R, A, B);
 
1967
        if (P::MultiplyRecursionLimit() >= 8 && N==8)
 
1968
                P::Multiply8(R, A, B);
 
1969
        else if (P::MultiplyRecursionLimit() >= 4 && N==4)
 
1970
                P::Multiply4(R, A, B);
2011
1971
        else if (N==2)
2012
 
                LowLevel::Multiply2(R, A, B);
 
1972
                P::Multiply2(R, A, B);
2013
1973
        else
 
1974
                DoRecursiveMultiply<P>(R, T, A, B, N, NULL);    // VC60 workaround: needs this NULL
 
1975
}
 
1976
 
 
1977
template <class P>
 
1978
void DoRecursiveMultiply(word *R, word *T, const word *A, const word *B, unsigned int N, const P *dummy)
 
1979
{
 
1980
        const unsigned int N2 = N/2;
 
1981
        int carry;
 
1982
 
 
1983
        int aComp = Compare(A0, A1, N2);
 
1984
        int bComp = Compare(B0, B1, N2);
 
1985
 
 
1986
        switch (2*aComp + aComp + bComp)
2014
1987
        {
2015
 
                const size_t N2 = N/2;
2016
 
                int carry;
2017
 
 
2018
 
                int aComp = Compare(A0, A1, N2);
2019
 
                int bComp = Compare(B0, B1, N2);
2020
 
 
2021
 
                switch (2*aComp + aComp + bComp)
2022
 
                {
2023
 
                case -4:
2024
 
                        LowLevel::Subtract(R0, A1, A0, N2);
2025
 
                        LowLevel::Subtract(R1, B0, B1, N2);
2026
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
2027
 
                        LowLevel::Subtract(T1, T1, R0, N2);
2028
 
                        carry = -1;
2029
 
                        break;
2030
 
                case -2:
2031
 
                        LowLevel::Subtract(R0, A1, A0, N2);
2032
 
                        LowLevel::Subtract(R1, B0, B1, N2);
2033
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
2034
 
                        carry = 0;
2035
 
                        break;
2036
 
                case 2:
2037
 
                        LowLevel::Subtract(R0, A0, A1, N2);
2038
 
                        LowLevel::Subtract(R1, B1, B0, N2);
2039
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
2040
 
                        carry = 0;
2041
 
                        break;
2042
 
                case 4:
2043
 
                        LowLevel::Subtract(R0, A1, A0, N2);
2044
 
                        LowLevel::Subtract(R1, B0, B1, N2);
2045
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
2046
 
                        LowLevel::Subtract(T1, T1, R1, N2);
2047
 
                        carry = -1;
2048
 
                        break;
2049
 
                default:
2050
 
                        SetWords(T0, 0, N);
2051
 
                        carry = 0;
2052
 
                }
2053
 
 
2054
 
                RecursiveMultiply(R0, T2, A0, B0, N2);
2055
 
                RecursiveMultiply(R2, T2, A1, B1, N2);
2056
 
 
2057
 
                // now T[01] holds (A1-A0)*(B0-B1), R[01] holds A0*B0, R[23] holds A1*B1
2058
 
 
2059
 
                carry += LowLevel::Add(T0, T0, R0, N);
2060
 
                carry += LowLevel::Add(T0, T0, R2, N);
2061
 
                carry += LowLevel::Add(R1, R1, T0, N);
2062
 
 
2063
 
                assert (carry >= 0 && carry <= 2);
2064
 
                Increment(R3, N2, carry);
 
1988
        case -4:
 
1989
                P::Subtract(R0, A1, A0, N2);
 
1990
                P::Subtract(R1, B0, B1, N2);
 
1991
                RecursiveMultiply<P>(T0, T2, R0, R1, N2);
 
1992
                P::Subtract(T1, T1, R0, N2);
 
1993
                carry = -1;
 
1994
                break;
 
1995
        case -2:
 
1996
                P::Subtract(R0, A1, A0, N2);
 
1997
                P::Subtract(R1, B0, B1, N2);
 
1998
                RecursiveMultiply<P>(T0, T2, R0, R1, N2);
 
1999
                carry = 0;
 
2000
                break;
 
2001
        case 2:
 
2002
                P::Subtract(R0, A0, A1, N2);
 
2003
                P::Subtract(R1, B1, B0, N2);
 
2004
                RecursiveMultiply<P>(T0, T2, R0, R1, N2);
 
2005
                carry = 0;
 
2006
                break;
 
2007
        case 4:
 
2008
                P::Subtract(R0, A1, A0, N2);
 
2009
                P::Subtract(R1, B0, B1, N2);
 
2010
                RecursiveMultiply<P>(T0, T2, R0, R1, N2);
 
2011
                P::Subtract(T1, T1, R1, N2);
 
2012
                carry = -1;
 
2013
                break;
 
2014
        default:
 
2015
                SetWords(T0, 0, N);
 
2016
                carry = 0;
2065
2017
        }
 
2018
 
 
2019
        RecursiveMultiply<P>(R0, T2, A0, B0, N2);
 
2020
        RecursiveMultiply<P>(R2, T2, A1, B1, N2);
 
2021
 
 
2022
        // now T[01] holds (A1-A0)*(B0-B1), R[01] holds A0*B0, R[23] holds A1*B1
 
2023
 
 
2024
        carry += P::Add(T0, T0, R0, N);
 
2025
        carry += P::Add(T0, T0, R2, N);
 
2026
        carry += P::Add(R1, R1, T0, N);
 
2027
 
 
2028
        assert (carry >= 0 && carry <= 2);
 
2029
        Increment(R3, N2, carry);
2066
2030
}
2067
2031
 
2068
2032
// R[2*N] - result = A*A
2069
2033
// T[2*N] - temporary work space
2070
2034
// A[N] --- number to be squared
2071
2035
 
2072
 
void RecursiveSquare(word *R, word *T, const word *A, size_t N)
 
2036
template <class P>
 
2037
void DoRecursiveSquare(word *R, word *T, const word *A, unsigned int N, const P *dummy=NULL);
 
2038
 
 
2039
template <class P>
 
2040
inline void RecursiveSquare(word *R, word *T, const word *A, unsigned int N, const P *dummy=NULL)
2073
2041
{
2074
2042
        assert(N && N%2==0);
2075
 
        if (LowLevel::SquareRecursionLimit() >= 8 && N==8)
2076
 
                LowLevel::Square8(R, A);
2077
 
        if (LowLevel::SquareRecursionLimit() >= 4 && N==4)
2078
 
                LowLevel::Square4(R, A);
 
2043
        if (P::SquareRecursionLimit() >= 8 && N==8)
 
2044
                P::Square8(R, A);
 
2045
        if (P::SquareRecursionLimit() >= 4 && N==4)
 
2046
                P::Square4(R, A);
2079
2047
        else if (N==2)
2080
 
                LowLevel::Square2(R, A);
 
2048
                P::Square2(R, A);
2081
2049
        else
2082
 
        {
2083
 
                const size_t N2 = N/2;
2084
 
 
2085
 
                RecursiveSquare(R0, T2, A0, N2);
2086
 
                RecursiveSquare(R2, T2, A1, N2);
2087
 
                RecursiveMultiply(T0, T2, A0, A1, N2);
2088
 
 
2089
 
                int carry = LowLevel::Add(R1, R1, T0, N);
2090
 
                carry += LowLevel::Add(R1, R1, T0, N);
2091
 
                Increment(R3, N2, carry);
2092
 
        }
 
2050
                DoRecursiveSquare<P>(R, T, A, N, NULL); // VC60 workaround: needs this NULL
 
2051
}
 
2052
 
 
2053
template <class P>
 
2054
void DoRecursiveSquare(word *R, word *T, const word *A, unsigned int N, const P *dummy)
 
2055
{
 
2056
        const unsigned int N2 = N/2;
 
2057
 
 
2058
        RecursiveSquare<P>(R0, T2, A0, N2);
 
2059
        RecursiveSquare<P>(R2, T2, A1, N2);
 
2060
        RecursiveMultiply<P>(T0, T2, A0, A1, N2);
 
2061
 
 
2062
        word carry = P::Add(R1, R1, T0, N);
 
2063
        carry += P::Add(R1, R1, T0, N);
 
2064
        Increment(R3, N2, carry);
2093
2065
}
2094
2066
 
2095
2067
// R[N] - bottom half of A*B
2097
2069
// A[N] - multiplier
2098
2070
// B[N] - multiplicant
2099
2071
 
2100
 
void RecursiveMultiplyBottom(word *R, word *T, const word *A, const word *B, size_t N)
 
2072
template <class P>
 
2073
void DoRecursiveMultiplyBottom(word *R, word *T, const word *A, const word *B, unsigned int N, const P *dummy=NULL);
 
2074
 
 
2075
template <class P>
 
2076
inline void RecursiveMultiplyBottom(word *R, word *T, const word *A, const word *B, unsigned int N, const P *dummy=NULL)
2101
2077
{
2102
2078
        assert(N>=2 && N%2==0);
2103
 
        if (LowLevel::MultiplyBottomRecursionLimit() >= 8 && N==8)
2104
 
                LowLevel::Multiply8Bottom(R, A, B);
2105
 
        else if (LowLevel::MultiplyBottomRecursionLimit() >= 4 && N==4)
2106
 
                LowLevel::Multiply4Bottom(R, A, B);
 
2079
        if (P::MultiplyBottomRecursionLimit() >= 8 && N==8)
 
2080
                P::Multiply8Bottom(R, A, B);
 
2081
        else if (P::MultiplyBottomRecursionLimit() >= 4 && N==4)
 
2082
                P::Multiply4Bottom(R, A, B);
2107
2083
        else if (N==2)
2108
 
                LowLevel::Multiply2Bottom(R, A, B);
 
2084
                P::Multiply2Bottom(R, A, B);
2109
2085
        else
2110
 
        {
2111
 
                const size_t N2 = N/2;
2112
 
 
2113
 
                RecursiveMultiply(R, T, A0, B0, N2);
2114
 
                RecursiveMultiplyBottom(T0, T1, A1, B0, N2);
2115
 
                LowLevel::Add(R1, R1, T0, N2);
2116
 
                RecursiveMultiplyBottom(T0, T1, A0, B1, N2);
2117
 
                LowLevel::Add(R1, R1, T0, N2);
2118
 
        }
 
2086
                DoRecursiveMultiplyBottom<P>(R, T, A, B, N, NULL);
 
2087
}
 
2088
 
 
2089
template <class P>
 
2090
void DoRecursiveMultiplyBottom(word *R, word *T, const word *A, const word *B, unsigned int N, const P *dummy)
 
2091
{
 
2092
        const unsigned int N2 = N/2;
 
2093
 
 
2094
        RecursiveMultiply<P>(R, T, A0, B0, N2);
 
2095
        RecursiveMultiplyBottom<P>(T0, T1, A1, B0, N2);
 
2096
        P::Add(R1, R1, T0, N2);
 
2097
        RecursiveMultiplyBottom<P>(T0, T1, A0, B1, N2);
 
2098
        P::Add(R1, R1, T0, N2);
2119
2099
}
2120
2100
 
2121
2101
// R[N] --- upper half of A*B
2124
2104
// A[N] --- multiplier
2125
2105
// B[N] --- multiplicant
2126
2106
 
2127
 
void RecursiveMultiplyTop(word *R, word *T, const word *L, const word *A, const word *B, size_t N)
 
2107
template <class P>
 
2108
void RecursiveMultiplyTop(word *R, word *T, const word *L, const word *A, const word *B, unsigned int N, const P *dummy=NULL)
2128
2109
{
2129
2110
        assert(N>=2 && N%2==0);
2130
2111
 
2131
2112
        if (N==4)
2132
2113
        {
2133
 
                LowLevel::Multiply4(T, A, B);
 
2114
                P::Multiply4(T, A, B);
2134
2115
                memcpy(R, T+4, 4*WORD_SIZE);
2135
2116
        }
2136
2117
        else if (N==2)
2137
2118
        {
2138
 
                LowLevel::Multiply2(T, A, B);
 
2119
                P::Multiply2(T, A, B);
2139
2120
                memcpy(R, T+2, 2*WORD_SIZE);
2140
2121
        }
2141
2122
        else
2142
2123
        {
2143
 
                const size_t N2 = N/2;
 
2124
                const unsigned int N2 = N/2;
2144
2125
                int carry;
2145
2126
 
2146
2127
                int aComp = Compare(A0, A1, N2);
2149
2130
                switch (2*aComp + aComp + bComp)
2150
2131
                {
2151
2132
                case -4:
2152
 
                        LowLevel::Subtract(R0, A1, A0, N2);
2153
 
                        LowLevel::Subtract(R1, B0, B1, N2);
2154
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
2155
 
                        LowLevel::Subtract(T1, T1, R0, N2);
 
2133
                        P::Subtract(R0, A1, A0, N2);
 
2134
                        P::Subtract(R1, B0, B1, N2);
 
2135
                        RecursiveMultiply<P>(T0, T2, R0, R1, N2);
 
2136
                        P::Subtract(T1, T1, R0, N2);
2156
2137
                        carry = -1;
2157
2138
                        break;
2158
2139
                case -2:
2159
 
                        LowLevel::Subtract(R0, A1, A0, N2);
2160
 
                        LowLevel::Subtract(R1, B0, B1, N2);
2161
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
 
2140
                        P::Subtract(R0, A1, A0, N2);
 
2141
                        P::Subtract(R1, B0, B1, N2);
 
2142
                        RecursiveMultiply<P>(T0, T2, R0, R1, N2);
2162
2143
                        carry = 0;
2163
2144
                        break;
2164
2145
                case 2:
2165
 
                        LowLevel::Subtract(R0, A0, A1, N2);
2166
 
                        LowLevel::Subtract(R1, B1, B0, N2);
2167
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
 
2146
                        P::Subtract(R0, A0, A1, N2);
 
2147
                        P::Subtract(R1, B1, B0, N2);
 
2148
                        RecursiveMultiply<P>(T0, T2, R0, R1, N2);
2168
2149
                        carry = 0;
2169
2150
                        break;
2170
2151
                case 4:
2171
 
                        LowLevel::Subtract(R0, A1, A0, N2);
2172
 
                        LowLevel::Subtract(R1, B0, B1, N2);
2173
 
                        RecursiveMultiply(T0, T2, R0, R1, N2);
2174
 
                        LowLevel::Subtract(T1, T1, R1, N2);
 
2152
                        P::Subtract(R0, A1, A0, N2);
 
2153
                        P::Subtract(R1, B0, B1, N2);
 
2154
                        RecursiveMultiply<P>(T0, T2, R0, R1, N2);
 
2155
                        P::Subtract(T1, T1, R1, N2);
2175
2156
                        carry = -1;
2176
2157
                        break;
2177
2158
                default:
2179
2160
                        carry = 0;
2180
2161
                }
2181
2162
 
2182
 
                RecursiveMultiply(T2, R0, A1, B1, N2);
 
2163
                RecursiveMultiply<P>(T2, R0, A1, B1, N2);
2183
2164
 
2184
2165
                // now T[01] holds (A1-A0)*(B0-B1), T[23] holds A1*B1
2185
2166
 
2186
 
                int c2 = LowLevel::Subtract(R0, L+N2, L, N2);
2187
 
                c2 += LowLevel::Subtract(R0, R0, T0, N2);
2188
 
                int t = (Compare(R0, T2, N2) == -1);
 
2167
                word c2 = P::Subtract(R0, L+N2, L, N2);
 
2168
                c2 += P::Subtract(R0, R0, T0, N2);
 
2169
                word t = (Compare(R0, T2, N2) == -1);
2189
2170
 
2190
2171
                carry += t;
2191
2172
                carry += Increment(R0, N2, c2+t);
2192
 
                carry += LowLevel::Add(R0, R0, T1, N2);
2193
 
                carry += LowLevel::Add(R0, R0, T3, N2);
 
2173
                carry += P::Add(R0, R0, T1, N2);
 
2174
                carry += P::Add(R0, R0, T3, N2);
2194
2175
                assert (carry >= 0 && carry <= 2);
2195
2176
 
2196
2177
                CopyWords(R1, T3, N2);
2198
2179
        }
2199
2180
}
2200
2181
 
2201
 
inline int Add(word *C, const word *A, const word *B, size_t N)
 
2182
inline word Add(word *C, const word *A, const word *B, unsigned int N)
2202
2183
{
2203
2184
        return LowLevel::Add(C, A, B, N);
2204
2185
}
2205
2186
 
2206
 
inline int Subtract(word *C, const word *A, const word *B, size_t N)
 
2187
inline word Subtract(word *C, const word *A, const word *B, unsigned int N)
2207
2188
{
2208
2189
        return LowLevel::Subtract(C, A, B, N);
2209
2190
}
2210
2191
 
2211
 
inline void Multiply(word *R, word *T, const word *A, const word *B, size_t N)
2212
 
{
2213
 
        RecursiveMultiply(R, T, A, B, N);
2214
 
}
2215
 
 
2216
 
inline void Square(word *R, word *T, const word *A, size_t N)
2217
 
{
2218
 
        RecursiveSquare(R, T, A, N);
2219
 
}
2220
 
 
2221
 
inline void MultiplyBottom(word *R, word *T, const word *A, const word *B, size_t N)
2222
 
{
2223
 
        RecursiveMultiplyBottom(R, T, A, B, N);
2224
 
}
2225
 
 
2226
 
inline void MultiplyTop(word *R, word *T, const word *L, const word *A, const word *B, size_t N)
2227
 
{
2228
 
        RecursiveMultiplyTop(R, T, L, A, B, N);
2229
 
}
2230
 
 
2231
 
static word LinearMultiply(word *C, const word *A, word B, size_t N)
 
2192
inline void Multiply(word *R, word *T, const word *A, const word *B, unsigned int N)
 
2193
{
 
2194
#ifdef SSE2_INTRINSICS_AVAILABLE
 
2195
        if (HasSSE2())
 
2196
                RecursiveMultiply<P4Optimized>(R, T, A, B, N);
 
2197
        else
 
2198
#endif
 
2199
                RecursiveMultiply<LowLevel>(R, T, A, B, N);
 
2200
}
 
2201
 
 
2202
inline void Square(word *R, word *T, const word *A, unsigned int N)
 
2203
{
 
2204
#ifdef SSE2_INTRINSICS_AVAILABLE
 
2205
        if (HasSSE2())
 
2206
                RecursiveSquare<P4Optimized>(R, T, A, N);
 
2207
        else
 
2208
#endif
 
2209
                RecursiveSquare<LowLevel>(R, T, A, N);
 
2210
}
 
2211
 
 
2212
inline void MultiplyBottom(word *R, word *T, const word *A, const word *B, unsigned int N)
 
2213
{
 
2214
#ifdef SSE2_INTRINSICS_AVAILABLE
 
2215
        if (HasSSE2())
 
2216
                RecursiveMultiplyBottom<P4Optimized>(R, T, A, B, N);
 
2217
        else
 
2218
#endif
 
2219
                RecursiveMultiplyBottom<LowLevel>(R, T, A, B, N);
 
2220
}
 
2221
 
 
2222
inline void MultiplyTop(word *R, word *T, const word *L, const word *A, const word *B, unsigned int N)
 
2223
{
 
2224
#ifdef SSE2_INTRINSICS_AVAILABLE
 
2225
        if (HasSSE2())
 
2226
                RecursiveMultiplyTop<P4Optimized>(R, T, L, A, B, N);
 
2227
        else
 
2228
#endif
 
2229
                RecursiveMultiplyTop<LowLevel>(R, T, L, A, B, N);
 
2230
}
 
2231
 
 
2232
static word LinearMultiply(word *C, const word *A, word B, unsigned int N)
2232
2233
{
2233
2234
        word carry=0;
2234
2235
        for(unsigned i=0; i<N; i++)
2245
2246
// A[NA] ---- multiplier
2246
2247
// B[NB] ---- multiplicant
2247
2248
 
2248
 
void AsymmetricMultiply(word *R, word *T, const word *A, size_t NA, const word *B, size_t NB)
 
2249
void AsymmetricMultiply(word *R, word *T, const word *A, unsigned int NA, const word *B, unsigned int NB)
2249
2250
{
2250
2251
        if (NA == NB)
2251
2252
        {
2287
2288
        Multiply(R, T, A, B, NA);
2288
2289
        CopyWords(T+2*NA, R+NA, NA);
2289
2290
 
2290
 
        size_t i;
 
2291
        unsigned i;
2291
2292
 
2292
2293
        for (i=2*NA; i<NB; i+=2*NA)
2293
2294
                Multiply(T+NA+i, T, A, B+i, NA);
2302
2303
// T[3*N/2] - temporary work space
2303
2304
// A[N] ----- an odd number as input
2304
2305
 
2305
 
void RecursiveInverseModPower2(word *R, word *T, const word *A, size_t N)
 
2306
void RecursiveInverseModPower2(word *R, word *T, const word *A, unsigned int N)
2306
2307
{
2307
2308
        if (N==2)
2308
2309
        {
2315
2316
        }
2316
2317
        else
2317
2318
        {
2318
 
                const size_t N2 = N/2;
 
2319
                const unsigned int N2 = N/2;
2319
2320
                RecursiveInverseModPower2(R0, T0, A0, N2);
2320
2321
                T0[0] = 1;
2321
2322
                SetWords(T0+1, 0, N2-1);
2333
2334
// M[N] --- modulus
2334
2335
// U[N] --- multiplicative inverse of M mod 2**(WORD_BITS*N)
2335
2336
 
2336
 
void MontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, size_t N)
 
2337
void MontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, unsigned int N)
2337
2338
{
2338
2339
        MultiplyBottom(R, T, X, U, N);
2339
2340
        MultiplyTop(T, T+N, X, R, M, N);
2351
2352
// U[N/2] - multiplicative inverse of M mod 2**(WORD_BITS*N/2)
2352
2353
// V[N] --- 2**(WORD_BITS*3*N/2) mod M
2353
2354
 
2354
 
void HalfMontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, const word *V, size_t N)
 
2355
void HalfMontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, const word *V, unsigned int N)
2355
2356
{
2356
2357
        assert(N%2==0 && N>=4);
2357
2358
 
2365
2366
#define X2              (X+N)
2366
2367
#define X3              (X+N+N2)
2367
2368
 
2368
 
        const size_t N2 = N/2;
 
2369
        const unsigned int N2 = N/2;
2369
2370
        Multiply(T0, T2, V0, X3, N2);
2370
2371
        int c2 = Add(T0, T0, X0, N);
2371
2372
        MultiplyBottom(T3, T2, T0, U, N2);
2499
2500
}
2500
2501
 
2501
2502
// for use by Divide(), corrects the underestimated quotient {Q1,Q0}
2502
 
static void CorrectQuotientEstimate(word *R, word *T, word *Q, const word *B, size_t N)
 
2503
static void CorrectQuotientEstimate(word *R, word *T, word *Q, const word *B, unsigned int N)
2503
2504
{
2504
2505
        assert(N && N%2==0);
2505
2506
 
2536
2537
// A[NA] -------- dividend
2537
2538
// B[NB] -------- divisor
2538
2539
 
2539
 
void Divide(word *R, word *Q, word *T, const word *A, size_t NA, const word *B, size_t NB)
 
2540
void Divide(word *R, word *Q, word *T, const word *A, unsigned int NA, const word *B, unsigned int NB)
2540
2541
{
2541
2542
        assert(NA && NB && NA%2==0 && NB%2==0);
2542
2543
        assert(B[NB-1] || B[NB-2]);
2580
2581
        BT[1] = TB[NB-1] + (BT[0]==0);
2581
2582
 
2582
2583
        // start reducing TA mod TB, 2 words at a time
2583
 
        for (size_t i=NA-2; i>=NB; i-=2)
 
2584
        for (unsigned i=NA-2; i>=NB; i-=2)
2584
2585
        {
2585
2586
                AtomicDivide(Q+i-NB, TA+i-2, BT);
2586
2587
                CorrectQuotientEstimate(TA+i-NB, TP, Q+i-NB, TB, NB);
2591
2592
        ShiftWordsRightByBits(R, NB, shiftBits);
2592
2593
}
2593
2594
 
2594
 
static inline size_t EvenWordCount(const word *X, size_t N)
 
2595
static inline unsigned int EvenWordCount(const word *X, unsigned int N)
2595
2596
{
2596
2597
        while (N && X[N-2]==0 && X[N-1]==0)
2597
2598
                N-=2;
2604
2605
// A[NA] -- number to take inverse of
2605
2606
// M[N] --- modulus
2606
2607
 
2607
 
unsigned int AlmostInverse(word *R, word *T, const word *A, size_t NA, const word *M, size_t N)
 
2608
unsigned int AlmostInverse(word *R, word *T, const word *A, unsigned int NA, const word *M, unsigned int N)
2608
2609
{
2609
2610
        assert(NA<=N && N && N%2==0);
2610
2611
 
2612
2613
        word *c = T+N;
2613
2614
        word *f = T+2*N;
2614
2615
        word *g = T+3*N;
2615
 
        size_t bcLen=2, fgLen=EvenWordCount(M, N);
 
2616
        unsigned int bcLen=2, fgLen=EvenWordCount(M, N);
2616
2617
        unsigned int k=0, s=0;
2617
2618
 
2618
2619
        SetWords(T, 0, 3*N);
2690
2691
// A[N] - input
2691
2692
// M[N] - modulus
2692
2693
 
2693
 
void DivideByPower2Mod(word *R, const word *A, size_t k, const word *M, size_t N)
 
2694
void DivideByPower2Mod(word *R, const word *A, unsigned int k, const word *M, unsigned int N)
2694
2695
{
2695
2696
        CopyWords(R, A, N);
2696
2697
 
2711
2712
// A[N] - input
2712
2713
// M[N] - modulus
2713
2714
 
2714
 
void MultiplyByPower2Mod(word *R, const word *A, size_t k, const word *M, size_t N)
 
2715
void MultiplyByPower2Mod(word *R, const word *A, unsigned int k, const word *M, unsigned int N)
2715
2716
{
2716
2717
        CopyWords(R, A, N);
2717
2718
 
2722
2723
 
2723
2724
// ******************************************************************
2724
2725
 
2725
 
InitializeInteger::InitializeInteger()
2726
 
{
2727
 
        if (!g_pAssignIntToInteger)
2728
 
        {
2729
 
#ifdef CRYPTOPP_X86ASM_AVAILABLE
2730
 
                SetPentiumFunctionPointers();
2731
 
#endif
2732
 
                g_pAssignIntToInteger = AssignIntToInteger;
2733
 
        }
2734
 
}
2735
 
 
2736
2726
static const unsigned int RoundupSizeTable[] = {2, 2, 2, 4, 4, 8, 8, 8, 8};
2737
2727
 
2738
 
static inline size_t RoundupSize(size_t n)
 
2728
static inline unsigned int RoundupSize(unsigned int n)
2739
2729
{
2740
2730
        if (n<=8)
2741
2731
                return RoundupSizeTable[n];
2745
2735
                return 32;
2746
2736
        else if (n<=64)
2747
2737
                return 64;
2748
 
        else return size_t(1) << BitPrecision(n-1);
 
2738
        else return 1U << BitPrecision(n-1);
2749
2739
}
2750
2740
 
2751
2741
Integer::Integer()
2793
2783
        if (ByteCount() > sizeof(long))
2794
2784
                return false;
2795
2785
 
2796
 
        unsigned long value = (unsigned long)reg[0];
2797
 
        value += SafeLeftShift<WORD_BITS, unsigned long>((unsigned long)reg[1]);
 
2786
        unsigned long value = reg[0];
 
2787
        value += SafeLeftShift<WORD_BITS, unsigned long>(reg[1]);
2798
2788
 
2799
2789
        if (sign==POSITIVE)
2800
2790
                return (signed long)value >= 0;
2806
2796
{
2807
2797
        assert(IsConvertableToLong());
2808
2798
 
2809
 
        unsigned long value = (unsigned long)reg[0];
2810
 
        value += SafeLeftShift<WORD_BITS, unsigned long>((unsigned long)reg[1]);
 
2799
        unsigned long value = reg[0];
 
2800
        value += SafeLeftShift<WORD_BITS, unsigned long>(reg[1]);
2811
2801
        return sign==POSITIVE ? value : -(signed long)value;
2812
2802
}
2813
2803
 
2814
 
Integer::Integer(BufferedTransformation &encodedInteger, size_t byteCount, Signedness s)
 
2804
Integer::Integer(BufferedTransformation &encodedInteger, unsigned int byteCount, Signedness s)
2815
2805
{
2816
2806
        Decode(encodedInteger, byteCount, s);
2817
2807
}
2818
2808
 
2819
 
Integer::Integer(const byte *encodedInteger, size_t byteCount, Signedness s)
 
2809
Integer::Integer(const byte *encodedInteger, unsigned int byteCount, Signedness s)
2820
2810
{
2821
2811
        Decode(encodedInteger, byteCount, s);
2822
2812
}
2826
2816
        BERDecode(bt);
2827
2817
}
2828
2818
 
2829
 
Integer::Integer(RandomNumberGenerator &rng, size_t bitcount)
 
2819
Integer::Integer(RandomNumberGenerator &rng, unsigned int bitcount)
2830
2820
{
2831
2821
        Randomize(rng, bitcount);
2832
2822
}
2837
2827
                throw Integer::RandomNumberNotFound();
2838
2828
}
2839
2829
 
2840
 
Integer Integer::Power2(size_t e)
 
2830
Integer Integer::Power2(unsigned int e)
2841
2831
{
2842
2832
        Integer r((word)0, BitsToWords(e+1));
2843
2833
        r.SetBit(e);
2884
2874
        return *this;
2885
2875
}
2886
2876
 
2887
 
bool Integer::GetBit(size_t n) const
 
2877
bool Integer::GetBit(unsigned int n) const
2888
2878
{
2889
2879
        if (n/WORD_BITS >= reg.size())
2890
2880
                return 0;
2892
2882
                return bool((reg[n/WORD_BITS] >> (n % WORD_BITS)) & 1);
2893
2883
}
2894
2884
 
2895
 
void Integer::SetBit(size_t n, bool value)
 
2885
void Integer::SetBit(unsigned int n, bool value)
2896
2886
{
2897
2887
        if (value)
2898
2888
        {
2906
2896
        }
2907
2897
}
2908
2898
 
2909
 
byte Integer::GetByte(size_t n) const
 
2899
byte Integer::GetByte(unsigned int n) const
2910
2900
{
2911
2901
        if (n/WORD_SIZE >= reg.size())
2912
2902
                return 0;
2914
2904
                return byte(reg[n/WORD_SIZE] >> ((n%WORD_SIZE)*8));
2915
2905
}
2916
2906
 
2917
 
void Integer::SetByte(size_t n, byte value)
 
2907
void Integer::SetByte(unsigned int n, byte value)
2918
2908
{
2919
2909
        reg.CleanGrow(RoundupSize(BytesToWords(n+1)));
2920
2910
        reg[n/WORD_SIZE] &= ~(word(0xff) << 8*(n%WORD_SIZE));
2921
2911
        reg[n/WORD_SIZE] |= (word(value) << 8*(n%WORD_SIZE));
2922
2912
}
2923
2913
 
2924
 
lword Integer::GetBits(size_t i, size_t n) const
 
2914
unsigned long Integer::GetBits(unsigned int i, unsigned int n) const
2925
2915
{
2926
 
        lword v = 0;
2927
 
        assert(n <= sizeof(v)*8);
 
2916
        assert(n <= sizeof(unsigned long)*8);
 
2917
        unsigned long v = 0;
2928
2918
        for (unsigned int j=0; j<n; j++)
2929
 
                v |= lword(GetBit(i+j)) << j;
 
2919
                v |= GetBit(i+j) << j;
2930
2920
        return v;
2931
2921
}
2932
2922
 
2950
2940
        std::swap(sign, a.sign);
2951
2941
}
2952
2942
 
2953
 
Integer::Integer(word value, size_t length)
 
2943
Integer::Integer(word value, unsigned int length)
2954
2944
        : reg(RoundupSize(length)), sign(POSITIVE)
2955
2945
{
2956
2946
        reg[0] = value;
2960
2950
template <class T>
2961
2951
static Integer StringToInteger(const T *str)
2962
2952
{
2963
 
        int radix;
 
2953
        word radix;
2964
2954
        // GCC workaround
 
2955
        // std::char_traits doesn't exist in GCC 2.x
2965
2956
        // std::char_traits<wchar_t>::length() not defined in GCC 3.2 and STLport 4.5.3
2966
2957
        unsigned int length;
2967
2958
        for (length = 0; str[length] != 0; length++) {}
2994
2985
 
2995
2986
        for (unsigned i=0; i<length; i++)
2996
2987
        {
2997
 
                int digit;
 
2988
                word digit;
2998
2989
 
2999
2990
                if (str[i] >= '0' && str[i] <= '9')
3000
2991
                        digit = str[i] - '0';
3032
3023
 
3033
3024
unsigned int Integer::WordCount() const
3034
3025
{
3035
 
        return (unsigned int)CountWords(reg, reg.size());
 
3026
        return CountWords(reg, reg.size());
3036
3027
}
3037
3028
 
3038
3029
unsigned int Integer::ByteCount() const
3053
3044
                return 0;
3054
3045
}
3055
3046
 
3056
 
void Integer::Decode(const byte *input, size_t inputLen, Signedness s)
 
3047
void Integer::Decode(const byte *input, unsigned int inputLen, Signedness s)
3057
3048
{
3058
3049
        StringStore store(input, inputLen);
3059
3050
        Decode(store, inputLen, s);
3060
3051
}
3061
3052
 
3062
 
void Integer::Decode(BufferedTransformation &bt, size_t inputLen, Signedness s)
 
3053
void Integer::Decode(BufferedTransformation &bt, unsigned int inputLen, Signedness s)
3063
3054
{
3064
3055
        assert(bt.MaxRetrievable() >= inputLen);
3065
3056
 
3076
3067
 
3077
3068
        reg.CleanNew(RoundupSize(BytesToWords(inputLen)));
3078
3069
 
3079
 
        for (size_t i=inputLen; i > 0; i--)
 
3070
        for (unsigned int i=inputLen; i > 0; i--)
3080
3071
        {
3081
3072
                bt.Get(b);
3082
3073
                reg[(i-1)/WORD_SIZE] |= word(b) << ((i-1)%WORD_SIZE)*8;
3084
3075
 
3085
3076
        if (sign == NEGATIVE)
3086
3077
        {
3087
 
                for (size_t i=inputLen; i<reg.size()*WORD_SIZE; i++)
 
3078
                for (unsigned i=inputLen; i<reg.size()*WORD_SIZE; i++)
3088
3079
                        reg[i/WORD_SIZE] |= word(0xff) << (i%WORD_SIZE)*8;
3089
3080
                TwosComplement(reg, reg.size());
3090
3081
        }
3091
3082
}
3092
3083
 
3093
 
size_t Integer::MinEncodedSize(Signedness signedness) const
 
3084
unsigned int Integer::MinEncodedSize(Signedness signedness) const
3094
3085
{
3095
3086
        unsigned int outputLen = STDMAX(1U, ByteCount());
3096
3087
        if (signedness == UNSIGNED)
3102
3093
        return outputLen;
3103
3094
}
3104
3095
 
3105
 
void Integer::Encode(byte *output, size_t outputLen, Signedness signedness) const
 
3096
unsigned int Integer::Encode(byte *output, unsigned int outputLen, Signedness signedness) const
3106
3097
{
3107
3098
        ArraySink sink(output, outputLen);
3108
 
        Encode(sink, outputLen, signedness);
 
3099
        return Encode(sink, outputLen, signedness);
3109
3100
}
3110
3101
 
3111
 
void Integer::Encode(BufferedTransformation &bt, size_t outputLen, Signedness signedness) const
 
3102
unsigned int Integer::Encode(BufferedTransformation &bt, unsigned int outputLen, Signedness signedness) const
3112
3103
{
3113
3104
        if (signedness == UNSIGNED || NotNegative())
3114
3105
        {
3115
 
                for (size_t i=outputLen; i > 0; i--)
 
3106
                for (unsigned int i=outputLen; i > 0; i--)
3116
3107
                        bt.Put(GetByte(i-1));
3117
3108
        }
3118
3109
        else
3119
3110
        {
3120
3111
                // take two's complement of *this
3121
 
                Integer temp = Integer::Power2(8*STDMAX((size_t)ByteCount(), outputLen)) + *this;
3122
 
                temp.Encode(bt, outputLen, UNSIGNED);
 
3112
                Integer temp = Integer::Power2(8*STDMAX(ByteCount(), outputLen)) + *this;
 
3113
                for (unsigned i=0; i<outputLen; i++)
 
3114
                        bt.Put(temp.GetByte(outputLen-i-1));
3123
3115
        }
 
3116
        return outputLen;
3124
3117
}
3125
3118
 
3126
3119
void Integer::DEREncode(BufferedTransformation &bt) const
3130
3123
        enc.MessageEnd();
3131
3124
}
3132
3125
 
3133
 
void Integer::BERDecode(const byte *input, size_t len)
 
3126
void Integer::BERDecode(const byte *input, unsigned int len)
3134
3127
{
3135
3128
        StringStore store(input, len);
3136
3129
        BERDecode(store);
3141
3134
        BERGeneralDecoder dec(bt, INTEGER);
3142
3135
        if (!dec.IsDefiniteLength() || dec.MaxRetrievable() < dec.RemainingLength())
3143
3136
                BERDecodeError();
3144
 
        Decode(dec, (size_t)dec.RemainingLength(), SIGNED);
 
3137
        Decode(dec, dec.RemainingLength(), SIGNED);
3145
3138
        dec.MessageEnd();
3146
3139
}
3147
3140
 
3148
 
void Integer::DEREncodeAsOctetString(BufferedTransformation &bt, size_t length) const
 
3141
void Integer::DEREncodeAsOctetString(BufferedTransformation &bt, unsigned int length) const
3149
3142
{
3150
3143
        DERGeneralEncoder enc(bt, OCTET_STRING);
3151
3144
        Encode(enc, length);
3152
3145
        enc.MessageEnd();
3153
3146
}
3154
3147
 
3155
 
void Integer::BERDecodeAsOctetString(BufferedTransformation &bt, size_t length)
 
3148
void Integer::BERDecodeAsOctetString(BufferedTransformation &bt, unsigned int length)
3156
3149
{
3157
3150
        BERGeneralDecoder dec(bt, OCTET_STRING);
3158
3151
        if (!dec.IsDefiniteLength() || dec.RemainingLength() != length)
3161
3154
        dec.MessageEnd();
3162
3155
}
3163
3156
 
3164
 
size_t Integer::OpenPGPEncode(byte *output, size_t len) const
 
3157
unsigned int Integer::OpenPGPEncode(byte *output, unsigned int len) const
3165
3158
{
3166
3159
        ArraySink sink(output, len);
3167
3160
        return OpenPGPEncode(sink);
3168
3161
}
3169
3162
 
3170
 
size_t Integer::OpenPGPEncode(BufferedTransformation &bt) const
 
3163
unsigned int Integer::OpenPGPEncode(BufferedTransformation &bt) const
3171
3164
{
3172
3165
        word16 bitCount = BitCount();
3173
3166
        bt.PutWord16(bitCount);
3174
 
        size_t byteCount = BitsToBytes(bitCount);
3175
 
        Encode(bt, byteCount);
3176
 
        return 2 + byteCount;
 
3167
        return 2 + Encode(bt, BitsToBytes(bitCount));
3177
3168
}
3178
3169
 
3179
 
void Integer::OpenPGPDecode(const byte *input, size_t len)
 
3170
void Integer::OpenPGPDecode(const byte *input, unsigned int len)
3180
3171
{
3181
3172
        StringStore store(input, len);
3182
3173
        OpenPGPDecode(store);
3190
3181
        Decode(bt, BitsToBytes(bitCount));
3191
3182
}
3192
3183
 
3193
 
void Integer::Randomize(RandomNumberGenerator &rng, size_t nbits)
 
3184
void Integer::Randomize(RandomNumberGenerator &rng, unsigned int nbits)
3194
3185
{
3195
 
        const size_t nbytes = nbits/8 + 1;
 
3186
        const unsigned int nbytes = nbits/8 + 1;
3196
3187
        SecByteBlock buf(nbytes);
3197
3188
        rng.GenerateBlock(buf, nbytes);
3198
3189
        if (nbytes)
3225
3216
class KDF2_RNG : public RandomNumberGenerator
3226
3217
{
3227
3218
public:
3228
 
        KDF2_RNG(const byte *seed, size_t seedSize)
 
3219
        KDF2_RNG(const byte *seed, unsigned int seedSize)
3229
3220
                : m_counter(0), m_counterAndSeed(seedSize + 4)
3230
3221
        {
3231
3222
                memcpy(m_counterAndSeed + 4, seed, seedSize);
3287
3278
                DEREncodeOctetString(seq, seed.begin(), seed.size());
3288
3279
                seq.MessageEnd();
3289
3280
 
3290
 
                SecByteBlock finalSeed((size_t)bq.MaxRetrievable());
 
3281
                SecByteBlock finalSeed(bq.MaxRetrievable());
3291
3282
                bq.Get(finalSeed, finalSeed.size());
3292
3283
                kdf2Rng.reset(new KDF2_RNG(finalSeed.begin(), finalSeed.size()));
3293
3284
        }
3463
3454
 
3464
3455
void PositiveAdd(Integer &sum, const Integer &a, const Integer& b)
3465
3456
{
3466
 
        int carry;
 
3457
        word carry;
3467
3458
        if (a.reg.size() == b.reg.size())
3468
3459
                carry = Add(sum.reg, a.reg, b.reg, a.reg.size());
3469
3460
        else if (a.reg.size() > b.reg.size())
3525
3516
        }
3526
3517
}
3527
3518
 
3528
 
// MSVC .NET 2003 workaround
3529
 
template <class T> inline const T& STDMAX2(const T& a, const T& b)
3530
 
{
3531
 
        return a < b ? b : a;
3532
 
}
3533
 
 
3534
3519
Integer Integer::Plus(const Integer& b) const
3535
3520
{
3536
 
        Integer sum((word)0, STDMAX2(reg.size(), b.reg.size()));
 
3521
        Integer sum((word)0, STDMAX(reg.size(), b.reg.size()));
3537
3522
        if (NotNegative())
3538
3523
        {
3539
3524
                if (b.NotNegative())
3579
3564
 
3580
3565
Integer Integer::Minus(const Integer& b) const
3581
3566
{
3582
 
        Integer diff((word)0, STDMAX2(reg.size(), b.reg.size()));
 
3567
        Integer diff((word)0, STDMAX(reg.size(), b.reg.size()));
3583
3568
        if (NotNegative())
3584
3569
        {
3585
3570
                if (b.NotNegative())
3623
3608
        return *this;
3624
3609
}
3625
3610
 
3626
 
Integer& Integer::operator<<=(size_t n)
 
3611
Integer& Integer::operator<<=(unsigned int n)
3627
3612
{
3628
 
        const size_t wordCount = WordCount();
3629
 
        const size_t shiftWords = n / WORD_BITS;
3630
 
        const unsigned int shiftBits = (unsigned int)(n % WORD_BITS);
 
3613
        const unsigned int wordCount = WordCount();
 
3614
        const unsigned int shiftWords = n / WORD_BITS;
 
3615
        const unsigned int shiftBits = n % WORD_BITS;
3631
3616
 
3632
3617
        reg.CleanGrow(RoundupSize(wordCount+BitsToWords(n)));
3633
3618
        ShiftWordsLeftByWords(reg, wordCount + shiftWords, shiftWords);
3635
3620
        return *this;
3636
3621
}
3637
3622
 
3638
 
Integer& Integer::operator>>=(size_t n)
 
3623
Integer& Integer::operator>>=(unsigned int n)
3639
3624
{
3640
 
        const size_t wordCount = WordCount();
3641
 
        const size_t shiftWords = n / WORD_BITS;
3642
 
        const unsigned int shiftBits = (unsigned int)(n % WORD_BITS);
 
3625
        const unsigned int wordCount = WordCount();
 
3626
        const unsigned int shiftWords = n / WORD_BITS;
 
3627
        const unsigned int shiftBits = n % WORD_BITS;
3643
3628
 
3644
3629
        ShiftWordsRightByWords(reg, wordCount, shiftWords);
3645
3630
        if (wordCount > shiftWords)
3651
3636
 
3652
3637
void PositiveMultiply(Integer &product, const Integer &a, const Integer &b)
3653
3638
{
3654
 
        size_t aSize = RoundupSize(a.WordCount());
3655
 
        size_t bSize = RoundupSize(b.WordCount());
 
3639
        unsigned aSize = RoundupSize(a.WordCount());
 
3640
        unsigned bSize = RoundupSize(b.WordCount());
3656
3641
 
3657
3642
        product.reg.CleanNew(RoundupSize(aSize+bSize));
3658
3643
        product.sign = Integer::POSITIVE;
3750
3735
        q = a;
3751
3736
        q >>= n;
3752
3737
 
3753
 
        const size_t wordCount = BitsToWords(n);
 
3738
        const unsigned int wordCount = BitsToWords(n);
3754
3739
        if (wordCount <= a.WordCount())
3755
3740
        {
3756
3741
                r.reg.resize(RoundupSize(wordCount));
3757
3742
                CopyWords(r.reg, a.reg, wordCount);
3758
3743
                SetWords(r.reg+wordCount, 0, r.reg.size()-wordCount);
3759
3744
                if (n % WORD_BITS != 0)
3760
 
                        r.reg[wordCount-1] %= (word(1) << (n % WORD_BITS));
 
3745
                        r.reg[wordCount-1] %= (1 << (n % WORD_BITS));
3761
3746
        }
3762
3747
        else
3763
3748
        {
4009
3994
        OID oid(seq);
4010
3995
        if (oid != ASN1::prime_field())
4011
3996
                BERDecodeError();
4012
 
        m_modulus.BERDecode(seq);
 
3997
        modulus.BERDecode(seq);
4013
3998
        seq.MessageEnd();
4014
 
        m_result.reg.resize(m_modulus.reg.size());
 
3999
        result.reg.resize(modulus.reg.size());
4015
4000
}
4016
4001
 
4017
4002
void ModularArithmetic::DEREncode(BufferedTransformation &bt) const
4018
4003
{
4019
4004
        DERSequenceEncoder seq(bt);
4020
4005
        ASN1::prime_field().DEREncode(seq);
4021
 
        m_modulus.DEREncode(seq);
 
4006
        modulus.DEREncode(seq);
4022
4007
        seq.MessageEnd();
4023
4008
}
4024
4009
 
4034
4019
 
4035
4020
const Integer& ModularArithmetic::Half(const Integer &a) const
4036
4021
{
4037
 
        if (a.reg.size()==m_modulus.reg.size())
 
4022
        if (a.reg.size()==modulus.reg.size())
4038
4023
        {
4039
 
                CryptoPP::DivideByPower2Mod(m_result.reg.begin(), a.reg, 1, m_modulus.reg, a.reg.size());
4040
 
                return m_result;
 
4024
                CryptoPP::DivideByPower2Mod(result.reg.begin(), a.reg, 1, modulus.reg, a.reg.size());
 
4025
                return result;
4041
4026
        }
4042
4027
        else
4043
 
                return m_result1 = (a.IsEven() ? (a >> 1) : ((a+m_modulus) >> 1));
 
4028
                return result1 = (a.IsEven() ? (a >> 1) : ((a+modulus) >> 1));
4044
4029
}
4045
4030
 
4046
4031
const Integer& ModularArithmetic::Add(const Integer &a, const Integer &b) const
4047
4032
{
4048
 
        if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
 
4033
        if (a.reg.size()==modulus.reg.size() && b.reg.size()==modulus.reg.size())
4049
4034
        {
4050
 
                if (CryptoPP::Add(m_result.reg.begin(), a.reg, b.reg, a.reg.size())
4051
 
                        || Compare(m_result.reg, m_modulus.reg, a.reg.size()) >= 0)
 
4035
                if (CryptoPP::Add(result.reg.begin(), a.reg, b.reg, a.reg.size())
 
4036
                        || Compare(result.reg, modulus.reg, a.reg.size()) >= 0)
4052
4037
                {
4053
 
                        CryptoPP::Subtract(m_result.reg.begin(), m_result.reg, m_modulus.reg, a.reg.size());
 
4038
                        CryptoPP::Subtract(result.reg.begin(), result.reg, modulus.reg, a.reg.size());
4054
4039
                }
4055
 
                return m_result;
 
4040
                return result;
4056
4041
        }
4057
4042
        else
4058
4043
        {
4059
 
                m_result1 = a+b;
4060
 
                if (m_result1 >= m_modulus)
4061
 
                        m_result1 -= m_modulus;
4062
 
                return m_result1;
 
4044
                result1 = a+b;
 
4045
                if (result1 >= modulus)
 
4046
                        result1 -= modulus;
 
4047
                return result1;
4063
4048
        }
4064
4049
}
4065
4050
 
4066
4051
Integer& ModularArithmetic::Accumulate(Integer &a, const Integer &b) const
4067
4052
{
4068
 
        if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
 
4053
        if (a.reg.size()==modulus.reg.size() && b.reg.size()==modulus.reg.size())
4069
4054
        {
4070
4055
                if (CryptoPP::Add(a.reg, a.reg, b.reg, a.reg.size())
4071
 
                        || Compare(a.reg, m_modulus.reg, a.reg.size()) >= 0)
 
4056
                        || Compare(a.reg, modulus.reg, a.reg.size()) >= 0)
4072
4057
                {
4073
 
                        CryptoPP::Subtract(a.reg, a.reg, m_modulus.reg, a.reg.size());
 
4058
                        CryptoPP::Subtract(a.reg, a.reg, modulus.reg, a.reg.size());
4074
4059
                }
4075
4060
        }
4076
4061
        else
4077
4062
        {
4078
4063
                a+=b;
4079
 
                if (a>=m_modulus)
4080
 
                        a-=m_modulus;
 
4064
                if (a>=modulus)
 
4065
                        a-=modulus;
4081
4066
        }
4082
4067
 
4083
4068
        return a;
4085
4070
 
4086
4071
const Integer& ModularArithmetic::Subtract(const Integer &a, const Integer &b) const
4087
4072
{
4088
 
        if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
 
4073
        if (a.reg.size()==modulus.reg.size() && b.reg.size()==modulus.reg.size())
4089
4074
        {
4090
 
                if (CryptoPP::Subtract(m_result.reg.begin(), a.reg, b.reg, a.reg.size()))
4091
 
                        CryptoPP::Add(m_result.reg.begin(), m_result.reg, m_modulus.reg, a.reg.size());
4092
 
                return m_result;
 
4075
                if (CryptoPP::Subtract(result.reg.begin(), a.reg, b.reg, a.reg.size()))
 
4076
                        CryptoPP::Add(result.reg.begin(), result.reg, modulus.reg, a.reg.size());
 
4077
                return result;
4093
4078
        }
4094
4079
        else
4095
4080
        {
4096
 
                m_result1 = a-b;
4097
 
                if (m_result1.IsNegative())
4098
 
                        m_result1 += m_modulus;
4099
 
                return m_result1;
 
4081
                result1 = a-b;
 
4082
                if (result1.IsNegative())
 
4083
                        result1 += modulus;
 
4084
                return result1;
4100
4085
        }
4101
4086
}
4102
4087
 
4103
4088
Integer& ModularArithmetic::Reduce(Integer &a, const Integer &b) const
4104
4089
{
4105
 
        if (a.reg.size()==m_modulus.reg.size() && b.reg.size()==m_modulus.reg.size())
 
4090
        if (a.reg.size()==modulus.reg.size() && b.reg.size()==modulus.reg.size())
4106
4091
        {
4107
4092
                if (CryptoPP::Subtract(a.reg, a.reg, b.reg, a.reg.size()))
4108
 
                        CryptoPP::Add(a.reg, a.reg, m_modulus.reg, a.reg.size());
 
4093
                        CryptoPP::Add(a.reg, a.reg, modulus.reg, a.reg.size());
4109
4094
        }
4110
4095
        else
4111
4096
        {
4112
4097
                a-=b;
4113
4098
                if (a.IsNegative())
4114
 
                        a+=m_modulus;
 
4099
                        a+=modulus;
4115
4100
        }
4116
4101
 
4117
4102
        return a;
4122
4107
        if (!a)
4123
4108
                return a;
4124
4109
 
4125
 
        CopyWords(m_result.reg.begin(), m_modulus.reg, m_modulus.reg.size());
4126
 
        if (CryptoPP::Subtract(m_result.reg.begin(), m_result.reg, a.reg, a.reg.size()))
4127
 
                Decrement(m_result.reg.begin()+a.reg.size(), m_modulus.reg.size()-a.reg.size());
 
4110
        CopyWords(result.reg.begin(), modulus.reg, modulus.reg.size());
 
4111
        if (CryptoPP::Subtract(result.reg.begin(), result.reg, a.reg, a.reg.size()))
 
4112
                Decrement(result.reg.begin()+a.reg.size(), 1, modulus.reg.size()-a.reg.size());
4128
4113
 
4129
 
        return m_result;
 
4114
        return result;
4130
4115
}
4131
4116
 
4132
4117
Integer ModularArithmetic::CascadeExponentiate(const Integer &x, const Integer &e1, const Integer &y, const Integer &e2) const
4133
4118
{
4134
 
        if (m_modulus.IsOdd())
 
4119
        if (modulus.IsOdd())
4135
4120
        {
4136
 
                MontgomeryRepresentation dr(m_modulus);
 
4121
                MontgomeryRepresentation dr(modulus);
4137
4122
                return dr.ConvertOut(dr.CascadeExponentiate(dr.ConvertIn(x), e1, dr.ConvertIn(y), e2));
4138
4123
        }
4139
4124
        else
4142
4127
 
4143
4128
void ModularArithmetic::SimultaneousExponentiate(Integer *results, const Integer &base, const Integer *exponents, unsigned int exponentsCount) const
4144
4129
{
4145
 
        if (m_modulus.IsOdd())
 
4130
        if (modulus.IsOdd())
4146
4131
        {
4147
 
                MontgomeryRepresentation dr(m_modulus);
 
4132
                MontgomeryRepresentation dr(modulus);
4148
4133
                dr.SimultaneousExponentiate(results, dr.ConvertIn(base), exponents, exponentsCount);
4149
4134
                for (unsigned int i=0; i<exponentsCount; i++)
4150
4135
                        results[i] = dr.ConvertOut(results[i]);
4155
4140
 
4156
4141
MontgomeryRepresentation::MontgomeryRepresentation(const Integer &m)    // modulus must be odd
4157
4142
        : ModularArithmetic(m),
4158
 
          m_u((word)0, m_modulus.reg.size()),
4159
 
          m_workspace(5*m_modulus.reg.size())
 
4143
          u((word)0, modulus.reg.size()),
 
4144
          workspace(5*modulus.reg.size())
4160
4145
{
4161
 
        if (!m_modulus.IsOdd())
 
4146
        if (!modulus.IsOdd())
4162
4147
                throw InvalidArgument("MontgomeryRepresentation: Montgomery representation requires an odd modulus");
4163
4148
 
4164
 
        RecursiveInverseModPower2(m_u.reg, m_workspace, m_modulus.reg, m_modulus.reg.size());
 
4149
        RecursiveInverseModPower2(u.reg, workspace, modulus.reg, modulus.reg.size());
4165
4150
}
4166
4151
 
4167
4152
const Integer& MontgomeryRepresentation::Multiply(const Integer &a, const Integer &b) const
4168
4153
{
4169
 
        word *const T = m_workspace.begin();
4170
 
        word *const R = m_result.reg.begin();
4171
 
        const size_t N = m_modulus.reg.size();
 
4154
        word *const T = workspace.begin();
 
4155
        word *const R = result.reg.begin();
 
4156
        const unsigned int N = modulus.reg.size();
4172
4157
        assert(a.reg.size()<=N && b.reg.size()<=N);
4173
4158
 
4174
4159
        AsymmetricMultiply(T, T+2*N, a.reg, a.reg.size(), b.reg, b.reg.size());
4175
4160
        SetWords(T+a.reg.size()+b.reg.size(), 0, 2*N-a.reg.size()-b.reg.size());
4176
 
        MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
4177
 
        return m_result;
 
4161
        MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
 
4162
        return result;
4178
4163
}
4179
4164
 
4180
4165
const Integer& MontgomeryRepresentation::Square(const Integer &a) const
4181
4166
{
4182
 
        word *const T = m_workspace.begin();
4183
 
        word *const R = m_result.reg.begin();
4184
 
        const size_t N = m_modulus.reg.size();
 
4167
        word *const T = workspace.begin();
 
4168
        word *const R = result.reg.begin();
 
4169
        const unsigned int N = modulus.reg.size();
4185
4170
        assert(a.reg.size()<=N);
4186
4171
 
4187
4172
        CryptoPP::Square(T, T+2*N, a.reg, a.reg.size());
4188
4173
        SetWords(T+2*a.reg.size(), 0, 2*N-2*a.reg.size());
4189
 
        MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
4190
 
        return m_result;
 
4174
        MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
 
4175
        return result;
4191
4176
}
4192
4177
 
4193
4178
Integer MontgomeryRepresentation::ConvertOut(const Integer &a) const
4194
4179
{
4195
 
        word *const T = m_workspace.begin();
4196
 
        word *const R = m_result.reg.begin();
4197
 
        const size_t N = m_modulus.reg.size();
 
4180
        word *const T = workspace.begin();
 
4181
        word *const R = result.reg.begin();
 
4182
        const unsigned int N = modulus.reg.size();
4198
4183
        assert(a.reg.size()<=N);
4199
4184
 
4200
4185
        CopyWords(T, a.reg, a.reg.size());
4201
4186
        SetWords(T+a.reg.size(), 0, 2*N-a.reg.size());
4202
 
        MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
4203
 
        return m_result;
 
4187
        MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
 
4188
        return result;
4204
4189
}
4205
4190
 
4206
4191
const Integer& MontgomeryRepresentation::MultiplicativeInverse(const Integer &a) const
4207
4192
{
4208
4193
//        return (EuclideanMultiplicativeInverse(a, modulus)<<(2*WORD_BITS*modulus.reg.size()))%modulus;
4209
 
        word *const T = m_workspace.begin();
4210
 
        word *const R = m_result.reg.begin();
4211
 
        const size_t N = m_modulus.reg.size();
 
4194
        word *const T = workspace.begin();
 
4195
        word *const R = result.reg.begin();
 
4196
        const unsigned int N = modulus.reg.size();
4212
4197
        assert(a.reg.size()<=N);
4213
4198
 
4214
4199
        CopyWords(T, a.reg, a.reg.size());
4215
4200
        SetWords(T+a.reg.size(), 0, 2*N-a.reg.size());
4216
 
        MontgomeryReduce(R, T+2*N, T, m_modulus.reg, m_u.reg, N);
4217
 
        unsigned k = AlmostInverse(R, T, R, N, m_modulus.reg, N);
 
4201
        MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
 
4202
        unsigned k = AlmostInverse(R, T, R, N, modulus.reg, N);
4218
4203
 
4219
4204
//      cout << "k=" << k << " N*32=" << 32*N << endl;
4220
4205
 
4221
4206
        if (k>N*WORD_BITS)
4222
 
                DivideByPower2Mod(R, R, k-N*WORD_BITS, m_modulus.reg, N);
 
4207
                DivideByPower2Mod(R, R, k-N*WORD_BITS, modulus.reg, N);
4223
4208
        else
4224
 
                MultiplyByPower2Mod(R, R, N*WORD_BITS-k, m_modulus.reg, N);
 
4209
                MultiplyByPower2Mod(R, R, N*WORD_BITS-k, modulus.reg, N);
4225
4210
 
4226
 
        return m_result;
 
4211
        return result;
4227
4212
}
4228
4213
 
4229
4214
NAMESPACE_END