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://users.pandora.be/nicvroom/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, "" ELSE PRINT #1, " "; PRINT #1, "" END IF END IF CASE IS = 3 IF io = 1 THEN PRINT #1, "" ELSE PRINT #1, "" PRINT #1, "
"; PRINT #1, USING "######"; c; PRINT #1, ""; PRINT #1, USING "####.##"; tot; PRINT #1, ""; PRINT #1, USING "##########"; c; PRINT #1, ""; PRINT #1, USING "#######.##"; tot; 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