Cipolla's algorithm: Difference between revisions
→{{header|FreeBASIC}}: added GMP version
m (→{{header|Sage}}: Replaced broken Sage code in previous edit, now updating for efficiency and errors) |
(→{{header|FreeBASIC}}: added GMP version) |
||
Line 151:
</pre>
=={{header|FreeBASIC}}==
===LongInt version===
Had a close look at the EchoLisp code for step 2.
Used the FreeBASIC code from the Miller-Rabin task for prime testing.
Line 171 ⟶ 172:
Function mul_mod(a As ULongInt, b As ULongInt, modulus As ULongInt) As ULongInt
' returns a * b mod modulus
Dim As ULongInt x, y =
While b > 0
Line 281 ⟶ 282:
If n = 2 Then Return fp2square(a, p, w2)
If (n And 1) = 0 Then
Return fp2square(fp2pow(a, n \ 2, p, w2), p , w2)
Else
Return fp2mul(a, fp2pow(a, n -1, p, w2), p, w2)
Line 361 ⟶ 362:
Find solution for n = 881398088036 and p = 1000000000039
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}}==
| |||