58 lines
1.1 KiB
Rexx
58 lines
1.1 KiB
Rexx
/* REXX (required by some interpreters) */
|
|
Numeric Digits 1000000
|
|
ttest ='[(10, 13), (56, 101), (1030, 10009), (44402, 100049)]'
|
|
Do While pos('(',ttest)>0
|
|
Parse Var ttest '(' n ',' p ')' ttest
|
|
r = tonelli(n, p)
|
|
Say "n =" n "p =" p
|
|
Say " roots :" r (p - r)
|
|
End
|
|
Exit
|
|
|
|
legendre: Procedure
|
|
Parse Arg a, p
|
|
return pow(a, (p - 1) % 2, p)
|
|
|
|
tonelli: Procedure
|
|
Parse Arg n, p
|
|
q = p - 1
|
|
s = 0
|
|
Do while q // 2 == 0
|
|
q = q % 2
|
|
s = s+1
|
|
End
|
|
if s == 1 Then
|
|
return pow(n, (p + 1) % 4, p)
|
|
Do z=2 To p
|
|
if p - 1 == legendre(z, p) Then
|
|
Leave
|
|
End
|
|
c = pow(z, q, p)
|
|
r = pow(n, (q + 1) / 2, p)
|
|
t = pow(n, q, p)
|
|
m = s
|
|
t2 = 0
|
|
Do while (t - 1) // p <> 0
|
|
t2 = (t * t) // p
|
|
Do i=1 To m
|
|
if (t2 - 1) // p == 0 Then
|
|
Leave
|
|
t2 = (t2 * t2) // p
|
|
End
|
|
y=2**(m - i - 1)
|
|
b = pow(c, y, p)
|
|
If b=10008 Then Trace ?R
|
|
r = (r * b) // p
|
|
c = (b * b) // p
|
|
t = (t * c) // p
|
|
m = i
|
|
End
|
|
return r
|
|
pow: Procedure
|
|
Parse Arg x,y,z
|
|
If y>0 Then
|
|
p=x**y
|
|
Else p=x
|
|
If z>'' Then
|
|
p=p//z
|
|
Return p
|