Cipolla's algorithm: Difference between revisions

m
→‎GMP version: a few small tweaks
m (→‎{{header|Sage}}: Changed - to +, also a few efficiency updates for dealing with large numbers)
m (→‎GMP version: a few small tweaks)
Line 364:
===GMP version===
{{libheader|GMP}}
<lang freebasic>' version 1012-04-2017
' compile with: fbc -s console
 
Line 370:
 
Type fp2
x As mpz_ptrMpz_ptr
y As mpz_ptrMpz_ptr
End Type
 
Line 383:
Data "34035243914635549601583369544560650254325084643201" ', 10^50 + 151
 
Function fp2mul(a As fp2, b As fp2, p As mpz_ptrMpz_ptr, w2 As mpz_ptrMpz_ptr) As fp2
 
Dim As fp2 r
r.x = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(r.x)
r.y = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(r.y)
 
mpz_mulMpz_mul (r.x, a.y, b.y)
mpz_mulMpz_mul (r.x, r.x, w2)
mpz_addmulMpz_addmul(r.x, a.x, b.x)
mpz_modMpz_mod (r.x, r.x, p)
mpz_mulMpz_mul (r.y, a.x, b.y)
mpz_addmulMpz_addmul(r.y, a.y, b.x)
mpz_modMpz_mod (r.y, r.y, p)
 
Return r
Line 401:
End Function
 
Function fp2square(a As fp2, p As mpz_ptrMpz_ptr, w2 As mpz_ptrMpz_ptr) As fp2
 
Return fp2mul(a, a, p, w2)
Line 407:
End Function
 
Function fp2pow(a As fp2, n As mpz_ptrMpz_ptr, p As mpz_ptrMpz_ptr, w2 As mpz_ptrMpz_ptr) As fp2
 
If mpz_cmp_uiMpz_cmp_ui(n, 0) = 0 Then
mpz_set_uiMpz_set_ui(a.x, 1)
mpz_set_uiMpz_set_ui(a.y, 0)
Return a
End If
If mpz_cmp_uiMpz_cmp_ui(n, 1) = 0 Then Return a
If mpz_cmp_uiMpz_cmp_ui(n, 2) = 0 Then Return fp2square(a, p, w2)
If mpz_tstbitMpz_tstbit(n, 0) = 0 Then
mpz_fdiv_q_2expMpz_fdiv_q_2exp(n, n, 1) ' even
Return fp2square(fp2pow(a, n, p, w2), p, w2)
Else
mpz_sub_uiMpz_sub_ui(n, n, 1) ' odd
Return fp2mul(a, fp2pow(a, n, p, w2), p, w2)
End If
Line 428:
' ------=< MAIN >=------
 
Dim As Long i, result
Dim As ZString Ptr zstr
Dim As String n_str, p_str
 
Dim As mpz_ptrMpz_ptr a, n, p, p2, w2, x1, x2
a = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(a)
n = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(n)
p = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(p)
p2 = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(p2)
w2 = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(w2)
x1 = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(x1)
x2 = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(x2)
 
Dim As fp2 answer
answer.x = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(answer.x)
answer.y = Allocate(Len(__mpz_struct__Mpz_struct)) : Mpz_init(answer.y)
 
For i = 1 To 8
Read n_str
mpz_set_strMpz_set_str(n, n_str, 10)
If i < 8 Then
Read p_str
mpz_set_strMpz_set_str(p, p_str, 10)
Else
p_str = "10^50 + 151" ' set up last n
mpz_set_strMpz_set_str(p, "1" + String(50, "0"), 10)
mpz_add_uiMpz_add_ui(p, p, 151)
End If
 
Print "Find solution for n = "; n_str; " and p = "; p_str
 
resultIf Mpz_tstbit(p, 0) = mpz_probab_prime_p0 OrElse Mpz_probab_prime_p(p, 20) = 0 Then
If mpz_tstbit(p, 0) = 0 OrElse result = 0 Then
Print p_str; "is not a odd prime"
Print
Line 467 ⟶ 466:
 
' p is checked and is a odd prime
result' =legendre mpz_legendre(n,symbol p)needs 'to need tobe 1
If resultMpz_legendre(n, p) <> 1 Then
' p is checked and is a odd prime
If result <> 1 Then
Print n_str; " is not a square in F"; p_str
Print
Line 475 ⟶ 473:
End If
 
mpz_set_uiMpz_set_ui(a, 1)
Do
Do
Do
mpz_add_uiMpz_add_ui(a, a, 1)
mpz_mulMpz_mul(w2, a, a)
mpz_subMpz_sub(w2, w2, n)
Loop Until mpz_legendreMpz_legendre(w2, p) = -1
 
mpz_setMpz_set(answer.x, a)
mpz_set_uiMpz_set_ui(answer.y, 1)
mpz_add_uiMpz_add_ui(p2, p, 1) ' p2 = p + 1
mpz_fdiv_q_2expMpz_fdiv_q_2exp(p2, p2, 1) ' p2 = p2 \ 2 (p2 shr 1)
 
answer = fp2pow(answer, p2, p, w2)
 
Loop Until mpz_cmp_uiMpz_cmp_ui(answer.y, 0) = 0
mpz_setMpz_set(x1, answer.x)
mpz_subMpz_sub(x2, p, x1)
mpz_powm_uiMpz_powm_ui(a, x1, 2, p)
mpz_powm_uiMpz_powm_ui(p2, x2, 2, p)
If mpz_cmpMpz_cmp(a, n) = 0 AndAlso mpz_cmpMpz_cmp(p2, n) = 0 Then Exit Do
Loop
 
Line 506 ⟶ 504:
Next
 
mpz_clearMpz_clear(x1) : mpz_clearMpz_clear(p2) : mpz_clearMpz_clear(p) : mpz_clearMpz_clear(a) : mpz_clearMpz_clear(n)
mpz_clearMpz_clear(x2) : mpz_clearMpz_clear(w2) : mpz_clearMpz_clear(answer.x) : mpz_clearMpz_clear(answer.y)
 
' empty keyboard buffer
457

edits