Here's another factoring method, but it will not be particularly quick on your slow computer as it is intended to be used only on highly-parallel quantum computers:
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
