PureBasic Forums - English

Code: Select all

Procedure.l gcd(a,b)
 While b <> 0
  t = b
  b = a % b
  a = t
 Wend
 ProcedureReturn a
EndProcedure

Procedure.l Shanks(a.l)
 Protected qm1.l
 Protected q.l
 Protected qp1.l
 Protected pm1.l
 Protected p.l
 Protected b.l
 Protected test.f
 
 pm1 = Int(Sqr(a))
 qm1 = 1
 q = a-pm1*pm1
 
 Repeat
  b = Int((Int(Sqr(a))+pm1)/q)
  p = b*q-pm1
  qp1 = qm1+b*(pm1-p)
  pm1 = p
  qm1 = q
  q = qp1
 Until Int(Sqr(q)) = Sqr(q)
 
 b = Int((Int(Sqr(a))-pm1)/Sqr(q))
 pm1 = b*Sqr(q)+pm1
 qm1 = Sqr(q)
 q = (a-pm1*pm1)/qm1
 
 Repeat
  b = Int((Int(Sqr(a))+pm1)/q)
  p = b*q-pm1
  qp1 = qm1+b*(pm1-p)
  pm1 = p
  qm1 = q
  q = qp1
 Until pm1 = p
 
 ProcedureReturn gcd(a,p)
EndProcedure

Code: Select all

Procedure.i MinPrimeFactor (n)
   ; in : integer >= 2
   ; out: smallest prime factor of n, or 1 if n is a prime number
   Protected k

   If n = 2
      ProcedureReturn 1
   EndIf

   If n % 2 = 0
      ProcedureReturn 2
   EndIf

   For k = 3 To Int(Sqr(n)) Step 2
      If n % k = 0
         ProcedureReturn k
      EndIf
   Next

   ProcedureReturn 1
EndProcedure

Code: Select all

; FactShor  AKJ  31-Mar-09
; Factorise using Shor's Algorithm
; http://en.wikipedia.org/wiki/Shor%27s_algorithm
; Adapted from my Amstrad CPC Locomotive Basic program of 12-Mar-05

; Peter Shor's factoring algorithm:
; 1. Pick a random number a < N
; 2. Compute gcd(a, N). This may be done using the Euclidean algorithm.
; 3. If gcd(a, N)<>1, then there is a non-trivial factor of N, so we are done.
; 4. Otherwise, find r, the period of the function:  f(x) = a^x mod N
;    i.e. the order r of a
;    or the smallest positive integer r for which f(x + r) = f(x) .
; 5. If r is odd, go back to step 1.
; 6. Let h = r/2 = half-period.
; 7. If a^h = -1 (mod N), go back to step 1.
; 8. gcd(a^h ±1, N) is a nontrivial factor of N. We are done.

; AKJ modification to step 4 [and step 1], above:
;    If a is prime, the period r can be found by setting x=0 .
;    That is, by finding r>0 such that f(r) = f(0) = 1 .
;    Thus r is the first non-trivial solution to a^r mod N = 1 .
;    This requires that step 1 must pick only prime numbers, but
;    it is simpler to merely loop through the first few prime numbers.
; Step 4 is best done on a quantum computer.


EnableExplicit

Declare   Factor(N.q)
Declare   Euclid(i.q, j.q)
Declare.i Period(a, N.q)
Declare.q Power(a.q, x, N.q)
Declare   Result(f.q)

;-Main program
Define N$, N.q
OpenConsole()
Repeat
  Print("Number to factor (or Enter to end) : ")
  N$ = Input()
  N=Val(N$)
  If N<1: Break: EndIf
  Factor(N)
  PrintN("")
ForEver
CloseConsole()
End


Procedure Factor(N.q)
Protected found.b, e.q, f.q, i, p, a, r, h, time
If N=1: ProcedureReturn Result(1): EndIf
PrintN("Factoring "+Str(N))
; First, try small prime factors
found = #False
; 614889782588491410 = product of all prime numbers 1..49
e = Euclid(614889782588491410, N)
If e>1: found = #True
  Restore Primes1
  For i = 1 To 15
    Read.q p ; Prime number
    If e%p=0: Result(p): EndIf
  Next i
EndIf
; 3749562977351496827 = product of all prime numbers 50..99
e = Euclid(3749562977351496827, N)
If e>1: found = #True
  Restore Primes2
  For i = 1 To 15
    Read.q p ; Prime number
    If e%p=0: Result(p): EndIf
  Next i
EndIf
If found: ProcedureReturn: EndIf
; No small factors exist, so use Shor's algorithm with AKJ modifications
Restore Primes1
For i = 1 To 25 ; Make many factorisation attempts
  Read.q a
  time = ElapsedMilliseconds()
  r = Period(a, N) ; Period for the trial base a
  time - ElapsedMilliseconds()
  If time: PrintN("Time = "+Str(-time)+" milliseconds"): EndIf
  h = r>>1 ; Half period r/2
  PrintN("Half period = "+Str(h))
  f = 0
  If r&1=0 ; If r is even
    f = Power(a, h, N)
    PrintN(Str(a)+"^(half period) = "+Str(f))
    If f<>N-1 ; I.e. If a^h<>-1 mod N
      Result(Euclid(N, f-1))
      Result(Euclid(N, f+1))
      ProcedureReturn
    EndIf ; Factor found
  EndIf
Next i
Result(0) ; Failed
EndProcedure


Procedure.q Euclid(i.q, j.q) ; Best if i>j
Repeat
  i % j: Swap i, j
Until j=0
ProcedureReturn i
EndProcedure


Procedure.i Period(a, N.q)
; Return the period r of a^x mod N
Protected r, f.q
PrintN("")
PrintN("Trying base "+Str(a))
r = 0 ; Period
f = 1 ; a^r mod N
Select a
Case 2 ; Faster than the default case
  Repeat
    Repeat
      r + 1: f + f
    Until f>=N
    f - N ; f%N
  Until f = 1
; Case 3 ; Slightly faster than the default case
;   Repeat
;     r + 1: f + f+f
;     While f>=N: f - N: Wend ; f%N
;   Until f=1
Default
  Repeat
    r + 1: f = (f*a)%N
  Until f=1
EndSelect
ProcedureReturn r
EndProcedure


Procedure.q Power(a.q, x, N.q)
; Returns a^x mod N
; Uses the binary representation of x for speed
Protected p.q=1
a%N
Repeat
  If x&1: p = (p*a)%N: EndIf ; If the current bit of x is a 1
  x>>1 ; x/2
  If x: a = (a*a)%N: EndIf
Until x=0
ProcedureReturn p
EndProcedure


Procedure Result(f.q)
If f<=0
  PrintN("Factor not found")
ElseIf f=1
  ; Do nothing
ElseIf f<100
  PrintN("Prime factor is "+Str(f))
Else
  PrintN("Factor (possibly composite) is "+Str(f))
EndIf
EndProcedure


DataSection
Primes1: ; 15
Data.q 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47
Primes2: ; 10
Data.q 53, 59, 61, 67, 71, 73, 79, 83, 89, 97
EndDataSection