Plotting a function in hi-res graphics mode on a Commodore 64

I wrote a program to plot a function in hi-res graphics mode on the C64. It is entirely in BASIC, so it is mind-numbingly slow, but it works, and can be used to plot some fairly interesting functions.

This mode only allows one background color and one foreground color for each 8×8 block of pixels. It would, I suppose, be possible to plot multiple functions in multiple colors, but the intersections would be strange, as each 8×8 block can only have two colors, including the background. An obvious next step is to use the multi-color bitmap mode, but which this would allow four colors in any given 8×8 block, it would reduce the resolution from 320×200 down to 160×200. I will try this soon.

Another obvious improvement will be to move much of the setup stuff (like erasing the graphics ram ($2000-$3fff)) down into machine language, which would speed things up significantly.

The Graphics Book for the Commodore 64 was indispensable.

Here is the plot of f(x) = sin(x) + sin(1.5*x + pi/2)

And here is the code:

 10 rem *** funcplot
 20 rem
 30 pi = 3.141592653
 40 yu = 2.2 : yd = -2.2
 50 xu = 6pi : xd = -6pi
 60 xr = xu - xd : yr = yu - yd
 70 xs = xr/320 : ys = yr/200
 80 def fn f(x) = sin(x) + sin(1.5*x + pi/2)
 100 rem set vic address
 110 v = 53248
 120 rem set graphic ram address
 130 ga = 8192
 140 rem set video ram address
 150 vr = 1024
 160 rem set border to black
 170 poke v+32,0
 500 gosub 20000 : rem turn on hires
 510 gosub 21000 : rem graphic ram area
 520 gosub 22000 : rem set color ram
 530 gosub 23000 : rem clr graphic ram
 1000 rem draw x axis
 1005 y=0 : yc=int(200-(y-yd)*(200/yr)) 
 1010 for xc = 0 to 319 
 1020 gosub 30000 : rem set (xc,yc) 
 1030 next xc 
 1100 rem draw y axis
 1110 x = 0
 1120 for yc = 1 to 200
 1130 xc = int((x-xd)*(320/xr))
 1140 gosub 30000 : rem set (xc,yc)
 1150 next yc
 1300 rem plot function
 1310 for x = xd to xu step xs
 1320 y = fn f(x)
 1330 xc = int((x-xd)*(320/xr))
 1340 yc = int(200-(y-yd)*(200/yr))
 1350 gosub 30000 : rem set (xc,yc)
 1360 next x
 1400 poke 198,0 : wait 198,1
 1450 gosub 24000
 1500 end
 20000 rem turn on hi res graphics
 20010 rem   1. set bits 5/6 of v+17
 20020 rem   2. clr bit 4 of v+22
 20030 poke v+17,peek(v+17) or (11*16)
 20040 poke v+22,peek(v+22) and (255-16)
 20050 return
 21000 rem set graphic ram area
 21010 rem   1. set bit 3 of v+24
 21020 poke v+24, peek(v+24) or 8
 21030 return
 22000 rem set color ram
 22010 rem   1. color ram is 1024-2023
 22020 rem   2. set background 1 - white
 22030 rem   3. set foreground 0 - black
 22040 co = 0*16 + 1
 22050 for i = vr to vr+1000
 22060 poke i,co
 22070 next i
 22080 return
 23000 rem clear graphic ram
 23010 rem   1. graphic ram is ga to
 23020 rem        ga + 8000
 23030 for i = ga to ga +8000
 23040 poke i,0
 23050 next i
 23060 return
 24000 rem turn graphics off
 24010 rem   1. clr bits 5/6 of v+17
 24020 rem   2. clr bit 4 of v+22
 24030 rem   3. clr bit 3 of v+24
 24040 poke v+17,peek(v+17) and (255-96)
 24050 poke v+22,peek(v+22) and (255-16)
 24060 poke v+24,peek(v+24) and (255-8)
 24070 return
 30000 rem set pixel
 30010 ra = 320*int(yc/8)+(yc and 7)
 30020 ba = 8*int(xc/8)
 30030 ma = 2^(7-(xc and 7))
 30040 ad = ga + ra + ba
 30050 poke ad,peek(ad) or  ma
 30060 return
 31000 rem clr pixel
 31010 ra = 320*int(yc/8)+(yc and 7)
 31020 ba = 8*int(xc/8)
 31030 ma = 255-2^(7-(xc and 7))
 31040 ad = ga + ra + ba
 31050 poke ad,peek(ad) or  ma
 31060 return
 

Exploring the Goldbach Conjecture using a Commodore 64

The Goldbach Conjecture states that every even integer greater than two can be expressed as a sum of two primes.

I programmed my C64 in Basic to ask the question “How many ways can a particular even integer be expressed as the sum of two primes?” So, for instance, the number 10 can be expressed as 3+7 or as 5+5, so two ways. The number 12 can be expressed as 5+7, so only one way. The number 22 can be expressed as 11+11 or as 19+3 or as 17+5, so three ways.

I left the program running for twelve hours, to calculate the Goldbach partitions of all even integers from six to one thousand. I’m very sure that there are much more efficient routines, which would run faster, but it is at least clear that it can be done.

Note that the Goldbach conjecture is not proven. The Wikipedia article here states that the conjecture has been tested by computers for numbers up to 4 x 10^18. But I only went as high as 1000.

Here is a plot of my results (not on the C64, but rather using gnuplot on Windows.

And here is the code:

 50 qm = 1000
 100 rem get primes up to 1000
 200 dim p%(400)
 220 p%(1) = 2
 230 pc = 1
 240 n = 3
 300 for i = 1 to pc
 310 w = n/p%(i)
 320 if w = int(w) then 400
 330 next i
 340 pc = pc + 1
 350 p%(pc) = n
 355 print n
 400 for i = 1 to 1: next i: n = n+1
 410 if n < qm then goto 300
 1000 rem get golbach numbers
 1010 dim gb%(501)
 1020 rem gb%(1) = 0 and gb%(2) = 0
 1030 rem because 2 and 4 have no
 1040 rem goldbach partitions
 1050 gb%(1) = 0 : gb%(2) = 0
 1060 q2 = int(qm/2)
 1100 for m = 3 to q2
 1110 n = m*2
 1120 ct = 0
 1130 for i = 1 to pc
 1140 if p%(i) > n then goto 1200
 1150 pr = n - p%(i)
 1160 pf = 0
 1162 for j = 1 to pc
 1164 if p%(j) = pr then pf = 1:goto 1195
 1166 if p%(j) > pr then goto 1200
 1168 next j
 1195 if pf = 1 then ct = ct + 1
 1200 for j = 1 to 1:nextj:next i
 1210 if (ct/2) <> int(ct/2) then ct = ct + 1
 1220 ct = int(ct/2)
 1250 print n,ct
 1260 gb%(m) = ct
 1300 next m
 2000 open 2,8,2,"@:goldlist,s,w"
 2010 for q = 1 to int(qm/2)
 2020 print#2,q*2,gb%(q)
 2030 next q
 2040 close 2
 

Commodore 64: Using the Newton/Raphson method and Halley’s method to find the zeros of a function

Click for good descriptions of Newton’s method and Halley’s method.

I decided to write my C64 Basic code for two different functions:

  • f(x) = x*x*x – 2*x – 5
  • f(x) = exp(x) – x*x

Here are the plots for these functions, with the roots indicated by a blue dot:

These methods require choosing a first guess for the root. In both cases I chose x = 2 as a first guess.

Here is the code:

 10 REM ===============================
 20 REM
 30 REM FIND THE ZEROS OF A FUNCTION
 35 REM USING NEWTON/RAPHSON METHOD
 40 REM AND HALLEY'S METHOD
 60 REM
 70 REM ===============================
 100 PRINT CHR$(147)
 200 DEF FN F(A) = A*A*A - 2*A - 5 
 210 DEF FN F1(A) = 3*A*A - 2 
 220 DEF FN F2(A) = 6*A
 230 E$ = "XXX - 2X - 5" 
 235 PRINT "--------------------------------------" 
 240 GOSUB 10000 
 245 PRINT "--------------------------------------" 
 250 GOSUB 20000 
 300 DEF FN F(A) = EXP(A) - A*A
 310 DEF FN F1(A) = EXP(A) - 2*A 
 320 DEF FN F2(A) = EXP(A) - 2  
 330 E$ = "EXP(X) - X*X"
 335 PRINT "--------------------------------------"
 340 GOSUB 10000
 345 PRINT "--------------------------------------"
 350 GOSUB 20000
 1000 END
 10000 REM NEWTON/RAPHSON
 10002 X = 2
 10005 S = 0 : DI = 10 : TL = 0.000001
 10010 IF DI < TL THEN GOTO 10100
 10020 G = FN F(X) : G1 = FN F1(X)
 10030 X1 = X - G/G1
 10040 DI = ABS(X - X1)
 10050 X = X1
 10060 S = S + 1
 10070 GOTO 10010
 10100 PRINT "NEWTON'S METHOD FOR: ";E$
 10110 PRINT "ROOT = ";X
 10120 PRINT "F(ROOT) = "; FN F(X)
 10130 PRINT "FOUND IN ";S;" STEPS"
 10200 RETURN
 20000 REM HALLEY
 20002 X = 2
 20005 S = 0 : DI = 10 : TL = 0.000001
 20010 IF DI < TL THEN GOTO 20100
 20020 G = FN F(X) : G1 = FN F1(X)
 20025 G2 = FN F2(X)
 20030 NU = 2 * G * G1
 20040 DE = 2 * G1 * G1 - G * G2
 20050 X1 = X - NU/DE
 20060 DI = ABS(X - X1)
 20070 X = X1
 20080 S = S + 1
 20090 GOTO 20010
 20100 PRINT "HALLEY'S METHOD FOR ";E$
 20110 PRINT "ROOT = ";X
 20120 PRINT "F(ROOT) = "; FN F(X)
 20130 PRINT "FOUND IN ";S;" STEPS"
 20200 RETURN

Here is the screen shot: