org 100h

ADDSHIFT equ 4

  aas           ; db 3Fh: float[si-3] = 0.5
PSCALE:         ;=-16256
  db 0x80,0xC0,13h  ; long "add al,13h"
  int  10h

  add  dh,al    ; cs = ss = constants, stack, scratch
  mov  fs,dx    ; teal screen
  add  dh,al
  mov  gs,dx    ; orange screen
  add  dh,al
  mov  es,dx    ; 32768 stored points (TODO: bytes are enough when normalized correctly)
  mov  bp,SHAPE+1 ; u16[bp]=shape, [bp+si]=scratch, [bp+si+4]=MIN, [bp+si+8]=MAX

  mov  dx,3c8h  ; palette: ttttoooo (teal * orange)
  salc          ; al=0
  out  dx,al
  inc  dx
P:or   ax,1100001111000011b
  push ax
  out  dx,al    ; R
  add  al,ah
  shr  al,1
  out  dx,al    ; G
  add  al,ah
  shr  al,2
  out  dx,al    ; B
  pop  ax
  inc  ax
  jnz  P

M2:
  fninit
  fldz
  fldz
  fldz          ; {c b a}
  dec  bp
;  fst    qword[bp+si+4]  ; zero MAX and MAX2(=-MIN): needed for normalization
  pop  bx
  xor  di,di

M:
  push di       ;=0 ; sp counts frames
  call I

; combine orange and blue screens and fade them out, store to vram
; di=0, cx=0, bx=al=junk

  push es

S:xor  ax,ax
  xchg al,[gs:di]
  xchg ah,[fs:di]
  shr  ah,4
  and  al,0xF0
ROTSPEED:       ;=0x68
  push 0xA000
WIDTH: equ $-1  ;=160
  add  al,ah
  pop  es
  stosb
  loop S

;S:mov  dl,0xF0
;  mov  al,[gs:di]
;  shr  byte[gs:di],1
;  and  dl,[fs:di]
;  shr  byte[fs:di],1
;  shr  al,4
;ROTSPEED:       ;=0x68
;  push 0xA000
;WIDTH: equ $-1  ;=160
;  add  al,dl
;  pop  es
;  stosb
;  loop S

  pop  es


; draw all 32768 points with perspective
; di=0(index into points)
  mov  [bp+si],sp
  fild   word[bp+si]
  fidiv  word[byte si-100h+ROTSPEED] ; {t/ROTSPEED . . .}
  fsincos              ; {s c . . .}
  fdiv   st1,st0       ; {s c/s . . .}
D:
;rotate
  fld    st0           ; {s s c/s . . .}
  fimul  word[es:di]   ; {s*(x=p[i]) s c/s . . .}
  fld    st0           ; {sx sx s c/s . . .}
  fmul   st3           ; {cx sx s c/s . . .}
  scasw
  fld    st2           ; {s cx sx s c/s . . .}
  fimul  word[es:di]   ; {s*(z=p[i]) cx sx s c/s . . .}
  fadd   st1,st0       ; {sz sz+cx sx s c/s . . .}
  fmul   st4           ; {cz sz+cx sx s c/s . . .}
  fsubrp st2,st0       ; {Z=sz+cx X=cz-sx s c/s . . .}
  scasw
  mov  al,[es:di+3]    ; (color=p[i+6]>>8): 0..255


;rotate
;  fild   word[es:di]   ; {x=p[i] s c/s . . .}
;  fmul   st1           ; {sx s c/s . . .}
;  fld    st0           ; {sx sx s c/s . . .}
;  fmul   st3           ; {cx sx s c/s . . .}
;  scasw
;  fild   word[es:di]   ; {z=p[i+2] cx sx s c/s . . .}
;  fmul   st3           ; {sz cx sx s c/s . . .}
;  fadd   st1,st0       ; {sz sz+cx sx s c/s . . .}
;  fmul   st4           ; {cz sz+cx sx s c/s . . .}
;  fsubrp st2,st0       ; {Z=sz+cx X=cz-sx s c/s . . .}
;  scasw
;  mov  al,[es:di+3]       ; (color=p[i+6]>>8): 0..255


;don't rotate
;  fild   word[di]
;  scasw
;  fild   word[di]
;  scasw
;  fild   word[di]


  fiadd  word[si]   ; {BIG+Z X s c . . .}
  fiadd  word[si]   ; {2*BIG+Z X s c . . .}
  fidivr word[byte si-100h+WIDTH]; {w=WIDTH/(Z+2*BIG) X s c . . .}

  fmul   st1,st0       ; {w X*w s c . . .}
  fimul  word[es:di]   ; {(Y=p[i+4])*w X*w s c/s . . .}
  fistp  word[bp+si]   ; Y*w {X*w s c . . .}
  imul bx,[bp+si],320  ; handle Y overflow
  fistp  word[bp+si]   ; X*w {s c . . .}, ignore X overflow
  jo   O
  add  bx,[bp+si]
  add  bx,100*320+160  ; bx = adr = (iy+100)*320 + (ix+160)

  add  al,128
  shr  al,ADDSHIFT     ; 0..31

C:add  [fs:bx],al
  jnc  E
  sub  [fs:bx],al
E:
  xor  al,255>>ADDSHIFT

  add  [gs:bx],al
  jnc  F
  sub  [gs:bx],al
F:

O:dec  di
  dec  di
  jnz  D               ; i++

  fcompp               ; s c/s {. . .}  get rid of rotation coefficients

; esc check
  in   al,60h
  dec  al
  jnz  M

; generate about 1000 new points
; bx=constant
; di=0
; c b a = a*sin(p0*c+p1)+sin(p2*b+p3) c b
I:mov  ch,1024/256+1
  imul bx,sp,1024+128 ; where to start: sp decreases by 2
  jz   M2  ; switch to next fractal after making a rotation

R:pusha

  mov  cx,2          ; {c b a}
H:fld    st1         ; {b c b a}              ; {c sin(p0*b+p1) c b a}
  fimul  word[bp+di] ; {p0*b c b a}           ; {p2*c sin(p0*b+p1) c b a}
  ror  word[bp+di],4
  fiadd  word[bp+di] ; {p0*b+p1 c b a}        ; {p2*c+p3 sin(p0*b+p1) c b a}
  ror  word[bp+di],4
  fidiv  word[byte si-100h+PSCALE]
  fsin               ; {sin(p0*b+p1) c b a}   ; {sin(p2*c+p3) sin(p0*b+p1) c b a}
  loop H

  fmul   st4         ; {a*sin(p2*c+p3) sin(p0*b+p1) c b a}
  faddp  st1,st0     ; {a*sin(p2*c+p3)+sin(p0*b+p1) c b a}
  ffree  st3         ; {d=a*sin(p2*c+p3)+sin(p0*b+p1) c b}  ; can preserve rest with two fxch

; [B]=MAX=max(d_i), [B+4]=MAX2=max(-d_i) = -min(d_i)
; branchless min: mid = (a+b)/2; radius = abs(a-b)/2; max = mid+radius

;;  fld    st0
;;  fchs               ; {-d d . .}
;;  mov  cl,2
;;L:add  bp,4
;;  fld    st1         ; {d -d d . .}                  ; {-d MAX -d d . .}
;;  fadd   dword[bp+si]; {d+MAX -d d . .}              ; {-d+MAX2 MAX -d d . .}
;;  fld    st2         ; {d d+MAX -d d . .}            ; {-d -d+MAX2 MAX -d d . .}
;;  fsub   dword[bp+si]; {d-MAX d+MAX -d d . .}        ; {-d-MAX2 -d+MAX2 MAX -d d . .}
;;  fabs               ; {abs(d-MAX) d+MAX -d d . .}   ; {abs(-d-MAX2) -d+MAX2 MAX -d d . .}
;;  faddp  st1,st0     ; {2*max(d,MAX) -d d . .}       ; {2*max(-d,MAX2) MAX -d d . .}
;;  fmul   dword[si-3]; {MAX=max(d,MAX) -d d . .}     ; {MAX2=max(-d,MAX2) MAX -d d . .}
;;  fst    dword[bp+si]; {MAX -d d . .}                ; {-MIN=MAX2 MAX -d d . .}
;;  loop L
;;
;;; normalize
;;  fadd   st1,st0     ; {-MIN MAX-MIN -d . . .}
;;  fsubrp st2,st0     ; {MAX-MIN d-MIN . . .}
;;  fdivp  st1,st0     ; {d(0..1)=(d-MIN)/(MAX-MIN) . . .}
;;  fadd   st0         ; {d(0..2) . . .}
;;  fimul  word[si] ; {d(0..2*BIG) . . .}
;;  fisub  word[si] ; {d(-BIG..BIG) . . .}
;;

; no normalization
  fld    st0           ; 0..1 -> -6000..6000
  fadd   st0
  fimul  word[si+12]
  fisub  word[si+12]
;;


  fistp  word[es:bx]
  popa

  inc  bx
  inc  bx
  loop R

  ret

;       db 0xfc, 0xbb, 0x77, 0xfc
;       db 0x9c, 0xfe, 0x65, 0xBD
;       db 0x58, 0x7C, 0xA1, 0x78
;       db 0x06, 0x5C, 0x8B, 0xEB
       db 0x7C, 0x08, 0x80, 0x5F
       db 0xC6
SHAPE: equ $-2