Wilson primes of order n: Difference between revisions

=={{header|GW-BASIC}}==

<lang gwbasic>10 PRINT "n: Wilson primes"

20 PRINT "--------------------"

30 FOR N = 1 TO 11

40 PRINT USING "##";N;

50 FOR P=2 TO 18

60 GOSUB 140

70 IF PT=0 THEN GOTO 100

80 GOSUB 230

90 IF WNPT=1 THEN PRINT P;

100 NEXT P

110 PRINT

120 NEXT N

130 END

140 REM tests if the number P is prime

150 REM result is stored in PT

160 PT = 1

170 IF P=2 THEN RETURN

175 IF P MOD 2 = 0 THEN PT=0:RETURN

180 J=3

190 IF J*J>P THEN RETURN

200 IF P MOD J = 0 THEN PT = 0: RETURN

210 J = J + 2

220 GOTO 190

230 REM tests if the prime P is a Wilson prime of order N

240 REM make sure it actually is prime first

250 REM RESULT is stored in WNPT

260 WNPT=0

270 IF P=2 AND N=2 THEN WNPT = 1: RETURN

280 IF N>P THEN WNPT=0: RETURN

290 PROD# = 1 : P2 = P*P

300 FOR I = 1 TO N-1

310 PROD# = (PROD#*I) : GOSUB 3000

320 NEXT I

330 FOR I = 1 TO P-N

340 PROD# = (PROD#*I) : GOSUB 3000

350 NEXT I

360 PROD# = (P2+PROD#-(-1)^N) : GOSUB 3000

370 IF PROD# = 0 THEN WNPT = 1: RETURN

380 WNPT = 0: RETURN

3000 REM PROD# MOD P2 fails if PROD#>32767 so brew our own modulus function

3010 PROD# = PROD# - INT(PROD#/P2)*P2

3020 RETURN</lang>