# Calculate the least non-negative remainder when an integer a # is divided by a positive integer b. sub mod { my (\$a, \$b) = @_; my \$c = \$a % \$b; return \$c if (\$a >= 0); return 0 if (\$c == 0); return (\$c + \$b); } # Calculate a^b (mod c), where a,b and c are integers and a,b>=0, c>=1 sub mpow { my (\$a, \$b, \$c) = @_; my (\$x, \$y, \$z) = (\$a, \$b, 1); while (\$y > 0) { while (\$y % 2 == 0) { \$y = \$y / 2; \$x = \$x**2 % \$c; } \$y--; \$z = mod(\$z * \$x, \$c); } return \$z; } # Shanks-Tonelli algorithm to calculate y^2 = a (mod p) for p an odd prime sub tonelli { my (\$a, \$p) = @_; my (\$b, \$e, \$g, \$h, \$i, \$m, \$n, \$q, \$r, \$s, \$t, \$x, \$y, \$z); \$q = \$p - 1; \$t = mpow(\$a, \$q/2, \$p); return 0 if (\$t = \$q); # a is a quadratic non-residue mod p \$s = \$q; \$e = 0; while (\$s % 2 == 0) { \$s = \$s/2; \$e++; } # p-1 = s * 2^e; \$x = mpow(\$a, ++\$s/2, \$p); \$b = mpow(\$a, \$s, \$p); return \$x if (\$b == 1); \$n = 2; RES: while (\$n >= 2) { \$t = mpow(\$n, \$q/2, \$p); last RES if (\$t == \$q); } continue { \$n++; } # n is the least quadratic non-residue mod p \$g = mpow(\$n, \$s, \$p); \$r = \$e; OUT: { do { \$y = \$b; \$m = 0; IN: while (\$m <= \$r-1) { last IN if (\$y == 1); \$y = \$y**2 % \$p; } continue { \$m++; } return \$x if (\$m == 0); \$z = \$r - \$m - 1; \$h = \$g; for (\$i = 1; \$i < \$z; \$i++) {\$h = \$h**2 % \$p;} \$x = (\$x * \$h) % \$p; \$b = (\$b * \$h**2) % \$p; \$g = \$h**2 % \$p; \$r = \$m; } while 1; } }