Program: Findprim.bas
Implementation Select: Implementation
Operation Select: Program: FINDPRIM.BAS
DEFDBL A-Z
DECLARE SUB XAMODN (a#, x#, n#, asave#, lc#)
DECLARE SUB FACTOR1 (level#, lc#, n#, n1#)
DECLARE FUNCTION GCD# (a1#, b1#, lc#)
DECLARE SUB CALCXN (ipos#, xn#)
DECLARE SUB HTML (id#)
DECLARE SUB FACTOR2 (n#, lc#, n1#, n2#)
DECLARE SUB FACTOR (level#, lc#, n1#, n#)
' PROGRAM FINDPRIM
' Factorization Program using Shor's Algorithm
' FINDPRIM.EXE is at url:http://www.nicvroom.be/findprim.zip
' Input should be two prime numbers,
' but you can try any combination.
'
' Revision 0 Date 9 December 2002
' Revision 1 Date 14 December 2002
' Revision 2 Date 17 December 2002
' Revision 3 Date 20 December 2002 Optimize Standard factorization
' Revision 4 Date 22 December 2002 Increased Array sizes
' Revision 5 Date 1 January 2003 Enter x added in step 1
' Revision 6 Date 5 March 2003 Step 4 improved with Fast Fourier
' Revision 7 Date 10 March 2003 Added Fermat's Sieve algorithm
' Revision 8 Date 13 March 2003 Step 6 improved with convergents calc.
' Revision 8.1 Date 17 March 2003 Improved displays
' Revision 8.2 Date 25 March 2003 Classical Discrete FFT added
' Revision 8.3 Date 1 April 2003 Subroutine CALCXM added
' Revision 8.4 Date 7 April 2003 Optimise p(c) calculation
' Revision 8.5 Date 14 April 2003 SIN COS functionality combined together
' Revision 8.6 Date 9 May 2003 loop count added step 5, 6: sub XAMODN
CONST Euler = 1: CONST nnxn = 11000:
CONST nsxn = 1200
DIM SHARED nxn(nnxn) AS SINGLE 'number of xn
DIM SHARED sxn(nsxn) AS SINGLE 'save xn
DIM SHARED maxlevel
DIM SHARED c, cold, tot, totold, io, r ' Subroutine HTML
DIM SHARED func$, limit ' COS / SIN functionality
DIM SHARED n, x ' CALCXN
RANDOMIZE (TIMER)
io = 0 ' 0 = No disc io 1 = disc io 2 = html io
filenm$ = "shorss.txt" ' file name when io = 1
pi = ATN(1) * 4
CLS
PRINT "Factorization Program using Shor's Algorithm "
PRINT "http://everest.uwyo.edu/~moorhouse/quantum/talk2.pdf"
IF io > 0 THEN COLOR 12: PRINT "Data save to file "; filenm$; " "; io
a1:
COLOR 7
DO
INPUT "Enter prime #1 ", pm1
INPUT "Enter prime #2 ", pm2
n = pm1 * pm2
IF n > nnxn THEN PRINT "Number "; n; "too large to factor "; nnxn
LOOP UNTIL n <> nnxn OR n = 0
IF n = 0 THEN END
COLOR 15
max1 = 2 * n * n ' 2n^2
max2 = 3 * n * n ' 3n^2
PRINT "prime #1 "; pm1; " prime #2 "; pm2;
COLOR 14: PRINT " factor "; n;
COLOR 15: PRINT " min "; max1; "max "; max2
x = 0
q = 1
DO
x = x + 1
q = q * 2
LOOP UNTIL q > max1
PRINT "Step 1, Initial state"
PRINT "x= "; x; " q= "; q
PRINT "Option Enter x"; x;
INPUT lx
IF lx > 0 THEN
x = lx
q = 2 ^ x
PRINT "x= "; x; " q= "; q
END IF
PRINT "Step 2, Apply Modular Exponentiation"
FOR i = 0 TO nnxn: nxn(i) = 0: NEXT i
FOR i = 0 TO nsxn: sxn(i) = 0: NEXT i
IF io > 0 THEN OPEN filenm$ FOR APPEND AS #1
' Limit loop count for "practical" reasons: execution time
qmax = 2 ^ 16
'IF q > 2 ^ 23 THEN
' q = 2 ^ 23
' COLOR 12: PRINT "Reduced x= "; x; " q= "; q
'END IF
q1 = q
c1 = 0 ' Display PRINT control
match = 0
type1 = -1 ' 1 2 3 or 4
ssign = -1
FOR a = 0 TO q - 1
IF a = 0 THEN
xn = 1
ELSE
xn = xn * x ' x^a+1 = x * x^a
xn = xn MOD n ' x^a = x^a MOD n
END IF
IF xn < nnxn THEN nxn(xn) = nxn(xn) + 1 ' number of xn = number of xn + 1
IF a > 0 AND c1 = 0 AND ssign = -1 THEN
' calculate stop sign
jmax = a - 1
' IF jmax > nsxn THEN jmax = nsxn
IF jmax > 10 THEN jmax = 10
FOR j = 0 TO jmax
IF sxn(j) = xn OR xn = 0 THEN ssign = xn: r = a - j: aa = j ' r = periodicity
NEXT j
END IF
IF a > 0 AND c1 = 0 AND xn = ssign THEN
COLOR 14: PRINT USING "####"; a; : COLOR 13: PRINT USING "####"; xn;
IF a > 50 THEN c1 = 1
ELSE
IF c1 = 0 AND (xn > 0 AND a < 100 OR a < 10) THEN
IF xn < 10000 THEN COLOR 14: PRINT USING "####"; a; : COLOR 15: PRINT USING "####"; xn;
END IF
END IF
IF match = 0 THEN
lc = 0
IF a <= nsxn THEN
sxn(a) = xn ' save xn at position i
ELSE
IF a = 2 * nsxn THEN
PRINT
PRINT "Array SXN too small "; nsxn
END IF
END IF
IF a > 0 THEN
IF xn = ssign THEN
IF xn = 1 THEN ' stop sign = 1
IF a MOD 2 = 0 THEN ' odd number of values
IF a = n - 1 THEN
type1 = 6
PRINT
PRINT n; "= Prime number T6"
GOTO a20
ELSE
PRINT
CALCXN a / 2, y
lc = 0
p = GCD(y + 1, n, lc)
q = GCD(y - 1, n, lc)
r = a ' periodicity
IF p <> n THEN
type1 = 1
PRINT "T1 Shor's Algorithm: Odd number of Numbers ";
COLOR 15: PRINT "Periodicity r = "; r
IF io > 0 THEN PRINT #1, "T1 Shor's Algorithm: Odd number of Numbers Periodicity: r = "; r
ELSE
type1 = 4
COLOR 12:
PRINT "T4 Error! Odd number of Numbers ";
COLOR 15: PRINT "Periodicity r = "; r
IF io > 0 THEN PRINT #1, "Error! Odd number of Numbers Periodicity: r = "; r
END IF
COLOR 15: PRINT "GCD ";
COLOR 13: PRINT "a ="; a / 2; " y ="; y; " p ="; p; "q ="; q;
COLOR 7: PRINT TAB(63); "loop count ="; lc
IF io > 0 THEN
PRINT #1, "GCD "; "a ="; a / 2; " y ="; y; " p ="; p; "q ="; q
END IF
END IF
match = a
IF (pm1 <> 1 AND pm2 <> 1) THEN GOTO a20
ELSE
type1 = 3
COLOR 12
r = a
PRINT : PRINT "T3 Error! Even number of Numbers ";
COLOR 15: PRINT "Periodicity r = "; r
a = INT(r / 2)
y = sxn(a)
lc = 0
p = GCD(y, n, lc)
q = n / p
COLOR 15: PRINT "GCD ";
COLOR 13: PRINT "a ="; a; " y ="; y; " p ="; p; "q ="; q;
COLOR 7: PRINT TAB(63); "loop count ="; lc
lc = 0
a = INT(r / 2) + 1
y = sxn(a)
p = GCD(y, n, lc)
q = n / p
COLOR 15: PRINT "GCD ";
COLOR 13: PRINT "a ="; a; " y ="; y; " p ="; p; "q ="; q;
COLOR 7: PRINT TAB(63); "loop count ="; lc
match = a
IF (pm1 <> 1 AND pm2 <> 1) THEN GOTO a20
END IF
ELSE ' Stop sign <> 1
match = a ' Error - no 1 found
COLOR 13
PRINT
IF xn = sxn(1) THEN ' aa=1
type1 = 2
y = sxn(1)
p = GCD(y, n, lc)
q = n / p
PRINT "T2 No stop sign = 1 detected ";
COLOR 15: PRINT "Periodicity r = "; r
COLOR 15: PRINT "GCD ";
COLOR 13: PRINT "a = 1 "; "y ="; y; " p ="; p; "q ="; q;
COLOR 7: PRINT TAB(63); "loop count ="; lc
y = sxn(2)
p = GCD(y, n, lc)
q = n / p
COLOR 15: PRINT "GCD ";
COLOR 13: PRINT "a = 2 "; "y ="; y; " p ="; p; "q ="; q
ELSE
type1 = 5 ' Try 2 and 18 2 and 54
y = sxn(aa)
p = GCD(y, n, lc)
q = n / p
COLOR 12: PRINT "T5 No stop sign = 1 detected ";
COLOR 15: PRINT "Periodicity r = "; r
COLOR 15: PRINT "GCD ";
COLOR 13: PRINT "a ="; aa; "y ="; y; " p ="; p; "q ="; q;
COLOR 7: PRINT TAB(63); "loop count ="; lc
END IF
IF io > 0 THEN
PRINT #1, "T2 No stop sign =1 detected Periodicity r ="; r
PRINT #1, "GCD "; "a ="; aa; " y ="; y; " p ="; p; "q ="; q
END IF
END IF
IF pm1 <> 1 AND pm2 <> 1 THEN GOTO a20
END IF
END IF
ELSE
IF a > qmax THEN GOTO a3 ' No correct frequency list
END IF
NEXT a
a3:
PRINT
IF type1 = -1 THEN
COLOR 12: PRINT "Error": GOTO a20
END IF
INPUT "Select any key"; a$
a4:
PRINT "Frequency list"
cnt = 0
FOR i = 0 TO n
IF nxn(i) > 0 THEN
IF nxn(i) > 1 AND q1 > qmax THEN nxn(i) = INT(q1 / r)
cnt = cnt + 1
COLOR 14: PRINT USING "#####"; i; :
IF INT(q1 / r) < 99999 THEN
COLOR 15: PRINT USING "#####"; nxn(i);
ELSE
COLOR 15: PRINT USING "########"; nxn(i);
IF cnt MOD 6 = 0 THEN PRINT
END IF
IF cnt MOD 160 = 0 THEN
PRINT : INPUT "Enter E to Exit"; a$: cnt = 0
IF a$ <> "" THEN GOTO a5
END IF
END IF
NEXT i
a5:
PRINT
IF type1 > 2 OR r < 2 THEN GOTO a20
' Only for type = 1 or 2
PRINT "Step 3, Observe Register 2";
DO
PRINT
COLOR 15: PRINT "Enter a "; : COLOR 14: PRINT "yellow number ";
COLOR 15: PRINT "from above list. Enter = Random "; :
COLOR 14: INPUT ; sel1
IF sel1 = 0 THEN
DO
ipos = INT(RND * r) + 1:
CALCXN ipos, xn
sel1 = xn
LOOP UNTIL sel1 < nnxn
END IF
LOOP UNTIL nxn(sel1) > 0 OR sel1 = 0
COLOR 14
' calculate a0
a0 = -1 ' calculate later
PRINT :
' PRINT sel1; : COLOR 15: PRINT nxn(sel1); : PRINT " a0"; a0
IF sel1 = 0 THEN
IF io > 0 THEN CLOSE #1
GOTO a20
END IF
COLOR 7: PRINT "Display new state a0, a1, a2 "
q = q1
c1 = 0: tab1 = 1
FOR a = 0 TO q - 1
IF a = 0 THEN
xn = 1
ELSE
xn = xn * x
xn = xn MOD n
END IF
IF xn = 0 THEN a = q - 1 ' Terminate loop
IF xn = sel1 THEN
IF a0 = -1 THEN a0 = a
IF c1 < 6 THEN
c1 = c1 + 1: da = a
IF c1 = 1 THEN COLOR 13 ELSE COLOR 14:
PRINT TAB(tab1);
PRINT USING "#######"; a; :
COLOR 15: PRINT USING "#####"; xn;
tab1 = tab1 + 12
IF c1 = 6 THEN
PRINT : tab1 = 1: da = 6 * r + 1: a = a + INT((q - da - a) / r) * r
' IF q > qmax THEN GOTO a10
END IF
ELSE
PRINT TAB(tab1);
IF a < 9999999 THEN
COLOR 14: PRINT USING "#######"; a; :
tab1 = tab1 + 12
ELSE
COLOR 14: PRINT USING "#############"; a; :
tab1 = tab1 + 20
END IF
COLOR 15: PRINT USING "#####"; xn;
END IF
END IF
NEXT a
a10:
PRINT
q = q1
' Step 4
COLOR 15: PRINT "Step 4, Discrete Fourier Transform of Register 1"
IF q1 > qmax THEN nxn(sel1) = INT(q1 / r)
COLOR 7: PRINT "SQR of "; q * nxn(sel1); " q = "; q; " M = "; nxn(sel1); " r"; r;
COLOR 14: PRINT " Register 2 = "; sel1;
COLOR 7: PRINT " a0"; a0
COLOR 15: INPUT "All FFT state values ? A = All; S = short list; Enter = skip "; a$
type2 = -1 ' skip
IF UCASE$(a$) = "A" THEN type2 = 1
IF UCASE$(a$) = "S" THEN type2 = 0
ctot = 0: cold = -1: maxtot = 0
func$ = "R"
IF io > 0 THEN
PRINT #1, "n="; n; "x="; x; " q="; q; " M="; nxn(sel1); " r"; r; " Register 2= "; sel1; " a0"; a0
END IF
IF type2 < 0 THEN GOTO find0
IF func$ = "" THEN
COLOR 7: INPUT "COS functionality ? C = Cos; S = Sin(SIN functionality); Enter = Yes "; func$
IF UCASE$(func$) = "S" THEN
func$ = "SIN"
ELSE
IF UCASE$(func$) = "R" THEN func$ = "R" ELSE func$ = "COS"
END IF
END IF
IF io > 0 THEN HTML 1
FOR c = 0 TO q1 - 1
tot = 0: Retot = 0: Imtot = 0
w0 = 2 * pi * a0 * c / q1
IF q1 < 10000 THEN
FOR d = 0 TO nxn(sel1) - 1
w = 2 * pi * r * c * d / q1
w = w + w0
Retot = Retot + COS(w)
Imtot = Imtot + SIN(w)
NEXT d
IF func$ = "COS" THEN
tot = Retot
ELSE
IF func$ = "SIN" THEN tot = Imtot ELSE tot = SQR(Retot * Retot + Imtot * Imtot)
END IF
ELSE
' instantaneous calculation
w = 2 * pi * r * c / q1
sw = nxn(sel1) * w
' IF func$ = "COS" THEN dw = 0 ELSE dw = -pi / 2
IF INT(r * c / q1) = r * c / q1 THEN x2 = nxn(sel1) / 2: sw = 0: ELSE x2 = SIN(sw / 2) / SIN(w)
dw = 0
x3 = x2 * COS(sw / 2 + w0 + dw)
x3 = x3 + x2 * COS(sw / 2 - w + w0 + dw)
dw = -pi / 2
x4 = x2 * COS(sw / 2 + w0 + dw)
x4 = x4 + x2 * COS(sw / 2 - w + w0 + dw)
tot = SQR(x3 * x3 + x4 * x4)
END IF
' Summary listing
IF type2 = 0 THEN
IF ABS(tot) > SQR(nxn(sel1)) THEN
IF c > cold + 2 AND cold <> -1 THEN
PRINT : tab1 = 1
ELSE
IF tot * totold < 0 THEN PRINT : tab1 = 1 ' +++ ---
END IF
IF c = q1 / 2 THEN COLOR 12 ELSE COLOR 14
IF q1 <= 8192 THEN
tab1 = tab1 + 13
IF tab1 > 80 THEN PRINT : tab1 = 1
PRINT USING "######"; c;
COLOR 15: PRINT USING "####.##"; tot;
ELSE
tab1 = tab1 + 20
IF tab1 > 80 THEN PRINT : tab1 = 1
PRINT USING "##########"; c;
COLOR 15: PRINT USING "#######.##"; tot;
END IF
IF io > 0 THEN HTML 2
cold = c: totold = tot
END IF
ELSE
IF ABS(tot) > maxtot AND ABS(tot) > SQR(nxn(sel1)) THEN
maxtot = ABS(tot)
ctot = 1
END IF
IF ctot > 0 THEN ctot = ctot + 1
IF c = q1 / 2 THEN COLOR 12 ELSE COLOR 14
PRINT USING "######"; c; :
COLOR 15: PRINT USING "######.##"; tot;
IF io > 0 THEN HTML 4
IF (c + 1) MOD 5 = 0 THEN
PRINT
IF ctot > 50 THEN
INPUT "Select E to End "; a$:
ctot = 0: maxtot = 0
' IF (c + 1) MOD 100 = 0 THEN INPUT "Select E to End "; a$:
IF UCASE$(a$) <> "" THEN GOTO find0
END IF
END IF
END IF
a$ = INKEY$: IF a$ <> "" THEN PRINT : GOTO find0
NEXT c
IF io > 0 THEN HTML 3
PRINT
find0:
INPUT "All probability values ? A = All; S = short list; Enter = skip "; a$
type2 = -1 ' skip
IF UCASE$(a$) = "A" THEN type2 = 1
IF UCASE$(a$) = "S" THEN type2 = 0
ctot = 0: cold = -1: maxtot = 0
IF type2 < 0 THEN GOTO find2
IF func$ = "" THEN
COLOR 7: INPUT "COS functionality ? C = Cos; S = Sin (SIN functionality); Enter = Yes "; func$
IF UCASE$(func$) = "S" THEN
func$ = "SIN"
ELSE
IF UCASE$(func$) = "R" THEN func$ = "R" ELSE func$ = "COS"
END IF
END IF
limit = 1 / r
IF io > 0 THEN HTML 1
grtot = 0
FOR c = 0 TO q1 - 1
' FOR c = 0 TO 20
tot = 0: Retot = 0: Imtot = 0
w0 = 0
IF q1 < 10000 THEN
FOR d = 0 TO nxn(sel1) - 1
w = 2 * pi * r * c * d / q1
w = w + w0
Retot = Retot + COS(w)
Imtot = Imtot + SIN(w)
NEXT d
IF func$ = "COS" THEN
tot = Retot
ELSE
IF func$ = "SIN" THEN tot = Imtot ELSE tot = SQR(Retot * Retot + Imtot * Imtot)
END IF
ELSE
' instantaneous calculation
w = 2 * pi * r * c / q1
sw = nxn(sel1) * w
IF func$ = "COS" THEN dw = 0 ELSE dw = -pi / 2
IF INT(r * c / q1) = r * c / q1 THEN x2 = nxn(sel1) / 2: sw = 0: ELSE x2 = SIN(sw / 2) / SIN(w)
dw = 0
x3 = x2 * COS(sw / 2 + w0 + dw)
x3 = x3 + x2 * COS(sw / 2 - w + w0 + dw)
dw = -pi / 2
x4 = x2 * COS(sw / 2 + w0 + dw)
x4 = x4 + x2 * COS(sw / 2 - w + w0 + dw)
tot = SQR(x3 * x3 + x4 * x4)
END IF
tot = tot * tot / (q * nxn(sel1))
tot = tot * 100
grtot = grtot + tot
' Summary listing
IF type2 = 0 THEN
IF ABS(tot) > limit THEN
IF c > cold + 2 AND cold <> -1 THEN
PRINT : tab1 = 1
ELSE
IF tot * totold < 0 THEN PRINT : tab1 = 1 ' +++ ---
END IF
IF c = q1 / 2 THEN COLOR 12 ELSE COLOR 14
IF r < 100 THEN
tab1 = tab1 + 14
IF tab1 > 80 THEN PRINT : tab1 = 1
PRINT USING "######"; c;
COLOR 15: PRINT USING "####.###"; tot;
ELSE
tab1 = tab1 + 16
IF tab1 > 80 THEN PRINT : tab1 = 1
PRINT USING "########"; c;
COLOR 15: PRINT USING "##.#####"; tot;
END IF
IF io > 0 THEN HTML 5
cold = c: totold = ptot
END IF
ELSE
IF ABS(tot) > maxtot AND ABS(tot) > limit THEN
maxtot = ABS(tot)
ctot = 1
END IF
IF ctot > 0 THEN ctot = ctot + 1
IF tot > limit THEN COLOR 12 ELSE COLOR 14
PRINT USING "######"; c; :
COLOR 15: PRINT USING "##.#####"; tot;
IF io > 0 THEN HTML 6
IF (c + 1) MOD 5 = 0 THEN
PRINT
IF ctot > 50 THEN
INPUT "Select E to End "; a$:
ctot = 0: maxtot = 0
' IF (c + 1) MOD 100 = 0 THEN INPUT "Select E to End "; a$:
IF UCASE$(a$) <> "" THEN GOTO find1
END IF
END IF
END IF
a$ = INKEY$: IF a$ <> "" THEN PRINT : GOTO find1
NEXT c
PRINT
find1:
IF io > 0 THEN HTML 3
COLOR 7: PRINT USING "Grand probability total = ###.###"; grtot
IF io > 0 THEN PRINT #1, USING "Grand probability total = ###.###"; grtot
IF r > 40 THEN INPUT a$
find2:
COLOR 15: PRINT "Most probably measured values ";
COLOR 14: PRINT "in yellow"
FOR d = 0 TO r
COLOR 14
IF d = r OR d = 0 THEN COLOR 15
IF q1 < 99999 THEN
PRINT USING "######.##"; q1 * d / r;
IF (d + 1) MOD 8 = 0 THEN PRINT
ELSE
IF q1 < 99999999 THEN
PRINT USING "#########"; q1 * d / r;
IF (d + 1) MOD 8 = 0 THEN PRINT
ELSE
PRINT USING "############"; q1 * d / r;
IF (d + 1) MOD 6 = 0 THEN PRINT
END IF
END IF
IF io > 0 THEN
PRINT #1, USING "######.##"; q1 * d / r;
IF (d + 1) MOD 8 = 0 THEN PRINT #1, ""
END IF
IF d MOD 150 = 0 AND d <> 0 THEN
COLOR 15: PRINT : INPUT "Select E to End "; a$:
IF UCASE$(a$) <> "" THEN GOTO find4
END IF
NEXT d
IF io > 0 THEN PRINT #1, ""
PRINT
find4:
PRINT "Step 5, Measure Register 1 "
n2 = 0: a$ = "": q = q1
DO
PRINT "Enter Measured value. Enter = Random ";
COLOR 14: INPUT a$
IF a$ <> "" THEN n2 = INT(VAL(a$))
LOOP UNTIL n2 >= 0 AND n2 < q
IF a$ = "" OR n2 = 0 THEN
d = (INT(RND * r) + 1):
IF d = r THEN d = r - 1
n2 = INT(q1 * d / r + .5)
END IF
find5:
COLOR 15: PRINT "Step 6, Continued Fraction Convergent of ";
COLOR 14: PRINT n2
COLOR 7
q = q1: r1 = 0
temp = n2
top = q
a = 0: i = 0: lc = 0
PRINT 0; "a ="; a; TAB(17); "p = "; temp; "/"; top; ""
DO
i = i + 1
bottom = top
top = temp
a = INT(bottom / top)
temp = bottom - a * top
nxn(i) = a
' Calculate convergent
FOR j = i TO 1 STEP -1
IF j = i THEN
n1 = 1: n2 = nxn(j)
ELSE
n3 = n2 * nxn(j) + n1
n1 = n2
n2 = n3
END IF
NEXT j
IF n2 < n THEN r1 = n2: lc = lc + 1
PRINT i; "a ="; a; TAB(17); "p = "; temp; "/"; top; TAB(46); "Convergent "; n1; "/"; n2
LOOP UNTIL temp = 0
IF r1 = 0 THEN
COLOR 12: PRINT " Error - no periodicity found - r1 = "; r1
ELSE
IF r1 > r THEN
COLOR 12: PRINT " Error - periodicity too large - r1 = "; r1
ELSE
IF r MOD r1 <> 0 THEN
COLOR 12: PRINT " Error - periodicity r1 "; r1; " no multiplicity of r"; r
ELSE
COLOR 15: PRINT "Possible values for periodicity r are multiples of r ="; r1; :
COLOR 7: PRINT TAB(63); "loop count ="; lc
IF q < qmax THEN q = q1 ELSE q = qmax
FOR a = 0 TO q - 1
IF a = 0 THEN
xn = 1
ELSE
xn = xn * x
xn = xn MOD n
END IF
IF a > 0 AND a <= r THEN
IF a MOD r1 = 0 THEN COLOR 14: PRINT USING "#####"; a; : COLOR 15: PRINT USING "#####"; xn;
END IF
NEXT a
XAMODN x, r, n, xsave, lc
' PRINT "x"; x, "a"; r; "n"; n; "xsave"; xsave, lc
COLOR 7: PRINT TAB(63); "loop count ="; lc
IF type1 = 1 THEN ' Type 1
a = r
' y = sxn(a / 2)
CALCXN a / 2, y
lc = 0
p = GCD(y + 1, n, lc)
q = GCD(y - 1, n, lc)
ELSE ' Type 2 - No stop sign
a = 2
y = sxn(1)
lc = 0
p = GCD(y, n, lc)
q = n / p
END IF
COLOR 15: PRINT "GCD n ="; n;
COLOR 13: PRINT "a ="; a / 2; " y ="; y; " p ="; p; "q ="; q;
COLOR 7: PRINT TAB(63); "loop count ="; lc
END IF
END IF
END IF
COLOR 15
PRINT "Repeat Step 5 and 6 ? Y = Yes N = No Enter = No ";
COLOR 14: INPUT a$
IF UCASE$(a$) = "Y" THEN GOTO find2
IF UCASE$(a$) = "N" OR a$ = "" GOTO a20
n2 = VAL(a$)
IF n2 > 0 THEN GOTO find5
a20:
COLOR 12: PRINT "This terminates Shor's Algorithm simulation"
COLOR 14
PRINT "Classical factorization of "; n
COLOR 15
n1 = 1 ' Start from number 2
level = 0
lc = 0 ' Loop count
IF n MOD 2 = 0 THEN
FACTOR level, lc, n, n1
ELSE
FACTOR1 level, lc, n, n1
END IF
IF maxlevel = 1 THEN
PRINT "Prime number ";
ELSE
IF maxlevel = 2 THEN
PRINT "Eureka ";
ELSE
COLOR 12: PRINT "Error ";
END IF
END IF
FOR i = 1 TO maxlevel: PRINT nxn(i); ""; : NEXT i
PRINT " Loop count = "; lc
COLOR 14
PRINT "Fermat's Sieve factorization of "; n
FACTOR2 n, lc, n1, n2
COLOR 15
SELECT CASE n1
CASE IS = 0
COLOR 12: PRINT "Error "; n1; n2;
CASE IS = 1
PRINT "Prime number"; n2;
CASE ELSE
PRINT "Eureka "; n1; n2;
END SELECT
PRINT " Loop count = "; lc
IF io > 0 THEN CLOSE #1
IF type1 > 2 OR (pm1 <> 1 AND pm2 <> 1) GOTO a30
type3 = 0
COLOR 7: PRINT "Classical Discrete FFT ? Y or N ";
INPUT a$
IF (VAL(a$) >= 1 AND VAL(a$) <= 3) THEN
type3 = VAL(a$)
ELSE
IF UCASE$(a$) <> "Y" THEN GOTO a30
type3 = 3
END IF
COLOR 15: PRINT "Classical Discrete Fast Fourier Transform "; type3
off1 = 0
IF type1 = 1 THEN off1 = 0
IF type1 = 2 THEN off1 = 1
mulf = INT(100 / r) + 1
IF mulf > 3 THEN mulf = 3
FOR w = 0 TO mulf * r
tot = 0
FOR c = 0 TO r - 1
c1 = c MOD r + off1
w1 = (c + w) MOD r + off1
IF type3 = 1 THEN tot = tot + sxn(c1) * sxn(w1)
IF type3 = 2 THEN
IF sxn(c1) = sxn(off1) THEN valx = 1 ELSE valx = 0
IF sxn(w1) = sxn(off1) THEN valy = 1 ELSE valy = 0
tot = tot + valx * valy
END IF
IF type3 = 3 THEN tot = tot + COS(sxn(c1) - sxn(w1))
NEXT c
COLOR 14: PRINT USING "####"; w; :
COLOR 15: IF w MOD r = 0 THEN COLOR 12
PRINT USING "######"; tot;
NEXT w
PRINT
COLOR 12: PRINT "Peridicity r = "; r
a30: GOTO a1
END
SUB CALCXN (ipos, xn)
IF ipos <= nsxn THEN
xn = sxn(ipos)
ELSE
xn = sxn(nsxn)
FOR a = nsxn + 1 TO ipos
xn = xn * x
xn = xn MOD n
NEXT a
END IF
' PRINT "CALCXN"; ipos; xn
END SUB
SUB FACTOR (level, lc, n, n1)
' Even and Odd numbers
' n is the number to factor
' n1 is the start value -1
' lc = loop count
' level = level i.e. how often called
level = level + 1
maxlevel = level
DO
n1 = n1 + 1
lc = lc + 1 ' loop count
IF n MOD n1 = 0 THEN ' test if n1 = prime number
n2 = n / n1 ' New value to factor
nxn(level) = n1 ' save number
n1 = n1 - 1 ' Initialize n1 with previous value.
IF n2 <> 1 THEN ' Stop when New value is 1.
FACTOR level, lc, n2, n1
END IF
level = level - 1
EXIT SUB
END IF
LOOP UNTIL n1 * n1 > n
nxn(level) = n
level = level - 1
END SUB
SUB FACTOR1 (level, lc, n, n1)
' Odd numbers
' n is the number to factor
' n1 is the start value -1
' lc = loop count
' level = level i.e. how often called
level = level + 1
maxlevel = level
DO
n1 = n1 + 2
lc = lc + 1 ' loop count
IF n MOD n1 = 0 THEN ' test if n1 = prime number
n2 = n / n1 ' New value to factor
nxn(level) = n1 ' save number
n1 = n1 - 2 ' Initialize n1 with previous value.
IF n2 <> 1 THEN ' Stop when New value is 1.
FACTOR1 level, lc, n2, n1
END IF
level = level - 1
EXIT SUB
END IF
LOOP UNTIL n1 * n1 > n
nxn(level) = n
level = level - 1
'
END SUB
SUB FACTOR2 (n, lc, n1, n2)
' Fermat's Sieve factorization of n
' lc = loop count
nn = SQR(n)
nn = INT(nn)
lc = 0
' PRINT lc, n
DO
nn = nn + 1
s = nn * nn - n
s1 = INT(SQR(s))
lc = lc + 1
' PRINT lc, nn, s, s1
LOOP UNTIL s1 * s1 = s OR nn + s1 > n
IF s1 * s1 = s THEN
n1 = nn - s1:
n2 = nn + s1
ELSE
n1 = 0: n2 = n
END IF
END SUB
FUNCTION GCD (a1, b1, lc)
' Greatest Common Denominator.
IF Euler = 1 THEN
a = a1: b = b1
IF b > a THEN c = a: a = b: b = c
GCD1:
IF b = 0 THEN GCD = a: EXIT FUNCTION
lc = lc + 1
res1 = INT(a / b): c = a - res1 * b
a = b: b = c: c = 0
GOTO GCD1
ELSE
a = a1: b = b1
x = 1
match = 0
DO
x = x + 1
res1 = a / x
IF INT(res1) * x = a THEN
res2 = b / x
IF INT(res2) * x = b THEN match = x
res2 = b / res1
IF INT(res2) * res1 = b THEN match = res1
END IF
LOOP UNTIL match > 0 OR (x * x > a AND x * x > b)
PRINT "GCD"; x; a; b
GCD = match
END IF
END FUNCTION
SUB GETPRIEM (a1, a2)
COLOR 7
al = a1
n = 1
a2 = 0
' PRINT a1, a2
DO
n = 1
DO
n = n + 2
y = INT(al / n)
' LOOP UNTIL al MOD n = 0 OR n * n > a1
' PRINT al; "n"; n; "y"; y; n * y; n * n
LOOP UNTIL n * y = al OR n * n > al
' INPUT a$
IF n * n > al THEN a1 = al: a2 = n: EXIT SUB ' PRINT "priem"; a1; n: END
al = al + 2
y = INT(al / 100)
IF y * 100 = a1 THEN PRINT a1
a$ = INKEY$: IF a$ = CHR$(ESC) THEN END
LOOP
END SUB
SUB HTML (id)
SELECT CASE id
CASE IS = 1
IF io = 1 THEN
PRINT func$; " Functionality "; TAB(20); "limit ="; limit
END IF
IF io = 2 THEN
PRINT #1, func$; " Functionality
"
PRINT #1, "limit ="; limit
PRINT #1, ""
PRINT #1, " "
END IF
CASE IS = 2
IF io = 1 THEN
IF c > cold + 5 AND cold <> -1 THEN
PRINT #1, ""
ELSE
IF tot * totold < 0 THEN ' +++ ---
PRINT #1, ""
END IF
END IF
IF q1 <= 8192 THEN
PRINT #1, USING "######"; c;
PRINT #1, USING "####.##"; tot;
ELSE
PRINT #1, USING "##########"; c;
PRINT #1, USING "#######.##"; tot;
END IF
ELSE
IF c > cold + 5 AND cold <> -1 THEN
PRINT #1, "
"
PRINT #1, " "
ELSE
IF tot * totold < 0 THEN ' +++ ---
PRINT #1, "
"
PRINT #1, " "
END IF
END IF
IF q1 <= 8192 THEN
PRINT #1, " | ";
PRINT #1, USING "######"; c;
PRINT #1, " | ";
PRINT #1, "";
PRINT #1, USING "####.##"; tot;
PRINT #1, " | "
ELSE
PRINT #1, " ";
PRINT #1, USING "##########"; c;
PRINT #1, " | ";
PRINT #1, "";
PRINT #1, USING "#######.##"; tot;
PRINT #1, " | "
END IF
END IF
CASE IS = 3
IF io = 1 THEN
PRINT #1, ""
ELSE
PRINT #1, "
"
PRINT #1, "
"; CHR$(13); CHR$(10)
END IF
CASE IS = 4
IF io = 1 THEN
PRINT #1, USING "######"; c; :
PRINT #1, USING "####.##"; tot;
IF (c + 1) MOD 8 = 0 THEN
PRINT #1, ""
END IF
ELSE
PRINT #1, " ";
PRINT #1, USING "#######"; c;
PRINT #1, " | ";
PRINT #1, "";
PRINT #1, USING "####.##"; tot;
PRINT #1, " | "
IF (c + 1) MOD 8 = 0 THEN
PRINT #1, " "
PRINT #1, " "
END IF
END IF
CASE IS = 5
IF io = 1 THEN
IF c > cold + 5 AND cold <> -1 THEN
PRINT #1, ""
ELSE
IF tot * totold < 0 THEN ' +++ ---
PRINT #1, ""
END IF
END IF
IF r < 100 THEN
PRINT #1, USING "######"; c;
PRINT #1, USING "####.###"; tot;
ELSE
PRINT #1, USING "########"; c;
PRINT #1, USING "##.#####"; tot;
END IF
ELSE
IF c > cold + 5 AND cold <> -1 THEN
PRINT #1, "
"
PRINT #1, " "
ELSE
IF tot * totold < 0 THEN ' +++ ---
PRINT #1, "
"
PRINT #1, " "
END IF
END IF
IF r < 100 THEN
PRINT #1, " | ";
PRINT #1, USING "######"; c;
PRINT #1, " | ";
PRINT #1, "";
PRINT #1, USING "####.###"; tot;
PRINT #1, " | "
ELSE
PRINT #1, " ";
PRINT #1, USING "########"; c;
PRINT #1, " | ";
PRINT #1, "";
PRINT #1, USING "##.#####"; tot;
PRINT #1, " | "
END IF
END IF
CASE IS = 6
IF io = 1 THEN
PRINT #1, USING "######"; c; :
PRINT #1, USING "##.#####"; tot;
IF (c + 1) MOD 8 = 0 THEN
PRINT #1, ""
END IF
ELSE
PRINT #1, " ";
PRINT #1, USING "#######"; c;
PRINT #1, " | ";
PRINT #1, "";
PRINT #1, USING "##.#####"; tot;
PRINT #1, " | "
IF (c + 1) MOD 8 = 0 THEN
PRINT #1, "
"
PRINT #1, " "
END IF
END IF
END SELECT
END SUB
SUB XAMODN (x, a, n, xsave, lc)
' Calculate x^a mod n
' Concepts from Grigoris Lionis - See sci.physics 8/5/03
DIM save(30), powerx(30)
i = 0: lc = 0
x1 = x: a1 = a
DO
IF i = 0 THEN
x1 = x1 MOD n
power = 1
ELSE
' x1 = x1 * x1 MOD n
ps = x1 * x1
qs = INT(ps / n)
x1 = ps - qs * n
lc = lc + 1
power = power * 2
END IF
save(i) = x1
powerx(i) = power
rest = a1 - power
i = i + 1
' PRINT i, power, x1, rest
LOOP UNTIL rest < power
xsave = x1
IF rest > 0 THEN
DO
i = 0
a1 = rest
DO
x1 = save(i)
power = powerx(i)
rest = a1 - power
i = i + 1
'PRINT i, power, x1, rest
'INPUT a$: IF UCASE$(a$) = "E" THEN END
LOOP UNTIL rest < power
' xsave = xsave * x1: xsave = xsave MOD n:
ps = xsave * x1
qs = INT(ps / n)
xsave = ps - qs * n
lc = lc + 1
LOOP UNTIL rest = 0
END IF
END SUB
Back to my home page Contents of This Document