DEFDBL A-Z
DIM col AS LONG, i AS LONG

SCREEN _NEWIMAGE(800, 600, 32)
cx = 400: cy = 300

WHILE INKEY$ <> CHR$(27)
  a = -7: b = 10
  CLS

  FOR u = 1 TO 7000
    xr = RND * 12 - 6: xi = RND * 12 - 6
    GOSUB ECfunc
    col = _RGB(128 + xr * 16, 128 - xi * 16, 128 + yi * 16)
    GOSUB plot
    yr = -yr: yi = -yi 'negative root
    col = _RGB(128 + xr * 16, 128 - xi * 16, 128 + yi * 16)
    GOSUB plot
  NEXT

  _DISPLAY
  _LIMIT 15
WEND
SYSTEM

ECfunc: 'inputs (xr,xi) outputs (yr,yi)
'y^2 = q = x^3 +a*x + b 'elliptic curve over complex x
qr = xr * (xr * xr - 3 * xi * xi + a) + b
qi = xi * (3 * xr * xr - xi * xi + a)
modq = SQR(qr * qr + qi * qi)
yr = 0.707107 * SQR(modq + qr): yi = qi / (2 * yr)
RETURN

plot: 'projects point to R2 and plots it
px = cx + 25 * xr
py = cy + 25 * yr
PSET (px, py), col
RETURN