Elliptic Curves Mesh

'Elliptic Curves Mesh (c)2021 @_sorceress

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

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

WHILE INKEY$ <> CHR$(27)
  a = 10 * SIN(TIMER * .43): b = 10 * SIN(TIMER * .23)
  'a = -7: b = 10
  t = TIMER / 4: sint = SIN(t): cost = COS(t)

  CLS
  'mesh plot
  FOR xi = -5 TO 5.1 STEP .2
    FOR xr = -5 TO 5 STEP .01
      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
  NEXT
  FOR xr = -5 TO 5.1 STEP .2
    FOR xi = -5 TO 5 STEP .01
      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
  NEXT

  'real subspace
  col = _RGB(255, 255, 0): xi = 0
  FOR xr = -10 TO 10 STEP .001
    GOSUB ECfunc
    IF ABS(yi) < ABS(yr) THEN GOSUB plot: yr = -yr: yi = -yi: GOSUB plot
  NEXT

  'axes
  xr = 10: xi = 0: yr = 0: yi = 0
  GOSUB plot: LINE (cx, cy)-(px, py), _RGB(255, 0, 0): LOCATE py \ 16 + 1, px \ 8 + 1: PRINT "Re(z)"
  xr = 0: xi = 10: yr = 0: yi = 0
  GOSUB plot: LINE (cx, cy)-(px, py), _RGB(0, 255, 0): LOCATE py \ 16 + 1, px \ 8 + 1: PRINT "Im(z)"
  xr = 0: xi = 0: yr = 10: yi = 0
  GOSUB plot: LINE (cx, cy)-(px, py), _RGB(0, 0, 255): LOCATE py \ 16 + 1, px \ 8 + 1: PRINT "Re(w)"
  xr = 0: xi = 0: yr = 0: yi = 10
  GOSUB plot: LINE (cx, cy)-(px, py), _RGB(127, 127, 127): LOCATE py \ 16 + 1, px \ 8 + 1: PRINT "Im(w)"

  _DISPLAY
  _LIMIT 10
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)
IF yr <> 0 THEN yi = qi / (2 * yr) ELSE yi = SQR(-qr)
RETURN

plot: 'projects point to R2 and plots it
gx = xi * cost - yi * sint
gy = xi * sint + yi * cost
px = 400 + 25 * (xr + gx * .5)
py = 300 - 25 * (yr + gy * .5)
PSET (px, py), col
RETURN