These serve only as vehicles of verification of the correctness of the algorithms, and test beds for experimenting. In a loop millions of function calls can be performed with random parameters, and compared to the result of the built-in extended GCD function, taking only a few minutes time.
for i = 1:1e6
x = floor(2^51*rand(1,2)); % 2^9 to see problems at small numbers
x(2) = max(bitor(x(2),1),3); % m = odd
b = XXXinv(x(1),x(2));
if modinv(x(1),x(2)) ~= b || b > x(2) || b < 0
disp([x,b])
end
end
Note that the parameters are at most 51-bit random integers. Intermediate results of the algorithms can get one bit longer than the parameters, and MATLAB stores integers in double precision floating point registers, having 52-bit precision.
Modular inverse based on the built in extended GCD algorithm:
function C = modinv(a,m)
%modinv(a,m) -> inverse of a mod m, based on Euclidean extended GCD alg
[G,C,D] = gcd(a,m);
if G~= 1, C = 0; return, end
if C < 0, C = C + m; end
Binary Right-shift algorithm, with signed arithmetic:
function R = RSSinv(a,m)
%RSSinv(a,m) -> inverse of a mod m, Right-Shift binary alg
% m = odd, Signed arithmetic
if ~mod(m,2) || a==0, R = 0; return, end
U = m; V = a;
R = 0; S = 1;
while V > 0
if ~mod(U,2), U = U/2; % U is even
if ~mod(R,2), R = R/2;
else R = (R+m)/2; end % R = full length
elseif ~mod(V,2), V = V/2; % V is even
if ~mod(S,2), S = S/2;
else S = (S+m)/2; end % S = full length
else % U, V odd
if U > V, U = U - V; R = R - S;
else V = V - U; S = S - R; end
end
end
if U > 1, R = 0; return,end
while R < 0, R = R + m; end % loops at most twice (-2m < R < 3m)
while R > m, R = R - m; end % R can grow longer than m
Binary Right-shift algorithm, with unsigned arithmetic:
function R = RSUinv(a,m)
%RSUinv(a,m) -> inverse of a mod m, Right-Shift binary alg, using Unsigned
% m = odd
if ~mod(m,2), R = 0; return, end
U = m; V = a;
R = 0; S = 1;
while V > 0
if ~mod(U,2) % U is even
U = U/2;
if ~mod(R,2), R = R/2;
else R = (R+m)/2; end % R = full length
elseif ~mod(V,2) % V is even
V = V/2;
if ~mod(S,2), S = S/2;
else S = (S+m)/2; end % S = full length
else % U, V odd
if U > V, U = U - V;
if R < S, R = R - S + m; else R = R - S; end
else, V = V - U;
if S < R, S = S - R + m; else S = S - R; end
end
end
end
if U > 1, R = 0; return, end
Shifting Euclidean algorithm with unsigned arithmetic:
function S = SEUinv(a,m)
%SEUinv(a,m) -> inverse of a mod m, Shifting Euclidean alg, using Unsigned
if a < m
U = m; V = a;
R = 0; S = 1;
else
V = m; U = a;
S = 0; R = 1;
end
while V > 1
f = bits(U) - bits(V); % U >= V, f >= 0
if U < bitshift(V,f), f = f - 1; end
U = U - bitshift(V,f); % always U >= 0
t = S;
for i = 1:f
t = t+t; % #adds <= bits(U)+bits(V)
if t > m, t = t - m; end
end
R = R - t; % check R < t beforhand
if R < 0, R = R + m; end % one of R,S gets large soon
if U < V
t = U; U = V; V = t; % swap(U,V)
t = R; R = S; S = t; % swap(R,S)
end
end
if V == 0, S = 0; end
It uses the simple function bits(n), telling the number of bits n is represented with:
function b = bits(n)
[x,b] = log2(n);
Shifting Euclidean algorithm with signed arithmetic:
function S = SESinv(a,m)
%SESinv(a,m) -> inverse of a mod m, Shifting Euclidean alg, using Signed
if a < m, U = m; V = a; R = 0; S = 1;
else V = m; U = a; S = 0; R = 1; end
while abs(V) > 1
f = bits(U) - bits(V);
if sign(U) == sign(V)
U = U - shift(V,f);
R = R - shift(S,f);
else
U = U + shift(V,f);
R = R + shift(S,f);
end
if abs(U) < abs(V)
t = U; U = V; V = t; % swap(U,V)
t = R; R = S; S = t; % swap(R,S)
end
end
if V== 0, S = 0;
elseif V < 0, S = -S; end
if S > m, S = S - m;
elseif S < 0, S = S + m; end
It uses also the simple function shift(n,i), shifting the signed integer n by i bits to the left:
function s = shift(n,i)
s = sign(n)*bitshift(abs(n),i);
Double Plus-minus Right-shift algorithm (RS2+−) using signed arithmetic:
function R = R2Sinv(a,m)
%R2Sinv(a,m) -> inverse of a mod m, Right-shift 2-fold +/- binary alg
% m = odd, using Signed arithmetic
if ~mod(m,2) || a==0, R = 0; return, end
U = m; V = a;
R = 0; S = 1;
Q = mod(m,4); m2 = 2*m;
while ~mod(V,2) % for each trailing 0 bit
V = V/2;
if ~mod(S,2), S = S/2;
elseif S > m, S = (S-m)/2;
else S = (S+m)/2; end
end
while 1 % U,V = odd
if U > V
if mod(U,4) == mod(V,4)
U = U - V; R = R - S;
else
U = U + V; R = R + S;
end
U = U/4; T = mod(R,4);
if T == 0, R = R/4;
elseif T == 2, R = (R+m2)/4;
elseif T == Q, R = (R-m)/4;
else R = (R+m)/4; end
while ~mod(U,2),U = U/2;
if ~mod(R,2),R = R/2;
elseif R > m,R = (R-m)/2;
else R = (R+m)/2; end
end % |R| < m
else
if mod(U,4) == mod(V,4)
V = V - U; S = S - R;
else
V = V + U; S = S + R;
end
if V == 0, break; end
V = V/4; T = mod(S,4);
if T == 0, S = S/4;
elseif T == 2, S = (S+m2)/4;
elseif T == Q, S = (S-m)/4;
else S = (S+m)/4; end
while ~mod(V,2),V = V/2;
if ~mod(S,2),S = S/2;
elseif S > m,S = (S-m)/2;
else S = (S+m)/2; end
end
end
end
if U > 1, R = 0; return, end
if R < 0, R = R + m; end
Algorithm
Steps/bit
|
Right-Shift
|
Left-Shift
|
Shift-Euclidean
|
RS1
|
RS+−
|
RS2+−
|
RSDH
|
RSDH+−
|
LS1
|
LS3
|
SE
|
SE3
|
Iterations
|
0.7045 n
|
0.6115 n
|
0.6115 n
|
0.7045 n
|
0.6115 n
|
0.7650 n
|
0.6646 n
|
0.7684 n
|
0.6744 n
|
UV shift cost
|
0.3531 n2
−1.2200 n
|
0.3065 n2
−1.1891 n
|
0.3065 n2
−1.1891 n
|
0.3531 n2
−1.2200 n
|
0.3065 n2
−1.1891 n
|
0.3834 n2
−0.8836 n
|
0.3967 n2
−0.8435 n
|
0.3101 n2
−1.0646 n
|
0.2708 n2
−0.8742 n
|
RS shift cost
|
1.0592 n2
−4.9984 n
|
1.2259 n2
−5.2592 n
|
0.9808 n2
−5.1720 n
|
0.9241 n2
−3.3945 n
|
0.8021 n2
−3.3794 n
|
0.5300 n2
−4.9665 n
|
0.5558 n2
−5.1855 n
|
0.3101 n2
−2.9784 n
|
0.2708 n2
−2.5787 n
|
Total shift cost
|
1.4123 n2
−6.2184 n
|
1.5324 n2
−6.4483 n
|
1.2873 n2
−6.3611 n
|
1.2772 n2
−4.6145 n
|
1.1086 n2
−4.5685 n
|
0.9134 n2
−5.8501 n
|
0.9525 n2
−6.0290 n
|
0.6202 n2
−4.0430 n
|
0.5416 n2
−3.4529 n
|
UV subtract cost
|
0.3531 n2
+0.2658 n
|
0.3065 n2
+0.2967 n
|
0.3065 n2
+0.2967 n
|
0.3531 n2
+0.2658 n
|
0.3065 n2
+0.2967 n
|
0.3835 n2
+0.4377 n
|
0.3331 n2
+0.5942 n
|
0.3851 n2
+0.4276 n
|
0.3380 n2
+0.4958 n
|
RS subtract cost
|
1.4123 n2
−4.8065 n
|
1.5325 n2
−4.8844 n
|
1.2873 n2
−4.5004 n
|
0.9241 n2
−1.4559 n
|
0.8021 n2
−0.7786 n
|
0.3834 n2
−1.0101 n
|
0.3331 n2
−0.9160 n
|
0.3851 n2
−1.0331 n
|
0.3380 n2
−0.7125 n
|
Total subtract cost
|
1.7654 n2
−4.5407 n
|
1.8390 n2
−4.5877 n
|
1.5938 n2
−4.2037 n
|
1.2772 n2
−1.1901 n
|
1.1086 n2
−0.4819 n
|
0.7669 n2
−0.5724 n
|
0.6662 n2
−0.3218 n
|
0.7702 n2
−0.6055 n
|
0.6760 n2
−0.2167 n
|
Complexity @
0 cost shift
|
1.7654 n2
|
1.8390 n2
|
1.5938 n2
|
1.2772 n2
|
1.1086 n2
|
0.7669 n2
|
0.6662 n2
|
0.7702 n2
|
0.6750 n2
|
Complexity @
¼ add cost shift
|
2.1185 n2
|
2.2221 n2
|
1.9156 n2
|
1.5965 n2
|
1.3858 n2
|
0.9953 n2
|
0.9043 n2
|
0.9253 n2
|
0.8114 n2
|
Complexity @
1 add cost shift
|
3.1777 n2
|
3.3714 n2
|
2.8811 n2
|
2.5544 n2
|
2.2172 n2
|
1.6803 n2
|
1.6187 n2
|
1.3904 n2
|
1.2176 n2
|
UV shifts by 1
|
0.3522 n
|
−
|
−
|
0.3522 n
|
−
|
0.1983 n
|
0.1977 n
|
0.2576 n
|
0.2143 n
|
UV shifts by 2
|
0.1761 n
|
0.3058 n
|
0.3058 n
|
0.1761 n
|
0.3058 n
|
0.2463 n
|
0.2388 n
|
0.1705 n
|
0.1573 n
|
UV shifts by 3
|
0.0881 n
|
0.1529 n
|
0.1529 n
|
0.0881 n
|
0.1529 n
|
0.1516 n
|
0.1778 n
|
0.0927 n
|
0.0831 n
|
Longer UV shifts
|
0.0881 n
|
0.1529 n
|
0.1529 n
|
0.0881 n
|
0.1529 n
|
0.1689 n
|
0.1772 n
|
0.0980 n
|
0.0857 n
|
RS shifts by 1
|
0.7925 n
|
0.7644 n
|
0.3364 n
|
0.6375 n
|
−
|
0.5202 n
|
0.5395 n
|
0.2576 n
|
0.2143 n
|
RS shifts by 2
|
0.1982 n
|
0.3440 n
|
0.4816 n
|
0.3188 n
|
0.5534 n
|
0.3142 n
|
0.3313 n
|
0.1705 n
|
0.1573 n
|
RS shifts by 3
|
0.0495 n
|
0.0860 n
|
0.1204 n
|
0.1594 n
|
0.2767 n
|
0.1280 n
|
0.1413 n
|
0.0927 n
|
0.0831 n
|
Longer RS shifts
|
0.0165 n
|
0.0287 n
|
0.0401 n
|
0.1594 n
|
0.2767 n
|
0.0952 n
|
0.0968 n
|
0.0980 n
|
0.0857 n
|
- /- 1/28/2017
Share with your friends: |