Cipolla's algorithm: Difference between revisions
Content added Content deleted
m (→{{header|Sage}}: Replaced broken Sage code in previous edit, now updating for efficiency and errors) |
(→{{header|FreeBASIC}}: added GMP version) |
||
Line 151: | Line 151: | ||
</pre> |
</pre> |
||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |
||
===LongInt version=== |
|||
Had a close look at the EchoLisp code for step 2. |
Had a close look at the EchoLisp code for step 2. |
||
Used the FreeBASIC code from the Miller-Rabin task for prime testing. |
Used the FreeBASIC code from the Miller-Rabin task for prime testing. |
||
Line 171: | Line 172: | ||
Function mul_mod(a As ULongInt, b As ULongInt, modulus As ULongInt) As ULongInt |
Function mul_mod(a As ULongInt, b As ULongInt, modulus As ULongInt) As ULongInt |
||
' returns a * b mod modulus |
' returns a * b mod modulus |
||
Dim As ULongInt x, y = |
Dim As ULongInt x, y = a mod modulus |
||
While b > 0 |
While b > 0 |
||
Line 281: | Line 282: | ||
If n = 2 Then Return fp2square(a, p, w2) |
If n = 2 Then Return fp2square(a, p, w2) |
||
If (n And 1) = 0 Then |
If (n And 1) = 0 Then |
||
Return fp2square(fp2pow(a, n\2, p, w2), p , w2) |
Return fp2square(fp2pow(a, n \ 2, p, w2), p , w2) |
||
Else |
Else |
||
Return fp2mul(a, fp2pow(a, n -1, p, w2), p, w2) |
Return fp2mul(a, fp2pow(a, n -1, p, w2), p, w2) |
||
Line 361: | Line 362: | ||
Find solution for n = 881398088036 and p = 1000000000039 |
Find solution for n = 881398088036 and p = 1000000000039 |
||
Solution found: x1 = 791399408049, x2 = 208600591990</pre> |
Solution found: x1 = 791399408049, x2 = 208600591990</pre> |
||
===GMP version=== |
|||
{{libheader|GMP}} |
|||
<lang freebasic>' version 10-04-2017 |
|||
' compile with: fbc -s console |
|||
#Include Once "gmp.bi" |
|||
Type fp2 |
|||
x As mpz_ptr |
|||
y As mpz_ptr |
|||
End Type |
|||
Data "10", "13" |
|||
Data "56", "101" |
|||
Data "8218", "10007" |
|||
Data "8219", "10007" |
|||
Data "331575", "1000003" |
|||
Data "665165880", "1000000007" |
|||
Data "881398088036", "1000000000039" |
|||
Data "34035243914635549601583369544560650254325084643201" ', 10^50 + 151 |
|||
Function fp2mul(a As fp2, b As fp2, p As mpz_ptr, w2 As mpz_ptr) As fp2 |
|||
Dim As fp2 r |
|||
r.x = Allocate(Len(__mpz_struct)) : Mpz_init(r.x) |
|||
r.y = Allocate(Len(__mpz_struct)) : Mpz_init(r.y) |
|||
mpz_mul (r.x, a.y, b.y) |
|||
mpz_mul (r.x, r.x, w2) |
|||
mpz_addmul(r.x, a.x, b.x) |
|||
mpz_mod (r.x, r.x, p) |
|||
mpz_mul (r.y, a.x, b.y) |
|||
mpz_addmul(r.y, a.y, b.x) |
|||
mpz_mod (r.y, r.y, p) |
|||
Return r |
|||
End Function |
|||
Function fp2square(a As fp2, p As mpz_ptr, w2 As mpz_ptr) As fp2 |
|||
Return fp2mul(a, a, p, w2) |
|||
End Function |
|||
Function fp2pow(a As fp2, n As mpz_ptr, p As mpz_ptr, w2 As mpz_ptr) As fp2 |
|||
If mpz_cmp_ui(n, 0) = 0 Then |
|||
mpz_set_ui(a.x, 1) |
|||
mpz_set_ui(a.y, 0) |
|||
Return a |
|||
End If |
|||
If mpz_cmp_ui(n, 1) = 0 Then Return a |
|||
If mpz_cmp_ui(n, 2) = 0 Then Return fp2square(a, p, w2) |
|||
If mpz_tstbit(n, 0) = 0 Then |
|||
mpz_fdiv_q_2exp(n, n, 1) ' even |
|||
Return fp2square(fp2pow(a, n, p, w2), p, w2) |
|||
Else |
|||
mpz_sub_ui(n, n, 1) ' odd |
|||
Return fp2mul(a, fp2pow(a, n, p, w2), p, w2) |
|||
End If |
|||
End Function |
|||
' ------=< MAIN >=------ |
|||
Dim As Long i, result |
|||
Dim As ZString Ptr zstr |
|||
Dim As String n_str, p_str |
|||
Dim As mpz_ptr a, n, p, p2, w2, x1, x2 |
|||
a = Allocate(Len(__mpz_struct)) : Mpz_init(a) |
|||
n = Allocate(Len(__mpz_struct)) : Mpz_init(n) |
|||
p = Allocate(Len(__mpz_struct)) : Mpz_init(p) |
|||
p2 = Allocate(Len(__mpz_struct)) : Mpz_init(p2) |
|||
w2 = Allocate(Len(__mpz_struct)) : Mpz_init(w2) |
|||
x1 = Allocate(Len(__mpz_struct)) : Mpz_init(x1) |
|||
x2 = Allocate(Len(__mpz_struct)) : Mpz_init(x2) |
|||
Dim As fp2 answer |
|||
answer.x = Allocate(Len(__mpz_struct)) : Mpz_init(answer.x) |
|||
answer.y = Allocate(Len(__mpz_struct)) : Mpz_init(answer.y) |
|||
For i = 1 To 8 |
|||
Read n_str |
|||
mpz_set_str(n, n_str, 10) |
|||
If i < 8 Then |
|||
Read p_str |
|||
mpz_set_str(p, p_str, 10) |
|||
Else |
|||
p_str = "10^50 + 151" ' set up last n |
|||
mpz_set_str(p, "1" + String(50, "0"), 10) |
|||
mpz_add_ui(p, p, 151) |
|||
End If |
|||
Print "Find solution for n = ";n_str; " and p = ";p_str |
|||
result = mpz_probab_prime_p(p, 20) |
|||
If mpz_tstbit(p, 0) = 0 OrElse result = 0 Then |
|||
Print p_str; "is not a odd prime" |
|||
Print |
|||
Continue For |
|||
End If |
|||
' p is checked and is a odd prime |
|||
result = mpz_legendre(n, p) ' need to 1 |
|||
' p is checked and is a odd prime |
|||
If result <> 1 Then |
|||
Print n_str; " is not a square in F"; p_str |
|||
Print |
|||
Continue For |
|||
End If |
|||
mpz_set_ui(a, 1) |
|||
Do |
|||
Do |
|||
Do |
|||
mpz_add_ui(a, a, 1) |
|||
mpz_mul(w2, a, a) |
|||
mpz_sub(w2, w2, n) |
|||
Loop Until mpz_legendre(w2, p) = -1 |
|||
mpz_set(answer.x, a) |
|||
mpz_set_ui(answer.y, 1) |
|||
mpz_add_ui(p2, p, 1) ' p2 = p + 1 |
|||
mpz_fdiv_q_2exp(p2, p2, 1) ' p2 = p2 \ 2 (p2 shr 1) |
|||
answer = fp2pow(answer, p2, p, w2) |
|||
Loop Until mpz_cmp_ui(answer.y, 0) = 0 |
|||
mpz_set(x1, answer.x) |
|||
mpz_sub(x2, p, x1) |
|||
mpz_powm_ui(a, x1, 2, p) |
|||
mpz_powm_ui(p2, x2, 2, p) |
|||
If mpz_cmp(a, n) = 0 AndAlso mpz_cmp(p2, n) = 0 Then Exit Do |
|||
Loop |
|||
zstr = Mpz_get_str(0, 10, x1) |
|||
Print "Solution found: x1 = "; *zstr; |
|||
zstr = Mpz_get_str(0, 10, x2) |
|||
Print ", x2 = "; *zstr |
|||
Print |
|||
Next |
|||
mpz_clear(x1) : mpz_clear(p2) : mpz_clear(p) : mpz_clear(a) : mpz_clear(n) |
|||
mpz_clear(x2) : mpz_clear(w2) : mpz_clear(answer.x) : mpz_clear(answer.y) |
|||
' empty keyboard buffer |
|||
While Inkey <> "" : Wend |
|||
Print : Print "hit any key to end program" |
|||
Sleep |
|||
End</lang> |
|||
{{out}} |
|||
<pre>Find solution for n = 10 and p = 13 |
|||
Solution found: x1 = 6, x2 = 7 |
|||
Find solution for n = 56 and p = 101 |
|||
Solution found: x1 = 37, x2 = 64 |
|||
Find solution for n = 8218 and p = 10007 |
|||
Solution found: x1 = 9872, x2 = 135 |
|||
Find solution for n = 8219 and p = 10007 |
|||
8219 is not a square in F10007 |
|||
Find solution for n = 331575 and p = 1000003 |
|||
Solution found: x1 = 855842, x2 = 144161 |
|||
Find solution for n = 665165880 and p = 1000000007 |
|||
Solution found: x1 = 524868305, x2 = 475131702 |
|||
Find solution for n = 881398088036 and p = 1000000000039 |
|||
Solution found: x1 = 208600591990, x2 = 791399408049 |
|||
Find solution for n = 34035243914635549601583369544560650254325084643201 and p = 10^50 + 151 |
|||
Solution found: x1 = 17436881171909637738621006042549786426312886309400, x2 = 82563118828090362261378993957450213573687113690751</pre> |
|||
=={{header|Go}}== |
=={{header|Go}}== |