    1REM  >  $.donmcd.calc.factors.lgefactor.IsPriMe2kp,  2.05.98.
    2MODE 8  : COLOUR 2 : REM 4 COLOUR YELLOW
    3FOR d = 1 TO 2
    4PRINT "BBC Basic 64 Progm  $.donmcd.calc.factors.lgefactor.IsPriMe2kp, 10/3/95, 5.05.98 "
    5   PRINT ' "**  RECOMMENDed  BBC BASIC 64.    *********************"
    6NEXT d
    7COLOUR 3 : REM 4 COLOUR WHITE.
    8PRINT ' " by Donald S  McDonald, 63/5 Hutchison Rd, NEWTOWN, Wellington 2"'" New Zealand.   Phone 64(NZ) +4(WN-) + 389-6820."
    9
   10PRINT ''"Solve for (odd) factors, p <= 2^16.5 = 92,681."'" = 200,000,000? of Large Numbers > 2E9 of special type .."'
   11
   12COLOUR 2
   13PRINT"   ABS(    C * a^b -i  )   <=   n  MOD p,     Progm ""IsPriMe2kp"""
   14PRINT"   multiplier * (base ^ index) + offset  <=  neighbourhood."
   15
   16REM @% = &2090A   :
   17    now = TIME  : REPORT
   18COLOUR 3
   19
   20PRINT '"Would you like to see some samples, examples?  Yy ( N) ? "
   21 x$ = GET$
   22 IF INSTR("Yy", x$) THEN PROCexamps
   23
   24c% = 0 : a% = 0 : b% = 0 : i% = 0 : n% = 0 :  REM not Resident variables.
   25PROCinput
   26now = TIME  :  @% = &90A
   27PROCbinary (b%)
   28ON ERROR REPORT :  PROCerror1
   29  @% = &90A
   30PROCmult
   31
   32 PROCerror2
   33END
   34
   35DEF PROCinput
   36COLOUR 2
   37PRINT ' "PROC input."'
   38COLOUR 3
   39INPUT LINE "enter text, formula.  Analysis of c, a, b, i, n  later : "; F$
   40
   41PRINT"      multiplier * (base ^ index) + offset  <=  neighbourhood."
   42PRINT"        2E9          2E9    2E9      2E9           2^5.  "'
   43INPUT "enter base ( expression a$ , max= 2E9, often 2 , 10, 2^10, 1E3. ) "; a$
   44a% = EVAL(a$)
   45PRINT "a$ = ";a$ " = " a%
   46REPEAT
   47   INPUT '"enter index (power b%, [default 1]) "; b$
   48   IF b$ = "" THEN b% = 1   ELSE  b% = EVAL(b$)
   49   IF b% < 2 THEN PRINT CHR$7 ' "Recomm. also use of Programs ""SERICalc4, PrimeCabi, a2p999"'"  //  Factors"" for simple primes"'"   and factors, to Integer Range 2^31,  b% <2."
   50   UNTIL b% >=1
   51INPUT "enter multiplier, c% (default = 1) "; c$
   52IF c$ = "" THEN c% = 1  ELSE c% = EVAL(c$)
   53PRINT; c%
   54INPUT "enter offset, +/- i% = EVAL i$ (default = -1) "; i$
   55IF i$ = "" THEN i% = -1 ELSE i% = EVAL(i$)
   56PRINT i$ " = "; i%
   57REPEAT
   58     INPUT "enter neighbourhood (sm. latitude) <= n%,  (default = 0 ) "; n$
   59     IF n$ = "" THEN n% = 0 ELSE n% = EVAL(n$)
   60     IF n% < 0 THEN PRINT "neighbourhood must not be negative, Re-enter."
   61     UNTIL n% >= 0
   62PRINT '' "find factors p, satisfying ,  Sint =  " F$ ' " ABS( ";c%" * "a%" ^ " b% " + "i% " )  <= "n%"  MOD p."
   63ENDPROC
   64
   65DEF PROCerror1
   66REM  @% = &2090A
   67PRINT ' " find factors p, satisfying Sint = " F$     ' " ABS( ";c%   " * "a%         " ^ " b%        " + "i%         " )  <= "n%     "  MOD p."
   68PRINT "highest factor tested,  P% = "; P%  "    cs. = "(TIME - now)"   at ERL.  =  ";ERL
   69PRINT CHR$7 "BBC Basic 64 Prog $..calc.normal.IsPriMe2kp,  e n d .  "
   70ENDPROC
   71DEF PROCerror2
   72 INPUT "RUN AGAIN  Enter G/go, else <RETurn> end. CLOSE Problem." ; get$
   73 IF get$<>"" AND INSTR("G/go", get$) THEN RUN
   74         CLOSE #0 :  *SPOOL
   75         PRINT "CLOSE  #0, ( ALL),  * SPOOL"
   76 
   77 STOP
   78ENDPROC
   79
   80DEF PROCbinary(b%)
   81PRINT '"PROC binary ( "; b%
   82  PRINT "Hex. "; ~b%
   83  h = 1  :  h2 = 2^h
   84  PRINT "bits h = ";h  " poss. area of improvement 1000^&F ?   2^h = "h2
   85DIM B%(LN(b%) / (h*LN2) +2)
   86B% = 0  :  B%(0) = b%
   87REPEAT
   88        B%(B%+1) = B%(B%) DIV h2
   89        B%(B%) = B%(B%) MOD h2
   90        PRINT; B%(B%+1) " ,  "B%(B%)
   91        B% += 1
   92        UNTIL B%(B%) = 0
   93
   94BB = B%
   95ENDPROC
   96
   97DEF PROCmult
   98
   99PRINT '"PROC mult."
  100T% = SQR( 2^31)
  101
  102
  103  INPUT "enter test divisors:  start$ [Return = 1,  .0 = quit] = "; start$
  104  IF start$ = "" THEN start = 1  ELSE start = EVAL( start$ )
  105  IF start = 0 THEN ENDPROC : start = 1   : REM  ??
  106 PRINT  "start " start$ " = "; start
  107 PRINT "test unique even prime (plus powers 2)  Spc = Yes,  Nn no.?"
  108 IF  GET$ = " " THEN PRINT "Yes.  " : PROCcf(2)  :
  109  INPUT "enter step$, generally even 2,  < Return = 2*b% > " ; step$
  110       IF step$ = ""  THEN step$ = "2*b%"
  111  step = EVAL( step$ )
  112  IF step = 0 THEN step = 2
  113  PRINT "step " step$ " = " ; step
  114  INPUT "enter minimum [Return = 3]  [N.B.** recommend BASIC 64.**] ";minimum$
  115  IF minimum$ = "" THEN  min = 3 ELSE  min = EVAL(minimum$)
  116  IF min < 3 THEN min = 3
  117  PRINT "minimum " minimum$ "  "  min  CHR$7
  118  IF min > start THEN start = start + step * INT(  (min-start-1) / step +1)
  119 now = TIME :
  120FOR P% = start TO 2^30 STEP step  :  REM  To ??
  121IF P% < 51 OR (P% MOD 3)*(P% MOD 5)*(P% MOD 7) > 0 THEN
  122     IF P% < 51 OR (P% MOD 11)*(P% MOD 13) > 0 THEN PROCcf(P%)
  123     ENDIF
  124NEXT P%
  125ENDPROC
  126
  127DEF FNM(X)
  128 IF X = X+1 THEN ERROR X, "exceeds accuracy of machine #//loss of accuracy"
  129 
  130X = X - P% * INT(X / P% +.5)
  131= X
  132:
  133
  134DEF PROCcf(P%)
  135LOCAL  POWER
  136  PRINT ",";  : REM centisec.  reduce output print.
  137       IF POS > 65 THEN PRINT ';P% " cs.= "TIME-now;
  138     A% = FNM(a%)
  139     Sint = A%  :  
  140     I% = FNM(i%)
  141     C% = FNM(c%)
  142     IF b% > 1 THEN
  143             IF BB < 2 THEN PRINT " BB = "; BB "  Sint = " Sint : STOP
  144     FOR  B% = BB-2 TO 0 STEP -1
  145        FOR bit = 1 TO h : Sint = FNM(Sint^2)  : NEXT bit
  146            REM    PRINT Sint "  ";
  147          Sint = FNM(FNM(A% ^ B%(B%)) * Sint) :  REM  PRINT ; Sint ",  ";
  148        NEXT B%
  149        ENDIF
  150
  151Sint = FNM( C% * Sint + I%)
  152IF ABS(Sint) <= n%  THEN
  153         PRINT ' CHR$7 "Sint = "; Sint "  factor = ";P% "  ";
  154         PPRIM = TRUE
  155         IF P% MOD 2 = 0 THEN  PPRIM = FALSE  :  PRINT " = "; 2 " * " P% DIV 2 ;
  156         FOR J% = 3 TO SQR P%
  157                IF P% MOD J% = 0 THEN PPRIM = FALSE : PRINT " = "; J% " * " P% DIV J% ;
  158                NEXT J%
  159         IF P% <= 3 OR PPRIM THEN
  160         PRINT "PRIME. (test dup. factors).  ";
  161         POWER = P%
  162         REPEAT
  163                POWER = POWER * P%
  164                IF POWER < T% THEN PROCcf(INT(POWER))
  165                UNTIL POWER >= T%  OR ABS(Sint) > n% : REM ??
  166         ENDIF
  167         PRINT  CHR$7 "  contin cs.= "; (TIME - now) "   " GET$
  168ENDIF
  169ENDPROC
  170  
  171DEF PROCexamps
  172PRINT '"Example.  10*1E4^7 - 1.   Factors of REP.UNIT 29 = 11,111,111, ..., 111."'" 9.3191.16763.43037.62003.yes.77843839397.last factor too big."
  173PRINT '"Example. W = 391581*2^216193 -1, Former World's Record largest known Prime"' "   No. Guinn Bk of Records 1991."'"4W + 3 factors 7* 6199 * 92219 * ... etc."
  174PRINT '"Example.  PI = 3 141 592 653 589 793 25-. =3*5*5* 66679* = 131*509* XXX"
  175PRINT "Ex. Bridge, 158 753 389 900.    Ex. RUBIK.  43 252 003 274 489 856 000."
  176PRINT "Ex. Bankcard no.,               Ex. ATM card no. (Both 16 digits)"
  177PRINT'"Ex. 35099 DIVIDES 1011,1213,1415,1617,1819.  (20 digits.)"
  178PRINT "Ex. Proves prime up to 2^33-9 = 8E9.  BASIC 64 c. 1E14?"
  179PRINT "Ex. 1,234,567,891*1E11 + 011,1213,1415  (21 digits.)"
  180PRINT '"Exs. 641 is factor of each following "
  181PRINT "2^32+1  (Fermat),  5^32+1,  1E32-1 (Rep.unit 32?)"
  182PRINT"27! + 1 prime? Yes. = 3*(2*11*13)^2*17*19*23 * ( 2^7*3^4*5^2*7 ) ^ 3 +1."
  183  PRINT "contin."GET$
  184PRINT"Maple V.2 uses a rather weak probabilistic test, try 2,152,302,898,747 ="'" 6763*10627*29947.  WelComBBS.SciMath2, 12/3/95."
  185PRINT"Largest known twin primes..   1,691,232*1001*1e4020+-1,  Dubner, Sci.Math"
  186 PRINT "Cole. 2^67-1, search minimum factor = 193.7 million"
  187 PRINT "70*858433+1  DIVIDES  2^858433 -1. (typo CDROM. M_859433 Is prime.)"
  188 PRINT "134MILL DIVIDES M239.  69,728,591 DIVIDES M6972859"
  189 PRINT "13,945,727 DIVIDES 2^6972863-1.  M-6972593 ISPRIME JUL.1999."
  190ENDPROC