org 0x7c00
        xor ax,ax
        xor bx,bx
        mov cx,0xff
        mov si,0x100
        mov di,0xfffe
        mov bp,0x091c

  mov al,13h
  int 10h
  push 0xA000
  pop es
  add dh,al
  mov gs,dx  ; z-buffer

PAL:
  mov dx,3C8h
  mov ax,cx
  out dx,al
  inc dx
  out dx,al
  shr al,1
  out dx,al
  shr al,1
  out dx,al
  loop PAL

M:
; clear/copy backbuffer
  xor di,di
CLR:
  salc
  xchg al,byte[gs:di]
  or al,al
  jnz CLNZ
  imul ax,di,-1  ; no tree: gradient
  sar ax,9
  or al,0x40
CLNZ:
  stosb
  or di,di
  jnz CLR

  mov sp,es

  in al,40h     ; pseudorandom seed
  or al,1

; generate patricles: struct Particle { u16 X,Y,Z,Active; }
G:imul ax,85    ; .X,.Y,.Z,.Active = random(-0.5..0)
  or ah,80h
  push ax
  or sp,sp
  js G

; loop until there's one particle left (on average)
  mov cl,255
L:push cx

; for each particle
;  mov si,N_PARTICLES*8 + 0x8000
  mov si,es
I:
  test word[si+6],sp ; active?
  jz EI

 ; find the nearest neighbor
  mov ax,65535  ; ax = best distance; bx = I
  mov bx,si
;  mov di,N_PARTICLES*8 + 0x8000
  mov di,es
J:
  cmp di,si ; i!=j?
  je EJ
  test word[di+6],sp ; active?
  jz EJ

  ; compute squared distance dX^2 + dY^2 + dZ^2
  push di
  push si
  push ax
  xor bp,bp  ; sum
  mov cl,3
DJ:
  lodsw
  sub ax,[di]
  
  imul ax    ; dx = dX^2 | dY^2 | dZ^2 (14-bit)
  add bp,dx  ; sum <= 48k
  scasw  ; di+=2
  loop DJ
  pop ax
  pop si
  pop di

  ; new best?
  cmp ax,bp
  jb EJ
  xchg ax,bp
  mov bx,di
EJ:
  sub di,8
  js J

 ; bx = nearest_neighbor
 ; kill this particle if the neighbor is very close
  mov [si+6],ax

 ; move this particle towards the nearest neighbor
  mov cl,3
MV:
  lodsw
  sub ax,[bx]

ISAR: ; bound the delta (poor man's normalization)
  sar ax,1
  imul dx,ax,-128
  jo ISAR

  imul bp,85 ; randomize
  mov dx,bp
  add al,dh

  sub [si-2],ax
  inc bx
  inc bx
  loop MV
  sub si,6

  sub word[si+2],-80  ; .Y += gravity

 ; draw this particle
 ; al=color ah=depth di=size bx=X bp=Y: [0 0] = center of screen
  mov ax,[si+4]  ; Z -> 0..32
  shr ax,11
  pop cx
  push cx
  cmp cl,0xD0
  jb PI
  jp PI
  add al,0xA0 ; bloom!
PI:
  call DRAW

EI:
  sub si,8
  js I
  pop cx

  in al,60h
  cmp al,1
  loopnz L
  jnz M
  int 20h

DRAW:  ; only si needs to be preserved
  mov bx,[si]  ; X
  sar bx,7
  movsx bp,byte[si+3]  ; Y

  mov di,9
  shr cx,5
  sub di,cx      ; size ~ iter

; al=color ah=depth di=size bx=X bp=Y
; draw square
  mov dx,di
PY:
  mov cx,di
PX:
  pusha
  add bp,dx
  add bx,cx
  imul di,bp,320
  lea di,[di+bx+128*320-30]
  cmp [gs:di],al
  jnb PE
  mov [gs:di],al
PE:
  popa

  loop PX
  dec dx
  jnz PY
  ret

        times 510-($-$$) db 0
        db 0x55,0xaa